Detecting relevant changes in the spatiotemporal mean function

Dette, Holger; Quanz, Pascal Article — Published Version Detecting relevant changes in the spatiotemporal mean function Journal of Time Series Analysis Provided in Cooperation with: John Wiley & Sons Suggested Citation: Dette, Holger; Quanz, Pascal (2023) : Detecting relevant changes in the spatiotemporal mean function, Journal of Time Series Analysis, ISSN 1467-9892, John Wiley & Sons, Ltd, Oxford, UK, Vol. 44, Iss. 5-6, pp. 505-532, https://doi.org/10.1111/jtsa.12674 This Version is available at: https://hdl.handle.net/10419/288140 Standard-Nutzungsbedingungen: Die Dokumente auf EconStor dürfen zu eigenen wissenschaftlichen Zwecken und zum Privatgebrauch gespeichert und kopiert werden. Sie dürfen die Dokumente nicht für öffentliche oder kommerzielle Zwecke vervielfältigen, öffentlich ausstellen, öffentlich zugänglich machen, vertreiben oder anderweitig nutzen. Sofern die Verfasser die Dokumente unter Open-Content-Lizenzen (insbesondere CC-Lizenzen) zur Verfügung gestellt haben sollten, gelten abweichend von diesen Nutzungsbedingungen die in der dort genannten Lizenz gewährten Nutzungsrechte. Terms of use: Documents in EconStor may be saved and copied for your personal and scholarly purposes. You are not to copy documents for public or commercial purposes, to exhibit the documents publicly, to make them publicly available on the internet, or to distribute or otherwise use the documents in public. If the documents have been made available under an Open Content Licence (especially Creative Commons Licences), you may exercise further usage rights as specified in the indicated licence. http://creativecommons.org/licenses/by-nc-nd/4.0/ JOURNAL OF TIME SERIES ANALYSIS J. Time Ser. Anal. 44: 505–532 (2023) Published online 11 January 2023 in Wiley Online Library (wileyonlinelibrary.com) DOI: 10.1111/jtsa.12674 SPECIAL ISSUE ARTICLE DETECTING RELEVANT CHANGES IN THE SPATIOTEMPORAL MEAN FUNCTION HOLGER DETTE AND PASCAL QUANZ Fakultät für Mathematik, Ruhr-Universität Bochum, Bochum, Germany For a spatiotemporal process {Xj(s,t)∣s∈S,t∈T}j=1,…,n, where Sdenotes the set of spatial locations and Tthe time domain, we consider the problem of testing for a change in the sequence of mean functions {𝜇j(s,t)∣s∈S,t∈T}j=1,…,n.In contrast to most of the literature, we are not interested in arbitrarily small changes but only in changes with a norm exceeding a given threshold. Asymptotically distribution free tests are proposed, which do not require the estimation of the long-run spatiotemporal covariance structure. In particular, we consider a fully functional approach and a test based on the cumulative sum paradigm, investigate the large sample properties of the corresponding test statistics and study their finite sample properties by means of simulation study. Received 16 March 2022; Accepted 21 November 2022 Keywords: Spatiotemporal process; functional data analysis; change point analysis; self-normalization; relevant hypotheses. JEL. 62M10; 62R10. 1. INTRODUCTION In many applications such as in the analysis of weatheror pollution-related data, measurements are obtained at different spatial locations over a certain time period at a high temporal frequency. Often there exists a natural segmentation of the time series such that it is reasonable to model at each spatial coordinate, say s, and on each segment, say j, the resulting data as a function, say t→Xj(s,t)of the time (on the corresponding segment). Typical examples are measurements at different geographical locations. For example, within the United States Climate Reference Network (USCRN) high resolution infrared surface temperature measurements at 126 stations in the US are publicly available on the website of the NOAA U.S. government agency. Here at each location sand each day jone observes the daily temperature curve t→Xj(s,t)(Diamond et al.,2013). Other examples include yearly curves at different locations over different years such as the daily mean temperature records from 1916 to 2018 in 40 representative Canadian cities, which are publicly available from the government of Canada website. In these applications data is typically modeled in the form Xj(s,t),s∈S,t∈T,j=1,…,n,(1.1) where Sis an arbitrary finite set and Tis a dense set (we will later consider an interval). A typical question in this context is if the sequence of mean functions {𝜇j(s,t)∣s∈S,t∈T}j=1,…,n,ofa spatiotemporal process {Xj(s,t)∣s∈S,t∈T}j=1,…,nhas changed over a specific time period. For a fixed location this corresponds to the meanwhile classical change point problem in functional data analysis (see, e.g.Berkes et al.,2009; Zhang et al.,2011; Aston and Kirch, 2012; Horváth and Kokoszka, 2012;Aueet al.,2018; Dette et al.,2020b, among many others). On the other hand, in the spatiotemporal context as considered in model (1.1) ∗Correspondence to: Holger Dette, Fakultät für Mathematik, Ruhr-Universität Bochum, Universitätsstraße 150, 44780 Bochum, Germany. E-mail: [email protected] © 2022 The Authors. Journal of Time Series Analysis published by John Wiley & Sons Ltd. This is an open access article under the terms of the Creative Commons Attribution-NonCommercial-NoDerivs License, which permits use and distribution in any medium, provided the original work is properly cited, the use is non-commercial and no modifications or adaptations are made. 506 H. DETTE AND P. QUANZ the change point problem is not well studied. Recently, Gromenko et al. (2017) proposed a test for the hypothesis of the existence of a change point in the mean function, say 𝜇j(s,t)=E[Xj(s,t)], in a sequence of independent observations. They formulated the null hypothesis and alternative in the form H0∶𝜇1=𝜇2=···=𝜇n and H1∶𝜇∶= 𝜇1=···=𝜇⌊n𝜗0⌋≠𝜇⌊n𝜗0⌋+1=···=𝜇n=∶ 𝜇+𝛿 for some 𝜗0∈(0,1)and combined the CUSUM principle with classical principal component analysis to construct a test for these hypotheses, which generalizes the approach of Berkes et al. (2009) to the spatiotemporal model (1.1). We also refer to the recent article of Zhao et al. (2021) who proposed change point analysis based on a composite likelihood criterion for a different spatiotemporal model. In contrast to this literature (and also to most of the literature on change point analysis for functional data), this article takes a different look at the change point problem. Our work is motivated by the observation that in many applications one might not be interested in arbitrarily ‘small’ changes in the mean function (in fact, one often does not believe that this function is completely constant over the whole time period for all locations). As an alternative, we therefore propose to test the hypothesis of the existence of a time point ⌊n𝜗0⌋such that the difference, say 𝛿, between the mean functions before and after this point in time is relevant. For this purpose, we define two measures of relevance. The first one corresponds to the fully functional approach as advocated in Aue et al. (2018)andis based on a norm of the difference function 𝛿. The second one is related to the FPCA approach as considered in Berkes et al. (2009) and Gromenko et al. (2017) and uses the norm of the projection of the difference on the leading principal components. The null hypothesis is then stated in the form that the squared norm is less or equal than a given threshold Δ>0, that is H0∶||𝛿||2≤Δ(see Section 2for details). We derive pivotal tests for both testing problems with neither requiring estimation of the long-run variance of the process {Xj}j=1,…,nnor the estimation of the covariance structure of the random field {Xj(s,t)∣s∈S,t∈T}. In Section 2we introduce the basic terminology and carefully define the two types of hypotheses considered in this article. Section 3is devoted to the fully functional approach, while the problem of testing relevant hypotheses by projections on the functional principal components is investigated in Section 4. Finally, in Section 5we illustrate our approach by means of a small simulation study and by the analysis of a data example. 2. RELEVANT CHANGES IN THE SPATIOTEMPORAL MEAN For an arbitrary finite set Slet L2(S×T)denote the set of all square integrable functions of the form f∶S×T→R with the common inner product ⟨f,g⟩=∑ s∈S∫T f(s,t)g(s,t)dt and the corresponding norm ||f||=⟨f,f⟩1∕2. From now on we put T=[0,1].Let{Xj}j∈Zbe a sequence of square integrable random functions on S×[0,1],where Xj=𝜇j+𝜂j,j∈Z,(2.1) {𝜂j}j∈Zis a centered error process and {𝜇j}j∈Zis a sequence of mean functions in L2(S×[0,1]). We assume that the mean functions are of the form 𝜇=𝜇1=···=𝜇⌊n𝜗0⌋,𝜇 ⌊n𝜗0⌋+1=···=𝜇n=𝜇+𝛿, wileyonlinelibrary.com/journal/jtsa © 2022 The Authors. J. Time Ser. Anal. 44: 505–532 (2023) Journal of Time Series Analysis published by John Wiley & Sons Ltd. DOI: 10.1111/jtsa.12674 CHANGE POINTS IN SPATIOTEMPORAL DATA 507 where 𝜇,𝛿 denote arbitrary deterministic but unknown elements in L2(S×[0,1])and 𝜗0∈(0,1)is a potential (unknown) change point. The case 𝛿≡0 corresponds to the situation of no change point. As explained in Section1, we are not interested in ‘small’ deviations before and after a potential change point and therefore consider the problem of monitoring the sequence for a relevant change in the mean function by testing the relevant hypotheses H0∶||𝛿||2≤Δversus H1∶||𝛿||2>Δ.(2.2) Here Δ>0 is a predefined threshold, which defines the difference before and after the time point ⌊n𝜗0⌋as relevant. Note that the case Δ=0 corresponds to the classical hypotheses (see Gromenko et al.,2017), but this case is not considered here. Our interest in hypotheses of the from (2.2) with Δ>0 stems from the fact that in applications it is often questionable to look for arbitrary small deviations. Instead, it is more reasonable to focus on (scientifically) relevant deviations, which are here defined by the threshold Δin (2.2). The choice of this threshold depends sensitively on the specific application (see Remark 3.4 for some discussion and Dette and Wied, 2014,for an example in the context of portfolio analysis based on multivariate data). We also note that for hypotheses of the form (2.2) the choice of the norm matters, as objects might be identified as close with respect to one norm (such as the L2-norm), while they might be considered as different with respect to another norm (such as the sup-norm). Moreover, we also mention that the null hypothesis and alternative in (2.2) can easily be changed, that is H0∶||𝛿||2>Δversus H1∶||𝛿||2≤Δ.(2.3) This formulation is attractive because it allows to decide for a non-relevant change (such that one can continue working under the assumption of a nearly constant mean function) at a controlled type I error. For real valued data, hypotheses of the form (2.2)and(2.3) have found considerable attention in the literature (see, e.g.themonographs of Chow and Liu, 1992; Wellek, 2010). This concept has also been used by (Liu et al.,2009; Gsteiger et al.,2011) and Dette et al. (2018) to establish the similarity of different parametric regression curves which are estimated from real valued data. In the context of functional data analysis, relevant hypotheses have been considered by (Fogarty and Small, 2014; Dette et al.,2020b) and Dette et al. (2020c) among others. A pivotal test for the hypotheses (2.2) (and as a consequence also for the hypotheses (2.3)) will be developed in Section 3. Recently, Gromenko et al. (2017) considered a different quantity to measure deviations of the difference 𝛿from the function 𝛿≡0, which is closely related to functional principal component analysis. More precisely, assume that ∶= {b1,b2,…}is a basis in L2(S×[0,1]) such that the linear span of is dense in L2(S×[0,1]).Then these authors proposed to test for a fixed order d∈Nwhether the sum of the squared scores ∑d k=1⟨𝛿,bk⟩2vanishes. In the context of testing relevant hypotheses, we are therefore interested in testing hypotheses of the form H0∶ d ∑ k=1⟨𝛿,bk⟩2≤Δversus H1∶ d ∑ k=1⟨𝛿,bk⟩2>Δ. A pivotal test for these hypotheses, where the basis functions are given by the eigenfunctions of a convex combination of the covariance kernels before and after the change point, will be developed in Section 4. We conclude this section by presenting several assumptions, which are required to prove the results in the next and the following sections. Assumption 2.1. (A1) The process {Xj}j∈Zin model (2.1) satisfies Xj={𝜇+𝜂(1) jj≤n𝜗0, 𝜇+𝛿+𝜂(2) j,j>n𝜗0, where {𝜂(1) j}j∈Zand {𝜂(2) j}j∈Zare stationary processes in L2(S×[0,1]). J. Time Ser. Anal. 44: 505–532 (2023) © 2022 The Authors. wileyonlinelibrary.com/journal/jtsa DOI: 10.1111/jtsa.12674 Journal of Time Series Analysis published by John Wiley & Sons Ltd. 508 H. DETTE AND P. QUANZ (A2) {𝜂(1) j}j∈Zand {𝜂(2) j}j∈Zform sequences of Bernoulli-shifts, that is, there exist a measurable space , measurable functions f1,f2∶∞→L2(S×[0,1])and a sequence of i.i.d., -valued and jointly (in (s,t,𝜔)) measurable random functions {𝜀j}j∈Z={𝜀j(s,t,𝜔)}j∈Zsuch that 𝜂(𝓁) j=f𝓁(𝜀j,𝜀j−1,…),𝓁=1,2 for all j∈Z. (A3) There exists a constant 𝜓∈(0,1)such that E‖‖‖𝜂(𝓁) j‖‖‖ 2+𝜓<∞,for𝓁=1,2. (A4) The sequences {𝜂(1) j}j∈Zand {𝜂(2) j}j∈Zcan be approximated by m-dependent sequences {𝜂(1) j,m}j∈Zand {𝜂(2) j,m}j∈Z, respectively, in the sense that for some 𝜅>2+𝜓 ∞ ∑ m=1(E‖‖‖𝜂(𝓁) 0−𝜂(𝓁) 0,m‖‖‖ 2+𝜓)1∕𝜅 <∞,𝓁=1,2, where 𝜂(𝓁) j,mis defined by 𝜂(𝓁) j,m=f𝓁(𝜀j,…,𝜀j−m+1,𝜺∗ j,m),𝓁=1,2 (2.4) with 𝜺∗ j,m=(𝜀∗ j,m,j−m,𝜀∗ j,m,j−m−1,…)and 𝜀∗ j,m,kare i.i.d. copies of 𝜀0and independent of {𝜀j}j∈Z. Assumption (A1) defines two segments before and after the change point. On each segment the corresponding spatiotemporal time series is assumed to be stationary, but there might be differences (in the mean and the error process) before and after the change point. The stochastic processes {𝜂(1) j}j∈Zand {𝜂(2) j}j∈Zin Assumption (A1) are not independent as they are defined by the same sequence of innovations {𝜀j}j∈Z(see Assumption (A2)). The dependence concept specified by condition (A4) has been used by (Aue et al.,2009; Hörmann and Kokoszka, 2010; Berkes et al.,2013; Horváth et al.,2014) among others. We also emphasize that our assumptions are different from those in Gromenko et al. (2017), who considered an independent and identically distributed error process {𝜂j}j∈Z. In particular, we allow for different long-run variances before and after the change point ⌊n𝜗0⌋. Moreover, these authors postulate separability in the spatiotemporal variance structure (in our case a long-run variance), which means that it factors into a purely spatial and a purely temporal component. This assumption simplifies the definition and the asymptotic analysis of their test statistics substantially. We will demonstrate below that, by using the concept of self-normalization, we can construct (asymptotically) pivotal test statistics for relevant hypotheses without any of these assumptions. 3. FULLY FUNCTIONAL DETECTION OF RELEVANT CHANGE POINTS We first consider a fully functional approach for testing the relevant hypotheses in (2.2). As in the case of the classical hypothesis H0∶||𝛿||=0, it is based on the CUSUM statistic, but it turns out that for relevant hypotheses it will be more difficult to obtain asymptotic quantiles of a corresponding test statistic. To be precise, we consider the common estimator for the unknown change point 𝜗0(see Hariz et al.,2007; Jandhyala et al.,2013, among many others) defined by  𝜗n∶= 1 narg max ⌊n𝜀⌋+1≤k≤n−⌊n𝜀⌋ k(n−k) n2‖‖‖‖‖‖ 1 k k ∑ j=1 Xj−1 n−k n ∑ j=k+1 Xj‖‖‖‖‖‖ 2 ,(3.1) wileyonlinelibrary.com/journal/jtsa © 2022 The Authors. J. Time Ser. Anal. 44: 505–532 (2023) Journal of Time Series Analysis published by John Wiley & Sons Ltd. DOI: 10.1111/jtsa.12674 CHANGE POINTS IN SPATIOTEMPORAL DATA 509 where 𝜀∈(0,1∕2)is a predefined constant. It can be shown by similar arguments as in Proposition 3.1 of Dette et al. (2020c) that, under Assumption 2.1, the estimator  𝜗nis consistent whenever ‖𝛿‖2>0and𝜗0∈(𝜀,1−𝜀), that is  𝜗n=𝜗0+oP(n−1∕2)(3.2) as n→∞. Next, for 𝜆∈[0,1],𝜗∈(𝜀,1−𝜀)we define the quantity Dn(𝜆,𝜗)∶= 1 ⌊n𝜗⌋⌊𝜆⌊n𝜗⌋⌋ ∑ j=1 Xj−1 n−⌊n𝜗⌋⌊n𝜗⌋+⌊𝜆(n−⌊n𝜗⌋)⌋ ∑ j=⌊n𝜗⌋+1 Xj∈L2(S×[0,1]),(3.3) where we also use the notation Dn(s,t,𝜆,𝜗)simultaneously to make the dependence on the spatial and temporal coordinate explicit. Note that Dnis well defined if 𝜀≥1 nand we will assume throughout this article that nis sufficiently large such that this condition is satisfied. If 𝜗=k∕nand 𝜆=1 the quantity Dn(1,k∕n)coincides with the expression in the squared norm in (3.1). Therefore, Dn(1, 𝜗n)is a natural estimator of the function 𝛿,which defines the difference before and after the change point. Consequently, it is reasonable to reject the null hypothesis in (2.2) for large values of the statistic ||Dn(1, 𝜗n)||. It will be shown later that √n(||Dn(1, 𝜗n)||2−||𝛿||2)converges weakly to a normal distribution with a complicated variance depending on a linear combination of the long-run variances of the processes {𝜂(1) j}j∈Zand {𝜂(2) j}j∈Z.To avoid its estimation, we will construct a pivotal statistic. Our main tool for this construction is the following result, which provides the weak convergence of the process {√n(||Dn(𝜆,  𝜗n)||2−𝜆2||𝛿||2)}𝜆∈[0,1]. For its statement we denote by K𝓁((s1,t1),(s2,t2)) ∶= ∑ h∈Z Cov(𝜂(𝓁) 0(s1,t1),𝜂(𝓁) h(s2,t2)) the long-run covariance kernel of the process {𝜂(𝓁) j}j∈Zfor 𝓁=1,2, which exists under Assumption 2.1 and by K((s1,t1),(s2,t2)) ∶= 1 𝜗0 K1((s1,t1),(s2,t2)) + 1 1−𝜗0 K2((s1,t1),(s2,t2)) (3.4) a scaled convex combination of these kernels. The following results can be proved by similar arguments as given in Dette et al. (2020b) and the details are omitted for the sake of brevity. Theorem 3.1. If Assumption 2.1 is satisfied and ||𝛿||>0, then √n{||Dn(𝜆,𝜗0)||2−𝜆2‖𝛿‖2}𝜆∈[0,1]⇝𝜏𝛿,𝜗0{𝜆B(𝜆)}𝜆∈[0,1],(3.5) √n{||Dn(𝜆,  𝜗n)||2−𝜆2‖𝛿‖2}𝜆∈[0,1]⇝𝜏𝛿,𝜗0{𝜆B(𝜆)}𝜆∈[0,1](3.6) as n→∞,where⇝denotes weak convergence in 𝓁∞([0,1]),{B(𝜆)}𝜆∈[0,1]denotes a standard Brownian motion and 𝜏2 𝛿,𝜗0∶= 4∑ s1,s2∈S∫∫ 𝛿(s1,t1)𝛿(s2,t2)K((s1,t1),(s2,t2))dt1dt2(3.7) with Kdefined in (3.4). J. Time Ser. Anal. 44: 505–532 (2023) © 2022 The Authors. wileyonlinelibrary.com/journal/jtsa DOI: 10.1111/jtsa.12674 Journal of Time Series Analysis published by John Wiley & Sons Ltd. 510 H. DETTE AND P. QUANZ This result leads to a very simple and pivotal test for the relevant hypotheses in (2.2). To be precise, we define  Dn∶= ||Dn(1, 𝜗n)||2,  Vn∶= (∫1 0(||Dn(𝜆,  𝜗n)||2−𝜆2||Dn(1, 𝜗n)||2)2d𝜈(𝜆))1∕2 ,(3.8) where  𝜗nis defined as in (3.1)and𝜈is a probability measure on the interval (0,1). We propose to reject the null hypothesis in (2.2) whenever  Dn>Δ+q1−𝛼(W) Vn,(3.9) where q1−𝛼(W)is the (1−𝛼)-quantile of distribution of the random variable W∶= B(1) (∫𝜆2(B(𝜆)−𝜆B(1))2d𝜈(𝜆))1∕2.(3.10) The following result shows that the decision rule (3.9) defines a consistent and asymptotic level 𝛼test for the hypotheses (2.3). Theorem 3.2. If Assumption 2.1 is satisfied, Δ>0and𝜗0∈(𝜀,1−𝜀),wehave lim n→∞P( Dn>Δ+q1−𝛼(W) Vn)=⎧ ⎪ ⎨ ⎪ ⎩ 0,if ‖𝛿‖2<Δ, 𝛼, if ‖𝛿‖2=Δ and 𝜏2 𝛿,𝜗0>0, 1,if ‖𝛿‖2>Δ. Remark 3.3. The distribution of the random variable Win (3.10) is symmetric. To see this, note that the numerator and denominator of Ware independent. This follows by using the L2representation B(t)=tZ0+∑k≥1Zksin(𝜋kt) 𝜋k of the Brownian motion for t∈[0,1]and comparing B(1)and B(t)−tB(1). Lastly, since B(1)is symmetric, the claim follows because (−B(1),(∫𝜆2(B(𝜆)−𝜆B(1))2d𝜈(𝜆))−1∕2) =(B(1),(∫𝜆2(B(𝜆)−𝜆B(1))2d𝜈(𝜆))−1∕2). Remark 3.4. (1) Note that the test (3.9) depends on the specification of the measure 𝜈on the interval (0,1), which has to be chosen in advance by the data analyst. However, in numerical experiments it turned out that this dependence does not have a significant influence on the rejection probabilities if the support of the measures has some distance to the boundaries 0 and 1 of this interval (see Section 5for some practical results). A heuristic explanation for this observation lies in the fact that the measure 𝜈appears in the definition of the statistic  Vn in (3.8) and in the quantiles of the random variable Win (3.10). Thus, intuitively there is a cancelation effect in the decision rule (3.9). (2) An important problem in applications is the choice of the threshold Δwhich is problem-specific. For this choice a careful discussion with experts from the field of application is recommended to understand in which difference they are really interested. Moreover, there are also several alternatives if this choice is difficult after these discussions. In particular, it follows from the proof of Theorem 3.1 that  Dn−||𝛿||2  Vn  −−−−−→W. wileyonlinelibrary.com/journal/jtsa © 2022 The Authors. J. Time Ser. Anal. 44: 505–532 (2023) Journal of Time Series Analysis published by John Wiley & Sons Ltd. DOI: 10.1111/jtsa.12674 CHANGE POINTS IN SPATIOTEMPORAL DATA 511 Consequently, an (asymptotic) (1−𝛼)confidence interval for the squared norm ||𝛿||2≥0 of the difference of the mean functions before and after the change point is given by [0, Dn+q1−𝛼(W) Vn].(3.11) Similarly, if it can be ruled out that the squared norm vanishes, a two sided interval for ||𝛿||2>0isgivenby (max{0, Dn−q1−𝛼∕2(W) Vn}, Dn+q1−𝛼∕2(W) Vn].(3.12) (3) It is also possible to test the relevant hypotheses in (2.2) for a finite number of thresholds Δ(1)<··· < Δ(L)simultaneously. In particular, acceptance of the null hypothesis with the threshold Δ(L0)also implies acceptance for all larger thresholds Δ(L0+1),…,Δ(L). Correspondingly, rejection for a Δ(L0)means rejection for all smaller thresholds. In this sense, evaluating the test for several thresholds is logically consistent for the user and it is possible to determine for fixed level 𝛼the largest threshold such that the null hypothesis is rejected. Remark 3.5. We expect that under additional assumptions the results of this article remain correct if the data is observed at discrete but sufficiently dense time points in the interval [0,1]and the methodology provided in this article is applied to continuous (with respect to t) curves obtained after smoothing the original data. To be a little more precise, assume, for example, that for each j=1,…,nand each location s∈Swe observe {Xj(s,t𝜈)∣𝜈=1,…,N},(3.13) where t0∶= 0<t1<···<tN−1<tN=1 defines a partition of the interval [0,1]such that maxN j=1|tj−tj−1|= O(1∕N)and Nconverges to infinity as n→∞.Forfixedjlet  Xj(s,t)denote the local linear estimator of the function (s,t)→Xj(s,t)from the data in (3.13). Then under additional smoothness assumptions on the error processes t→𝜂(𝓁)(s,t)and the mean functions t→𝜇(s,t)(s∈S,𝓁=1,2) and for a sufficiently dense partition of the interval [0,1]the estimation error maxs∈Smaxt∈[0,1]∣Xj(s,t)− Xj(s,t)∣is of order oP(n−1∕2)(uniformly with respect to j). Moreover, the procedures defined in this and the following section remain valid if they are applied to the smoothed data { Xj(s,t)∣s∈S,t∈[0,1]}j=1,…,n. Remark 3.6. As pointed out by a referee, there are numerous applications where there might exist more than one change point. For Euclidean data the literature in this context investigates multiple change points by considering classical hypotheses (see,e.g.Dette et al.,2020a, and the references therein), that is: one is interested in detecting arbitrarily small changes. On the other hand, the problem of testing for multiple relevant change points is challenging and has not been addressed rigorously so far. In the situation where 0 <𝜗 1<𝜗 2<···<𝜗 K<1 are unknown change points and Xj=𝜇+𝛿𝓁−1+𝜂(𝓁) jif ⌊n𝜗𝓁−1⌋+1≤j≤⌊n𝜗𝓁⌋;1≤𝓁≤K+1 (3.14) (with 𝛿0=0;𝜗0=0,𝜗K+1=1) Dette et al. (2020c) described a test for the hypotheses HL2 0∶1 K K ∑ 𝓁=1||𝛿𝓁−𝛿𝓁−1||2≤Δversus HL2 1∶1 K K ∑ 𝓁=1||𝛿𝓁−𝛿𝓁−1||2>Δ (the cumulative sum of changes is relevant) for L2([0,1])-valued data. This approach can be directly extended to the situation of spatiotemporal data as considered in this article. However, for other hypotheses J. Time Ser. Anal. 44: 505–532 (2023) © 2022 The Authors. wileyonlinelibrary.com/journal/jtsa DOI: 10.1111/jtsa.12674 Journal of Time Series Analysis published by John Wiley & Sons Ltd. 512 H. DETTE AND P. QUANZ such as Hmax 0∶K max 𝓁=1||𝛿𝓁−𝛿𝓁−1||2≤Δversus Hmax 1∶K max 𝓁=1||𝛿𝓁−𝛿𝓁−1||2>Δ (there exists at least one relevant change point) the situation is more complicated and will be a topic of future research. 4. TESTING RELEVANT HYPOTHESES BY CUMULATIVE SUMS OF SQUARED SCORES We address the problem of detecting relevant changes in the mean of a stationary functional time series by estimating scores. For functional data this approach has been successfully used by several authors in the context of testing classical hypotheses in the one-sample, two-sample and change point problem (see Benko et al.,2009; Berkes et al.,2009; Zhang and Shao, 2015, among others) and it has been generalized to spatiotemporal data by Gromenko et al. (2017). To the best authors’ knowledge tests for relevant hypotheses have not been constructed using this approach. To be precise, consider model (2.1) and note that in this scenario, it is possible that the covariance function also changes at the point ⌊n𝜗0⌋. Therefore, we denote by c(1)and c(2)the covariance kernels correspondingtothesamplesX1,…,X⌊n𝜗0⌋and X⌊n𝜗0⌋+1,…,Xnbefore and after the change point, respectively, and define c𝜗0∶= 𝜗0c(1)+(1−𝜗0)c(2)(4.1) as a convex combination of these two kernels. We denote by 𝜏1≥𝜏2≥…the ordered eigenvalues of the operator having covariance kernel c𝜗0with the corresponding orthonormal eigenfunctions w1,w2,…in L2(S×[0,1]).For afixedintegerd∈Nwe are interested in testing the (relevant) hypotheses H0∶ d ∑ k=1⟨𝛿,wk⟩2≤Δversus H1∶ d ∑ k=1⟨𝛿,wk⟩2>Δ,(4.2) where Δ>0 is a predefined threshold. We will briefly comment on the choice of din Section 5. Note that by Parseval’s identity ||𝛿||2=∑∞ k=1⟨𝛿,wk⟩2, and therefore - similar as for testing classical hypotheses – a test for the hypotheses (4.2) can also be used for the hypotheses (2.2). We refer to Remark 4.4 for a more detailed discussion of this approach in the context of testing relevant hypotheses. For the statements in this section we require the following assumptions. Assumption 4.1. The process {Xj}j∈Zin model (2.1) satisfies conditions (A1) and (A2) of Assumption 2.1. Furthermore, {Xj}j∈Zsatisfies (A3’) There exists a constant 𝜓∈(0,1)such that E‖‖‖𝜂(𝓁) j‖‖‖ 4+𝜓<∞,𝓁=1,2. (A4’) The sequences {𝜂(1) j}j∈Zand {𝜂(2) j}j∈Zcan be approximated by m-dependent sequences {𝜂(1) j,m}j∈Zand {𝜂(2) j,m}j∈Z, respectively, in the sense that for some 𝜅>4+𝜓 ∞ ∑ m=1(E‖‖‖𝜂(𝓁) 0−𝜂(𝓁) 0,m‖‖‖ 4+𝜓)1∕𝜅 <∞,𝓁=1,2, where 𝜂(𝓁) j,mis defined (2.4). wileyonlinelibrary.com/journal/jtsa © 2022 The Authors. J. Time Ser. Anal. 44: 505–532 (2023) Journal of Time Series Analysis published by John Wiley & Sons Ltd. DOI: 10.1111/jtsa.12674 CHANGE POINTS IN SPATIOTEMPORAL DATA 519 Figure 4. Empirical rejection probabilities of the test (3.9) for the hypotheses (2.2) with Δ=0.15. The error process is given by an fMA(1) model and different measures in the statistic (3.8) are considered. Left panel: 𝜈4; middle panel: 𝜈9; right panel: 𝜈19 Figure 5. Empirical rejection probabilities of the test (4.11) for the hypotheses (4.2) with Δ=0.15. The difference 𝛿between the mean functions is given by (5.1) and different error processes are considered. Left panel: scaled Brownian motion (d=15) in (5.2); middle panel: non-separable process (d=11) in (5.3); right panel: fMA(1) process (d=11) in (5.5) 5.1.2. Relevant Changes by Cumulative Sums of Squared Scores We briefly illustrate the finite sample properties of the test (4.11) for the hypotheses (4.2),whereweusethesame scenarios as before. The test requires the choice of the number of principal functional components and we choose the parameter dsuch that 95% of the variance in the data will be explained. This results in d=15, d=11 and d=11 functional principal components for the models (5.2), (5.3), and (5.5) respectively. The corresponding rejection probabilities are displayed in Figure 5and we observe the qualitative behavior predicted by Theorem 4.3. We also observe that the test (4.11) is conservative in the case of the non-separable process (5.3), while the nominal level at the boundary of the hypotheses ∑d k=1⟨𝛿,bk⟩2=Δis very well approximated for the scaled Brownian motion (5.2). Finally, we investigate the impact of the measure 𝜈in the scaling factor (4.10), where we again restrict ourselves to the case of an fMA(1) process and the uniform distributions 𝜈kin (5.7)fork=4,9 and 19. The corresponding results are shown in Figure 6and demonstrate that the test (4.11) is not very sensitive with respect to this choice, either. J. Time Ser. Anal. 44: 505–532 (2023) © 2022 The Authors. wileyonlinelibrary.com/journal/jtsa DOI: 10.1111/jtsa.12674 Journal of Time Series Analysis published by John Wiley & Sons Ltd. 520 H. DETTE AND P. QUANZ Figure 6. Empirical rejection probabilities of the test (4.11) for the hypotheses (4.2) with Δ=0.15. The error process is given by an fMA(1) model and different measures in the statistic (3.8) are considered. Left panel: 𝜈4; middle panel: 𝜈9; right panel: 𝜈19 Table I. Results of the test (3.9) for the Canadian weather data for different nominal level and different thresholds Δ∕𝛼10% 5% 1% 0.004 Reject Reject Accept 0.005 Reject Accept Accept ⋮⋮ ⋮ ⋮ 2.974 Reject Accept Accept 2.975 Accept Accept Accept Note: The value 0.004 represents the maximal threshold in (2.3) such that the null hypothesis of no relevant change is rejected at nominal level 5%, whereas the value 2.974 represents the maximal threshold in (2.3) such that the null hypothesis is rejected at nominal level 10%. 5.2. Data Example We conclude this article with an application of the two test procedures in a real data example. For this purpose, we use Canadian weather data, which consists of daily measurements at 40 representative Canadian cities. Thus, we observe yearly curves over different years (1916 −2018) at different locations. The data can be downloaded from the government of Canada website: https://climate.weather.gc.ca/historical_data/search_historic_ data_e.html. The available data contains several different measurements such as maximum/minimum temperature or precipitation amount. For the sake of brevity, we focus on the average daily temperatures. Due to missing values in the reported temperature data, four stations were chosen such that a large amount of available data overlaps and only little parts had to be interpolated or removed: Calgary International Airport, Alberta (ID: 2205), Medicine Hat Airport, Alberta (ID: 2273), Indian Head CDA, Saskatchewan (ID: 2925) and Ottawa CDA, Ontario (ID: 4333). In the notation of the previous sections this means S={1,2,3,4}and the sample size is given by n=113, which corresponds to the period from the years 1891 to 2007 without the years 1910, 1911, 1993 and 1995. The change point estimator in (3.1) yields  𝜗n=0.708, which approximately corresponds to the year 1969. We emphasize that this result depends on the definition of the segmentation and we have used the common yearly segmentation starting on January 1st and ending on December 31st (see also Besse et al.,2000; Sharipov et al.,2016; Gromenko et al.,2017, who used this segmentation as well). If one uses a different yearly segmentation the results may change slightly. For example, the yearly segmentation ranging from July 1st to June 30th yields  𝜗n=0.699 as estimate for the change point corresponding to the years 1968/69. wileyonlinelibrary.com/journal/jtsa © 2022 The Authors. J. Time Ser. Anal. 44: 505–532 (2023) Journal of Time Series Analysis published by John Wiley & Sons Ltd. DOI: 10.1111/jtsa.12674 CHANGE POINTS IN SPATIOTEMPORAL DATA 521 Figure 7. The function ein (5.8) for j=1,…,100 Table II. Results of the test (4.11) with d=51 principal components for the Canadian weather data for different nominal level and different thresholds Δ∕𝛼10% 5% 1% 0.531 Reject Reject Reject 0.532 Reject Reject Accept ⋮⋮ ⋮ ⋮ 3.301 Reject Reject Accept 3.302 Reject Accept Accept ⋮⋮ ⋮ ⋮ 4.467 Reject Accept Accept 4.468 Accept Accept Accept Note: The value 4.467 represents the maximal threshold in (4.2) such that the null hypothesis of no relevant change is rejected at nominal level 10%, the value 3.301 is the maximal threshold for the nominal level 5%and the value 0.531 is the maximal threshold for the nominal level 1%. Next we apply the testing procedures introduced in Section 3and 4starting with the fully functional approach. The results of the test (3.9) for different thresholds and different nominal levels are given in Table I. We observe that Δ=2.974 is the largest threshold such that the test (3.9) rejects the null hypothesis at nominal level 𝛼=0.1. Because there are four stations, this corresponds to an average effect of 2.974∕4≈0.74. Finally, we note that the one-sided confidence interval for ||𝛿||≥0in(3.11)isgivenby[0,4.397], while the two-sided interval for ||𝛿||>0 in (3.12) is obtained as (0.067,4.675]. Lastly, we consider the test based on cumulative sums of squared scores developed in Section 4. In this case, a choice for the value of dmust be made and we display in Figure 7the ratios j→e(j)∶= ∑j k=1 𝜆k trace (c 𝜗n),(5.8) where  𝜆1, 𝜆2,…are the estimated eigenvalues of the covariance operator corresponding to the covariance function in (4.1). We observe that the eigenvalues are slowly decreasing and we choose d=51, which results in a value of e(51)≈0.849 of explained variance. The results of the test (4.11) for the relevant hypotheses (4.2) are shown in Table II for different values of Δand 𝛼. We observe that the maximal threshold in (4.2) such that the null hypothesis of no relevant change is rejected at nominal level 10%, is given by 4.467. Finally, one and two-sided confidence J. Time Ser. Anal. 44: 505–532 (2023) © 2022 The Authors. wileyonlinelibrary.com/journal/jtsa DOI: 10.1111/jtsa.12674 Journal of Time Series Analysis published by John Wiley & Sons Ltd. 522 H. DETTE AND P. QUANZ intervals for the quantity (∑d k=1⟨𝛿,wk⟩2)1∕2are obtained in same way as described in Remark 3.4 and are given by [0,3.254]and (1.877,3.396], respectively (d=51). ACKNOWLEDGEMENTS The authors would like to thank two unknown referees for their constructive comments on an earlier version of this article. This research has been supported in part by the German Research Foundation (DFG), project number 45723897. Open Access funding enabled and organized by Projekt DEAL. DATA AVAILABILITY STATEMENT The data that support the findings of this study are available at https://climate.weather.gc.ca/historical_data/search_ historic_data_e.html. REFERENCES Aston JA, Kirch C. 2012. Detecting and estimating changes in dependent functional data. Journal of Multivariate Analysis 109:204–220. Aue A, Hörmann S, Horváth L, Reimherr M. 2009. Break detection in the covariance structure of multivariate time series models. The Annals of Statistics 37(6B):4046–4087. Aue A, Rice G, Sönmez O. 2018. Detecting and dating structural breaks in functional data without dimension reduction. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 80(3):509–529. Aue A, Dette H, Rice G. 2023. Two-sample tests for relevant differences in the eigenfunctions of covariance operators. Statistica Sinica 33:1–27. Benko M, Härdle W, Kneip A. 2009. Common functional principal components. The Annals of Statistics 37(1):1–34. Berkes I, Gabrys R, Horváth L, Kokoszka P. 2009. Detecting changes in the mean of functional observations. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 71(5):927–946. Berkes I, Horváth L, Rice G. 2013. Weak invariance principles for sums of dependent random functions. Stochastic Processes and their Applications 123(2):385–403. Besse PC, Cardot H, Stephenson DB. 2000. Autoregressive forecasting of some functional climatic variations. Scandinavian Journal of Statistics 27(4):673–687. Chow S-C, Liu P-J. 1992. Design and Analysis of Bioavailability and Bioequivalence Studies Marcel Dekker, New York. Dette H, Wied D. 2014. Detecting relevant changes in time series models. Journal of the Royal Statistical Society 78(2):371–394. Dette H, Möllenhoff K, Volgushev S, Bretz F. 2018. Equivalence of regression curves. Journal of the American Statistical Association 113:711–729. Dette H, Eckle T, Vetter M. 2020a. Multiscale change point detection for dependent data. Scandinavian Journal of Statistics 47(4):1243–1274. Dette H, Kokot K, Aue A. 2020b. Functional data analysis in the Banach space of continuous functions. Annals of Statistics 48(2):1168–1192. Dette H, Kokot K, Volgushev S. 2020c. Testing relevant hypotheses in functional time series via self-normalization. Journal of the Royal Statistical Society (B) 82(3):629–660. Diamond HJ, Karl TR, Palecki MA, Baker CB, Bell JE, Leeper RD, Easterling DR, Lawrimore JH, Meyers TP, Helfert MR, Goodge G, Thorne PW. 2013. U.S. Climate Reference Network after One Decade of Operations: Status and Assessment. Bulletin of the American Meteorological Society 94(4):485–498. Fogarty CB, Small DS. 2014. Equivalence testing for functional data with an application to comparing pulmonary function devices. Annals of Applied Statistics 8(4):2002–2026. Gromenko O, Kokoszka P, Reimherr M. 2017. Detection of change in the spatiotemporal mean function. Journal of the Royal Statistical Society (B) 79(1):29–50. Gsteiger S, Bretz F, Liu W. 2011. Simultaneous confidence bands for nonlinear regression models with application to population pharmacokinetic analyses. Journal of Biopharmaceutical Statistics 21(4):708–725. Hariz SB, Wylie JJ, Zhang Q. 2007. Optimal rate of convergence for nonparametric change-point estimators for nonstationary sequences. The Annals of Statistics 35:1802–1826. Hörmann S, Kokoszka P. 2010. Weakly dependent functional data. The Annals of Statistics 38(3):1845–1884. Horváth L, Kokoszka P. 2012. Inference for Functional Data with Applications Springer-Verlag, New York. wileyonlinelibrary.com/journal/jtsa © 2022 The Authors. J. Time Ser. Anal. 44: 505–532 (2023) Journal of Time Series Analysis published by John Wiley & Sons Ltd. DOI: 10.1111/jtsa.12674 CHANGE POINTS IN SPATIOTEMPORAL DATA 523 Horváth L, Kokoszka P, Rice G. 2014. Testing stationarity of functional time series. Journal of Econometrics 179(1):66–82. Jandhyala V, Fotopoulos S, MacNeill I, Liu P. 2013. Inference for single and multiple change-points in time series. Journal of Time Series Analysis 34(4):423–446. Liu W, Bretz F, Hayter AJ, Wynn HP. 2009. Assessing non-superiority, non-inferiority of equivalence when comparing two regression models over a restricted covariate region. Biometrics 65(4):1279–1287. Sharipov O, Tewes J, Wendler M. 2016. Sequential block bootstrap in a Hilbert space with application to change point analysis. Canadian Journal of Statistics 44(3):300–322. Wellek S. 2010. Testing Statistical Hypotheses of Equivalence and Noninferiority Chapman and Hall/CRC, New York. Zhang X, Shao X. 2015. Two sample inference for the second-order property of temporally dependent functional data. Bernoulli 21(2):909–929. Zhang X, Shao X, Hayhoe K, Wuebbles DJ. 2011. Testing the structural stability of temporally dependent functional observations and application to climate projections. Electronic Journal of Statistics 5:1765–1796. Zhao Z, Ma TF, Ng WL, Yau CY (2021). A composite likelihood-based approach for change-point detection in spatio-temporal process. arXiv:1904.06340. APPENDIX : PROOFS A.1. Some Preliminary Results We start with some preparations and present several results which are used in the proofs. We define for s∈Sand t∈[0,1] Sn(s,t,𝜆)∶=1 n⌊𝜆n⌋ ∑ j=1(Xj(s,t)−E[Xj(s,t)]) and state the following result, which can be obtained by generalizing Theorem 1.1 in Berkes et al. (2013)tothe space L2(S×[0,1]). Theorem A.1. If Assumption 2.1 is satisfied, there exists a sequence of Gaussian processes ({Γn(s,t,𝜆)∣s∈S,0≤𝜆, t≤1})n∈Nsuch that sup 0≤𝜆≤1∑ s∈S∫(√nSn(s,t,𝜆)−Γ n(s,t,𝜆))2dt=sup 0≤𝜆≤1‖‖‖√nSn(⋅,⋅,𝜆)−Γ n(⋅,⋅,𝜆)‖‖‖ 2=oP(1) and {Γn(s,t,𝜆)∣s∈S,0≤𝜆,t≤1} ={Γ(s,t,𝜆)∣s∈S,0≤𝜆,t≤1}, where Γ(s,t,𝜆)∶= ∞ ∑ i=1√𝜆i𝜙i(s,t)Wi(𝜆), 𝜆i,𝜙iare the eigenvalues and eigenvectors of the covariance operator of Xjand {Wi}i∈Nare independent standard Brownian motions. Our next auxiliary result quantifies the difference between the processes {Dn(𝜆,𝜗0)}𝜆∈[0,1]and {Dn(𝜆,  𝜗n)}𝜆∈[0,1] as n→∞. Lemma A.2. If Assumption 2.1 is satisfied, then (1) sup0≤𝜆≤1‖‖Dn(𝜆,𝜗)‖‖=OP(1)for 𝜗∈{𝜗0, 𝜗n}, (2) sup0≤𝜆≤1||Dn(𝜆,  𝜗n)−Dn(𝜆,𝜗0)||=oP(n−1∕2). J. Time Ser. Anal. 44: 505–532 (2023) © 2022 The Authors. wileyonlinelibrary.com/journal/jtsa DOI: 10.1111/jtsa.12674 Journal of Time Series Analysis published by John Wiley & Sons Ltd. 524 H. DETTE AND P. QUANZ Proof. We can assume that 1 n≤𝜆because, by definition, Dn(𝜆,𝜗0)≡Dn(𝜆,  𝜗n)≡0if0≤𝜆<1 n, and both assertions are trivially true. For a proof of part (1) we note that with 𝜅from (A4) we have sup 1 n≤𝜆≤1 1 √⌊𝜆n⌋‖‖‖‖‖‖⌊𝜆n⌋ ∑ j=1 𝜂(𝓁) j‖‖‖‖‖‖ =OP(log1∕𝜅(n)),𝓁=1,2,(A1) which follows from Lemma B.1 in the online supplementary material of Aue et al. (2018). Let 𝜆∈[1∕n,1].It suffices to show the assertion for 𝜗=𝜗0because the second part of the Lemma implies the statement for 𝜗= 𝜗n. Then ‖‖Dn(𝜆,𝜗0)‖‖≤‖‖‖ Dn(𝜆,𝜗0)‖‖‖+𝜆‖𝛿‖+oP(1)=‖‖‖ Dn(𝜆,𝜗0)‖‖‖+OP(1), where the process { Dn(𝜆,𝜗)}𝜆∈[0,1]in L2(S×[0,1])is defined by  Dn(𝜆,𝜗)∶= 1 ⌊n𝜗⌋⌊𝜆⌊n𝜗⌋⌋ ∑ j=1 𝜂(1) j−1 n−⌊n𝜗⌋⌊n𝜗⌋+⌊𝜆(n−⌊n𝜗⌋)⌋ ∑ j=⌊n𝜗⌋+1 𝜂(2) j. Note that  Dn(𝜆,𝜗)is the centered version of Dn(𝜆,𝜗)defined in (3.3)andthat Dn(𝜆,𝜗0)+𝜆𝛿 = Dn(𝜆,𝜗0)+OP(n−1)(A2) (note that 𝜗0∈(𝜀,1−𝜀)and that 𝜀≥1 nif nis sufficiently large). We have sup 1 n≤𝜆≤1‖‖ Dn(𝜆,𝜗0)‖‖≤sup 1 n≤𝜆≤1 1 ⌊n𝜗0⌋‖‖‖‖‖‖⌊𝜆⌊n𝜗0⌋⌋ ∑ j=1 𝜂(1) j‖‖‖‖‖‖ +sup 1 n≤𝜆≤1 1 n−⌊n𝜗0⌋‖‖‖‖‖‖⌊n𝜗0⌋+⌊𝜆(n−⌊n𝜗0⌋)⌋ ∑ j=⌊n𝜗0⌋+1 𝜂(2) j‖‖‖‖‖‖, where the first term can be estimated as follows sup 1 n≤𝜆≤1 1 ⌊n𝜗0⌋‖‖‖‖‖‖⌊𝜆⌊n𝜗0⌋⌋ ∑ j=1 𝜂(1) j‖‖‖‖‖‖≤1 √⌊n𝜗0⌋sup 1 n≤𝜆≤1 1 √⌊𝜆⌊n𝜗0⌋⌋‖‖‖‖‖‖⌊𝜆⌊n𝜗0⌋⌋ ∑ j=1 𝜂(1) j‖‖‖‖‖‖ ≤1 √⌊n𝜗0⌋sup 1 n≤𝜆≤1 1 √⌊𝜆n⌋‖‖‖‖‖‖⌊𝜆n⌋ ∑ j=1 𝜂(1) j‖‖‖‖‖‖ =OP(log1∕𝜅(n) √n)=oP(1). (here in the second inequality, we expanded the set over which we take the supremum and in the last step, we used the estimate (A1)). A similar argument provides the same rate for the second term, which proves part (1) of Lemma A.2. For the proof of the second assertion we need the following (slightly more general) statement from the beginning of Section B.1 in Dette et al. (2020c), which states that the sequence of processes {Γn(s,t,𝜆)∣s∈S,0≤t,𝜆≤1}n∈Nin Theorem A.1 satisfies sup 𝜈,𝜆∈[0,1]∶ |𝜈−𝜆|≤𝜅n‖‖Γn(⋅,⋅,𝜈)−Γ n(⋅,⋅,𝜆)‖‖2=oP(1).(A3) wileyonlinelibrary.com/journal/jtsa © 2022 The Authors. J. Time Ser. Anal. 44: 505–532 (2023) Journal of Time Series Analysis published by John Wiley & Sons Ltd. DOI: 10.1111/jtsa.12674 CHANGE POINTS IN SPATIOTEMPORAL DATA 525 for any and positive sequence (𝜅n)n∈Nwith 𝜅n→0. By adding and subtracting 𝜆𝛿 we have Dn(𝜆,  𝜗n)−Dn(𝜆,𝜗0)=  Dn(𝜆,  𝜗n)− Dn(𝜆,𝜗0)+OP(n−1),where  Dn(𝜆,  𝜗n)− Dn(𝜆,𝜗0)= 1 ⌊n 𝜗n⌋⌊𝜆⌊n 𝜗n⌋⌋ ∑ j=1 𝜂(1) j−1 ⌊n𝜗0⌋⌊𝜆⌊n𝜗0⌋⌋ ∑ j=1 𝜂(1) j +1 n−⌊n𝜗0⌋⌊n𝜗0⌋+⌊𝜆(n−⌊n𝜗0⌋)⌋ ∑ j=⌊n𝜗0⌋+1 𝜂(2) j−1 n−⌊n 𝜗n⌋⌊n 𝜗n⌋+⌊𝜆(n−⌊n 𝜗n⌋)⌋ ∑ j=⌊n 𝜗n⌋+1 𝜂(2) j. For the first difference we obtain by a similar argument as in the proof of part (1) that sup 1 n≤𝜆≤1√n‖‖‖‖‖‖ 1 ⌊n 𝜗n⌋⌊𝜆⌊n 𝜗n⌋⌋ ∑ j=1 𝜂(1) j−1 ⌊n𝜗0⌋⌊𝜆⌊n𝜗0⌋⌋ ∑ j=1 𝜂(1) j‖‖‖‖‖‖ ≤n ⌊n 𝜗n⌋sup 1 n≤𝜆≤1‖‖‖‖‖‖ 1 √n⌊𝜆⌊n 𝜗n⌋⌋ ∑ j=1 𝜂(1) j−Γ n(⋅,⋅,𝜆⌊n 𝜗n⌋ n)‖‖‖‖‖‖ +n ⌊n𝜗0⌋sup 1 n≤𝜆≤1‖‖‖‖‖‖ 1 √n⌊𝜆⌊n𝜗0⌋⌋ ∑ j=1 𝜂(1) j−Γ n(⋅,⋅,𝜆⌊n𝜗0⌋ n)‖‖‖‖‖‖ +n ⌊n 𝜗n⌋sup 1 n≤𝜆≤1‖‖‖‖‖‖ Γn(⋅,⋅,𝜆⌊n 𝜗n⌋ n)−Γ n(⋅,⋅,𝜆 𝜗n)‖‖‖‖‖‖ +n ⌊n𝜗0⌋sup 1 n≤𝜆≤1‖‖‖‖‖Γn(⋅,⋅,𝜆⌊n𝜗0⌋ n)−Γ n(⋅,⋅,𝜆𝜗0)‖‖‖‖‖ +sup 1 n≤𝜆≤1‖‖‖‖‖ n ⌊n 𝜗n⌋Γn(⋅,⋅,𝜆 𝜗n)− n ⌊n𝜗0⌋Γn(⋅,⋅,𝜆𝜗0)‖‖‖‖‖ ≤sup 𝜈,𝜆∈[0,1]∶ |𝜈−𝜆|≤1 n‖‖Γn(⋅,⋅,𝜈)−Γ n(⋅,⋅,𝜆)‖‖ +O(1)sup 𝜈,𝜆∈[0,1]∶ |𝜈−𝜆|≤| 𝜗n−𝜗0|‖‖Γn(⋅,⋅,𝜈)−Γ n(⋅,⋅,𝜆)‖‖+|𝜗0− 𝜗n| 𝜗0 𝜗n OP(1)+oP(1) =oP(1), wherewehaveused(A3) in the last equality. Assertion (2) of Lemma A.2 now follows by using a similar argument for the second difference. ◾ We conclude our preparations recalling the definition of the inner product in (4.4) and state a lemma regarding the weak convergence of the process (Z(1) n,Z(2) n)⊤∶= {(Z(1) n(𝜆),Z(2) n(𝜆))⊤}𝜆∈[0,1](A4) ∶= {(1 √n⌊𝜆n⌋ ∑ j=1(⟨𝜂(k) j⊗𝜂(k) j−c(k),𝜉k⟩2+⟨𝜂(k) j,𝜁k⟩))k=1,2}𝜆∈[0,1] , J. Time Ser. Anal. 44: 505–532 (2023) © 2022 The Authors. wileyonlinelibrary.com/journal/jtsa DOI: 10.1111/jtsa.12674 Journal of Time Series Analysis published by John Wiley & Sons Ltd. 526 H. DETTE AND P. QUANZ where {𝜂(1) j}j∈Zand {𝜂(2) j}j∈Zare centered processes in L2(S×[0,1])and 𝜁1,𝜁2∈L2(S×[0,1]),𝜉1,𝜉2∈ L2((S×[0,1])2)are given functions. We emphasize that we consider the process (Z(1) n,Z(2) n)⊤with different parameters 𝜁1,𝜁2∈L2(S×[0,1]),𝜉1,𝜉2∈L2((S×[0,1])2)in the proofs of the results of Section 3and 4. The proof of the following result is similar to the proof of Lemma B.1 in Dette et al. (2020c) and therefore omitted. Lemma A.3. Let 𝜁1,𝜁2∈L2(S×[0,1]),𝜉1,𝜉2∈L2((S×[0,1])2)be fixed but arbitrary functions and let {𝜂(1) j}j∈Z, {𝜂(2) j}j∈Zdenote centered processes satisfying Assumption 4.1. Then the process defined in (A4) converges weakly in 𝓁∞([0,1])2,thatis (Z(1) n,Z(2) n)⊤⇝Σ1∕2(B1,B2)⊤,(A5) where B1,B2are independent Brownian motions and Σ=(Σ kl)k,l=1,2is a 2×2 matrix with entries Σkl =∑ h∈Z(Cov(⟨𝜂(k) 0⊗𝜂(k) 0,𝜉k⟩2,⟨𝜂(l) h⊗𝜂(l) h,𝜉l⟩2)+Cov(⟨𝜂(k) 0⊗𝜂(k) 0,𝜉k⟩2,⟨𝜂(l) h,𝜁l⟩) +Cov(⟨𝜂(k) 0,𝜁k⟩,⟨𝜂(l) h⊗𝜂(l) h,𝜉l⟩2)+Cov(⟨𝜂(k) 0,𝜁k⟩,⟨𝜂(l) h,𝜁l⟩)).(A6) Moreover, in the case 𝜉1≡𝜉2≡0 Assumption 2.1 instead of Assumption 4.1 is sufficient for the weak convergence in (A5). A.2. Proof of Theorem 4.1 We will show that sup 0≤𝜆≤1√𝜆‖‖‖‖ 𝜗nc(1) 𝜆−𝜗0E[𝜂(1) 0⊗𝜂(1) 0]‖‖‖‖2=OP(log2∕𝜅(n) √n),(A7) sup 0≤𝜆≤1√𝜆‖‖‖‖(1− 𝜗n)c(2) 𝜆−(1−𝜗0)E[𝜂(2) 0⊗𝜂(2) 0]‖‖‖‖2=OP(log2∕𝜅(n) √n).(A8) For the sake of brevity, we restrict ourselves to a proof of (A7); the proof of (A8) follows by similar arguments. First, we assume that 𝜆<⌊n 𝜗n⌋−1,thenc(1) 𝜆=0and √𝜆‖‖‖‖ 𝜗nc(1) 𝜆−𝜗0E[𝜂(1) 0⊗𝜂(1) 0]‖‖‖‖2<𝜗0 √⌊n 𝜗n⌋|||E[𝜂(1) 0⊗𝜂(1) 0]|||=OP(n−1∕2) uniformly in 𝜆∈[0,⌊n 𝜗n⌋−1).If⌊n 𝜗n⌋−1≤𝜆≤1, we use the representation c(1) 𝜆=1 ⌊𝜆⌊n 𝜗n⌋⌋⌊𝜆⌊n 𝜗n⌋⌋ ∑ j=1 Xj⊗Xj−1 ⌊𝜆⌊n 𝜗n⌋⌋2⌊𝜆⌊n 𝜗n⌋⌋ ∑ i,j=1 Xi⊗Xj(A9) and consider the first term. For this we get 1 ⌊𝜆⌊n 𝜗n⌋⌋⌊𝜆⌊n 𝜗n⌋⌋ ∑ j=1 Xj⊗Xj =1 ⌊𝜆⌊n 𝜗n⌋⌋⌊𝜆⌊n 𝜗n⌋⌋∧⌊n𝜗0⌋ ∑ j=1 Xj⊗Xj+1 ⌊𝜆⌊n 𝜗n⌋⌋ ⌊𝜆⌊n 𝜗n⌋⌋ ∑ j=⌊𝜆⌊n 𝜗n⌋⌋∧⌊n𝜗0⌋+1 Xj⊗Xj wileyonlinelibrary.com/journal/jtsa © 2022 The Authors. J. Time Ser. Anal. 44: 505–532 (2023) Journal of Time Series Analysis published by John Wiley & Sons Ltd. DOI: 10.1111/jtsa.12674 CHANGE POINTS IN SPATIOTEMPORAL DATA 527 =1 ⌊𝜆⌊n 𝜗n⌋⌋⌊𝜆⌊n 𝜗n⌋⌋∧⌊n𝜗0⌋ ∑ j=1(𝜂(1) j⊗𝜂(1) j+𝜇⊗𝜂 (1) j+𝜂(1) j⊗𝜇+𝜇⊗𝜇 ) +1 ⌊𝜆⌊n 𝜗n⌋⌋ ⌊𝜆⌊n 𝜗n⌋⌋ ∑ j=⌊𝜆⌊n 𝜗n⌋⌋∧⌊n𝜗0⌋+1(𝜂(2) j⊗𝜂(2) j+𝜇⊗𝜂 (2) j+𝜂(2) j⊗𝜇+𝜇⊗𝜇),(A10) where we use the notation 𝜇 ∶= 𝜇+𝛿. The second term in (A9) can be rewritten as follows: 1 ⌊𝜆⌊n 𝜗n⌋⌋2⌊𝜆⌊n 𝜗n⌋⌋ ∑ i,j=1 Xi⊗Xj =1 ⌊𝜆⌊n 𝜗n⌋⌋2⌊𝜆⌊n 𝜗n⌋⌋∧⌊n𝜗0⌋ ∑ i,j=1(𝜂(1) i⊗𝜂(1) j+𝜇⊗𝜂 (1) j+𝜂(1) i⊗𝜇+𝜇⊗𝜇 ) +1{⌊n𝜗0⌋≤⌊𝜆⌊n 𝜗n⌋⌋} ⌊𝜆⌊n 𝜗n⌋⌋2⎧ ⎪ ⎨ ⎪ ⎩ ⌊n𝜗0⌋ ∑ i=1⌊𝜆⌊n 𝜗n⌋⌋ ∑ j=⌊n𝜗0⌋+1(𝜂(1) i⊗𝜂(2) j+𝜇⊗𝜂 (2) j+𝜂(1) i⊗𝜇+𝜇⊗𝜇) +⌊𝜆⌊n 𝜗n⌋⌋ ∑ i=⌊n𝜗0⌋+1⌊n𝜗0⌋ ∑ j=1(𝜂(2) i⊗𝜂(1) j+𝜂(2) i⊗𝜇+𝜇⊗𝜂 (1) j+𝜇⊗𝜇 ) +⌊𝜆⌊n 𝜗n⌋⌋ ∑ i,j=⌊n𝜗0⌋+1(𝜂(2) i⊗𝜂(2) j+𝜂(2) i⊗𝜇+𝜇⊗𝜂 (2) j+𝜇⊗𝜇)⎫ ⎪ ⎬ ⎪ ⎭ . We now consider the cases ⌊𝜆⌊n 𝜗n⌋⌋<⌊n𝜗0⌋and ⌊𝜆⌊n 𝜗n⌋⌋≥⌊n𝜗0⌋separately. For this purpose, we define the set Λn∶= {⌊n 𝜗n⌋−1≤𝜆≤1|⌊𝜆⌊n 𝜗n⌋⌋<⌊n𝜗0⌋}and note that sup 1 n≤𝜆≤1√𝜆‖‖‖‖ 𝜗nc(1) 𝜆−𝜗0E[𝜂(1) 0⊗𝜂(1) 0]‖‖‖‖2 =max{sup 𝜆∈Λn√𝜆‖‖‖‖ 𝜗nc(1) 𝜆−𝜗0E[𝜂(1) 0⊗𝜂(1) 0]‖‖‖‖2,sup 𝜆∈ΛC n√𝜆‖‖‖‖ 𝜗nc(1) 𝜆−𝜗0E[𝜂(1) 0⊗𝜂(1) 0]‖‖‖‖2}. We consider each case individually, starting with the first term, where the supremum is taken over the set Λn, that is, the inequality ⌊𝜆⌊n 𝜗n⌋⌋<⌊n𝜗0⌋holds. Observing that the sum in (A10) vanishes, we obtain in this case that c(1) 𝜆=1 ⌊𝜆⌊n 𝜗n⌋⌋⌊𝜆⌊n 𝜗n⌋⌋ ∑ j=1 𝜂(1) j⊗𝜂(1) j−1 ⌊𝜆⌊n 𝜗n⌋⌋2⌊𝜆⌊n 𝜗n⌋⌋ ∑ i,j=1 𝜂(1) i⊗𝜂(1) j.(A11) In the subsequent discussion, we will repeatedly make use of the following inequality, without explicitly mentioning it. It is obtained by expanding the set over which the supremum is calculated and substituting J. Time Ser. Anal. 44: 505–532 (2023) © 2022 The Authors. wileyonlinelibrary.com/journal/jtsa DOI: 10.1111/jtsa.12674 Journal of Time Series Analysis published by John Wiley & Sons Ltd. 528 H. DETTE AND P. QUANZ 𝜆′=𝜆⌊n 𝜗n⌋ n: sup 𝜆∈Λn 1 √⌊𝜆⌊n 𝜗n⌋⌋‖‖‖‖‖‖⌊𝜆⌊n 𝜗n⌋⌋ ∑ j=1 𝜂(1) j‖‖‖‖‖‖≤sup 1 n≤𝜆′≤⌊n 𝜗n⌋ n 1 √⌊𝜆′n⌋‖‖‖‖‖‖⌊𝜆′n⌋ ∑ j=1 𝜂(1) j‖‖‖‖‖‖≤sup 1 n≤𝜆≤1 1 √⌊𝜆n⌋‖‖‖‖‖‖⌊𝜆n⌋ ∑ j=1 𝜂(1) j‖‖‖‖‖‖. Regarding the second term in (A11)wehavewith𝜅from (A4’) sup 𝜆∈Λn√𝜆 𝜗n‖‖‖‖‖‖ 1 ⌊𝜆⌊n 𝜗n⌋⌋2⌊𝜆⌊n 𝜗n⌋⌋ ∑ i,j=1 𝜂(1) i⊗𝜂(1) j‖‖‖‖‖‖2 =sup 𝜆∈Λn 𝜆 𝜗n ⌊𝜆⌊n 𝜗n⌋⌋ 1 √𝜆‖‖‖‖‖‖‖‖ 1 √⌊𝜆⌊n 𝜗n⌋⌋⌊𝜆⌊n 𝜗n⌋⌋ ∑ j=1 𝜂(1) j‖‖‖‖‖‖‖‖ 2 ≤OP(n−1∕2)(sup 1 n≤𝜆≤1 1 √⌊𝜆n⌋‖‖‖‖‖‖⌊𝜆n⌋ ∑ j=1 𝜂(1) j‖‖‖‖‖‖)2 =OP(log2∕𝜅(n) √n), by (A1). This implies sup 𝜆∈Λn√𝜆‖‖‖‖ 𝜗nc(1) 𝜆−𝜗0E[𝜂(1) 0⊗𝜂(1) 0]‖‖‖‖2 ≤sup 𝜆∈Λn√𝜆‖‖‖‖‖‖  𝜗n ⌊𝜆⌊n 𝜗n⌋⌋⌊𝜆⌊n 𝜗n⌋⌋ ∑ j=1(𝜂(1) j⊗𝜂(1) j−E[𝜂(1) 0⊗𝜂(1) 0])‖‖‖‖‖‖2 +| 𝜗n−𝜗0|E[𝜂(1) 0⊗𝜂(1) 0]+OP(log2∕𝜅(n) √n) ≤OP(n−1∕2)sup 1 n≤𝜆≤1 1 √⌊𝜆n⌋‖‖‖‖‖‖⌊𝜆n⌋ ∑ j=1(𝜂(1) j⊗𝜂(1) j−E[𝜂(1) 0⊗𝜂(1) 0])‖‖‖‖‖‖2 +OP(log2∕𝜅(n) √n) ≤OP(log2∕𝜅(n) √n), whereweused(3.2) and again a variation of (A1), this time for L2((S×[0,1])2)-valued random variables. It remains to calculate the supremum over the set ΛC n.If⌊𝜆⌊n 𝜗n⌋⌋ ≥⌊n𝜗0⌋is true, then a tedious but straightforward calculation gives c(1) 𝜆=1 ⌊𝜆⌊n 𝜗n⌋⌋⌊n𝜗0⌋ ∑ j=1 𝜂(1) j⊗𝜂(1) j−1 ⌊𝜆⌊n 𝜗n⌋⌋2⌊n𝜗0⌋ ∑ i,j=1 𝜂(1) i⊗𝜂(1) j +1 ⌊𝜆⌊n 𝜗n⌋⌋⌊𝜆⌊n 𝜗n⌋⌋ ∑ j=⌊n𝜗0⌋+1 (𝜂(2) j+𝛿)⊗(𝜂(2) j+𝛿) wileyonlinelibrary.com/journal/jtsa © 2022 The Authors. J. Time Ser. Anal. 44: 505–532 (2023) Journal of Time Series Analysis published by John Wiley & Sons Ltd. DOI: 10.1111/jtsa.12674