Bivariate change point detection: joint detection of changes in expectation and variance

TL;DR

This paper introduces a nonparametric, joint bivariate method for detecting change points in both expectation and variance of univariate data sequences, improving inference over univariate methods.

Contribution

It proposes a novel bivariate moving sum approach for joint detection of changes in mean and variance, supported by a multi-scale analysis and theoretical validation.

Findings

Effective detection of change points in expectation and variance.

Improved inference accuracy compared to univariate methods.

Validated through theoretical analysis and simulations.

Abstract

A method for change point detection is proposed. We consider a univariate sequence of independent random variables with piecewise constant expectation and variance, apart from which the distribution may vary periodically. We aim to detect change points in both expectation and variance. For that, we propose a statistical test for the null hypothesis of no change points and an algorithm for change point detection. Both are based on a bivariate moving sum approach that jointly evaluates the mean and the empirical variance. The joint consideration helps improve inference as compared to separate univariate approaches. We infer on the strength and the type of changes with confidence. Nonparametric methodology supports the analysis of diverse data. Additionally, a multi-scale approach addresses complex patterns in change points and effects. We demonstrate the performance through theoretical results and simulation studies. A companion R-package jcp (available on CRAN) is discussed.

Tables

Table 1. Table 1: Three ways of testing.

test statistic	asymp. level	rejection area $ℛ$ w.r.t. ${\left(J_{h,t}\right)}_t$	assumption
$M_h\left(I\right)$ (9)	$=\alpha$	${𝒞}^c=\left\{{(x,y)}^{𝖳}\in {ℝ}^2\right|d_I\left({(x,y)}^{𝖳}\right)>Q\}$	symmetry
$M_h\left(\hat{\Gamma }\right)$ (11)	$=\alpha$	${\widehat{ℰ}}^c=\hat{\Gamma }\cdot {𝒞}^c$, time dependent	/
$M_h\left(\infty \right)$ (12)		${𝒮}^c=\left\{{(x,y)}^{𝖳}\in {ℝ}^2\right|d_{\infty }\left({(x,y)}^{𝖳}\right)>Q\}$	/

Table 2. Table 2: Performance evaluation. 1000 1000 1000 simulations, T = 1000 𝑇 1000 T=1000 , C = { 200 , 260 , 500 , 720 , 810 } 𝐶 200 260 500 720 810 C=\{200,260,500,720,810\} , α = 5 % 𝛼 percent 5 \alpha=5\% , H = { 50 , 60 , … , 200 } 𝐻 50 60 … 200 H=\{50,60,\ldots,200\} , μ = 11 , 13 , 10 , 8 , 5 , 5 𝜇 11 13 10 8 5 5 \mu=11,13,10,8,5,5 , σ = 1 , 3 , 3 , 3 , 4 , 1.3 𝜎 1 3 3 3 4 1.3 \sigma=1,3,3,3,4,1.3 . N ​ ( μ , σ 2 ) 𝑁 𝜇 superscript 𝜎 2 N(\mu,\sigma^{2}) with ℛ = 𝒞 c ℛ superscript 𝒞 𝑐 \mathscr{R}=\mathscr{C}^{c} , g ​ a ​ m ​ m ​ a ​ ( s , λ ) 𝑔 𝑎 𝑚 𝑚 𝑎 𝑠 𝜆 gamma(s,\lambda) with ℛ = ℰ c ℛ superscript ℰ 𝑐 \mathscr{R}=\mathscr{E}^{c} , u ​ n ​ i ​ f ​ ( a , b ) 𝑢 𝑛 𝑖 𝑓 𝑎 𝑏 unif(a,b) and for p = 2 , 6 , 10 𝑝 2 6 10 p=2,6,10 periodic processes as in Section  3 with ℛ = 𝒞 c ℛ superscript 𝒞 𝑐 \mathscr{R}=\mathscr{C}^{c} . Table: | C ^ | ^ 𝐶 |\hat{C}| and | κ ^ u ( 10 ) | , ( | κ ^ u ( 25 ) | , | κ ^ u ( 5 ) | ) superscript subscript ^ 𝜅 𝑢 10 superscript subscript ^ 𝜅 𝑢 25 superscript subscript ^ 𝜅 𝑢 5 |\hat{\kappa}_{u}^{(10)}|,(|\hat{\kappa}_{u}^{(25)}|,|\hat{\kappa}_{u}^{(5)}|) for u = 1 , … , 5 𝑢 1 … 5 u=1,\ldots,5 .

	$\left|\hat{C}\right|$	$\left|{\hat{\kappa }}_1^{\left(d\right)}\right|$	$\left|{\hat{\kappa }}_2^{\left(d\right)}\right|$	$\left|{\hat{\kappa }}_3^{\left(d\right)}\right|$	$\left|{\hat{\kappa }}_4^{\left(d\right)}\right|$	$\left|{\hat{\kappa }}_5^{\left(d\right)}\right|$
$normal$	4963	957 (996, 853)	845 (919, 705)	698 (883, 531)	854 (974, 687)	943 (994, 841)
$gamma$	4598	970 (999, 854)	853 (917, 710)	662 (879, 513)	811 (935, 687)	476 (584, 381)
$unif$	4908	989 (1000, 916)	844 (906, 706)	707 (893, 536)	869 (953, 713)	987 (1000, 956)
$p=2$	4953	979 (1000, 920)	840 (922, 702)	682 (879, 509)	867 (977, 711)	973 (999, 899)
$p=6$	4958	981 (1000, 899)	863 (931, 731)	701 (884, 551)	856 (975, 691)	973 (998, 890)
$p=10$	4944	988 (1000, 904)	856 (928, 744)	686 (884, 561)	872 (973, 755)	971 (1000, 905)

Equations

$$\mathbb{E}[X_m]=\mathbb{E}[X_1]=:\mu$$

$${\mu }^{⟨k⟩}:=\frac{1}{p}{\sum }_{m=1}^p\mu (m{)}^{⟨k⟩}$$

$${\nu }^2:=\frac{1}{p}{\sum }_{m=1}^p\nu (m{)}^2={\mu }^{\{4\}}-{\sigma }^4,$$

$$\sqrt{n}\left[\left(\begin{array}{c}\hat{\mu }\\ {\hat{\sigma }}^2\end{array}\right)-\left(\begin{array}{c}\mu \\ {\sigma }^2\end{array}\right)\right]\stackrel{d}{⟶}{N}_2(0,\Sigma )\text{with}\Sigma :=\left(\begin{array}{cc}{\sigma }^2& \rho \sigma \nu \\ \rho \sigma \nu & {\nu }^2\end{array}\right),$$

$$X_{1,1},\dots ,X_{1,n{c}_1},X_{2,n{c}_1+1},\dots ,X_{2,n{c}_2},\dots ,X_{|C|+1,n{c}_{|C|}+1},\dots ,X_{|C|+1,nT},$$

$$E_{h,t}^{(n)}:=\frac{{\hat{\mu }}_r-{\hat{\mu }}_{\ell }}{[({\hat{\sigma }}_r^2+{\hat{\sigma }}_{\ell }^2)/(nh){]}^{1/2}}\text{and}{V}_{h,t}^{(n)}:=\frac{{\hat{\sigma }}_r^2-{\hat{\sigma }}_{\ell }^2}{[({\hat{\nu }}_r^2+{\hat{\nu }}_{\ell }^2)/(nh){]}^{1/2}}.$$

$${\hat{\mu }}_j^{⟨k⟩}:=\frac{1}{nh}{\sum }_{i\in I_j}{X}_i^k\text{and}{\hat{\mu }}_j^{\{k\}}:=\frac{1}{nh}{\sum }_{i\in I_j}(X_i-{\hat{\mu }}_j^{⟨1⟩}{)}^k.$$

$$\Gamma :=\left(\begin{array}{cc}1& \rho \\ \rho & 1\end{array}\right)=AD{A}^t:=\left(\begin{array}{cc}1/\sqrt{2}& 1/\sqrt{2}\\ 1/\sqrt{2}& -1/\sqrt{2}\end{array}\right)\cdot \left(\begin{array}{cc}1+\rho & 0\\ 0& 1-\rho \end{array}\right)\cdot \left(\begin{array}{cc}1/\sqrt{2}& 1/\sqrt{2}\\ 1/\sqrt{2}& -1/\sqrt{2}\end{array}\right).$$

$${\mathcal{L}}_{h,t}:={\left(\frac{(W_{t+h}-W_t)-(W_t-W_{t-h})}{\sqrt{2h}},\frac{({\mathcal{W}}_{t+h}-{\mathcal{W}}_t)-({\mathcal{W}}_t-{\mathcal{W}}_{t-h})}{\sqrt{2h}}\right)}^{\mathsf{T}},$$

$$M_h^{(n)}(\Gamma ):=\underset{t\in [h,T-h]}{\max }{d}_{\Gamma }(J_{h,t}^{(n)})\stackrel{d}{⟶}\underset{t\in [h,T-h]}{\max }{d}_I({\mathcal{L}}_{h,t}),$$

$$\hat{\rho }:=\frac{{\hat{\mu }}_r^{\{3\}}+{\hat{\mu }}_{\ell }^{\{3\}}}{({\hat{\sigma }}_r^2+{\hat{\sigma }}_{\ell }^2{)}^{1/2}\cdot ({\hat{\nu }}_r^2+{\hat{\nu }}_{\ell }^2{)}^{1/2}},$$

$$M_h^{(n)}(\hat{\Gamma }):=\underset{t}{\max }{d}_{\hat{\Gamma }}(J_{h,t}^{(n)}).$$

$$M_h^{(n)}(\infty ):=\underset{t}{\max }{d}_{\infty }(J_{h,t}^{(n)}).$$

$$J_{h,c}^{(n)}-{\Delta }_{h,c}^{(n)}\cdot j_{h,c}^{(n)}\stackrel{d}{⟶}{N}_2(0,{\Gamma }_c),$$

$$j_{h,c}^{(n)}:=\left(\begin{array}{c}\frac{{\mu }_2-{\mu }_1}{[({\sigma }_2^2+{\sigma }_1^2)/(nh){]}^{1/2}}\\ \frac{{\sigma }_2^2-{\sigma }_1^2}{[({\nu }_2^2+{\nu }_1^2)/(nh){]}^{1/2}}\end{array}\right)\text{and}{\Delta }_{h,c}^{(n)}:=\left(\begin{array}{cc}{\left(\frac{{\sigma }_2^2+{\sigma }_1^2}{{\hat{\sigma }}_r^2+{\hat{\sigma }}_{\ell }^2}\right)}^{1/2}& 0\\ 0& {\left(\frac{{\nu }_2^2+{\nu }_1^2}{{\hat{\nu }}_r^2+{\hat{\nu }}_{\ell }^2}\right)}^{1/2}\end{array}\right),$$

$${\Gamma }_c:=\left(\begin{array}{cc}1& {\rho }_c\\ {\rho }_c& 1\end{array}\right)\text{with}{\rho }_c:=\frac{{\mu }_2^{\{3\}}+{\mu }_1^{\{3\}}}{({\sigma }_2^2+{\sigma }_1^2{)}^{1/2}\cdot ({\nu }_2^2+{\nu }_1^2{)}^{1/2}}.$$

$${\left[{\tilde{D}}_{h,t}\cdot \left(J_{h,t}^{(n)}-{\Delta }_{h,t}^{(n)}\cdot j_{h,t}^{(n)}\right)\right]}_t\stackrel{d}{⟶}({\mathcal{L}}_{h,t}^{\ast }{)}_t.$$

$${\tilde{\vartheta }}_{\ell }$$

$${\tilde{\mu }}_{\ell }$$

$${\tilde{\sigma }}_{\ell }^2$$

$${\tilde{\nu }}_{\ell }^2$$

$$\displaystyle\qquad\qquad+\frac{t-c}{h}[\mu_{2}^{\{4\}}-4\mu_{2}^{\{3\}}(\tilde{\mu}_{\ell}-\mu_{2})+6\sigma_{2}^{2}(\tilde{\mu}_{\ell}-\mu_{2})^{2}+(\tilde{\mu}_{\ell}-\mu_{2})^{4}]\Big{)}-\tilde{\sigma}_{\ell}^{4},$$

$$e_{h,t}^{(n)}:=\frac{{\tilde{\mu }}_r-{\tilde{\mu }}_{\ell }}{[({\tilde{\sigma }}_r^2+{\tilde{\sigma }}_{\ell }^2)/(nh){]}^{1/2}}\text{and}{v}_{h,t}^{(n)}:=\frac{{\tilde{\sigma }}_r^2-{\tilde{\sigma }}_{\ell }^2}{[({\tilde{\nu }}_r^2+{\tilde{\nu }}_{\ell }^2)/(nh){]}^{1/2}},$$

$${\Delta }_{h,t}^{(n)}:=\left(\begin{array}{cc}{\Delta }_{h,t,1}^{(n)}& 0\\ 0& {\Delta }_{h,t,2}^{(n)}\end{array}\right)\text{and}{\tilde{D}}_{h,t}:=\left(\begin{array}{cc}{\tilde{D}}_{h,t,1}& 0\\ 0& {\tilde{D}}_{h,t,2}\end{array}\right),$$

$$({\Gamma }^{-1/2}\cdot J_{h_1,t}^{(n)},\dots ,{\Gamma }^{-1/2}\cdot J_{h_w,t}^{(n)}{)}_t\stackrel{d}{⟶}({\mathcal{L}}_{h_1,t},\dots ,{\mathcal{L}}_{h_w,t}{)}_t,$$

$$M^{(n)}(\Gamma ):=\underset{h\in H}{\max }{M}_h^{(n)}(\Gamma )=\underset{h\in H}{\max }\underset{t\in [h,T-h]}{\max }{d}_{\Gamma }(J_{h,t}^{(n)})\stackrel{d}{⟶}\underset{h\in H}{\max }\underset{t\in [h,T-h]}{\max }{d}_I({\mathcal{L}}_{h,t}),$$

$${\left(\frac{1}{n}{\sum }_{i=1}^{\lfloor nt\rfloor }{X}_i^k\right)}_t⟶(t{\mu }^{⟨k⟩}{)}_t.$$

$$\displaystyle\frac{1}{n}\sum\nolimits_{i=1}^{n}X_{i}^{k}=\Big{(}\frac{p}{n}\sum\nolimits_{l=0}^{\lfloor n/p\rfloor-1}Z_{l}^{\langle k\rangle}\Big{)}+\Big{(}\frac{1}{n}\sum\nolimits_{i=\lfloor n/p\rfloor p+1}^{n}X_{i}^{k}\Big{)}\longrightarrow\mu^{\langle k\rangle}.$$

$${\mu }^{\{2\}}=\frac{1}{p}{\sum }_{m=1}^p\mu (m{)}^{\{2\}}=\frac{1}{p}{\sum }_{m=1}^p\mu (m{)}^{⟨2⟩}-\frac{1}{p}{\sum }_{m=1}^p\mu (m{)}^2={\mu }^{⟨2⟩}-{\mu }^2,$$

$${\nu }^2=\frac{1}{p}{\sum }_{m=1}^p\nu (m{)}^{\{4\}}=\frac{1}{p}{\sum }_{m=1}^p\mu (m{)}^{\{4\}}-\frac{1}{p}{\sum }_{m=1}^p(\sigma (m{)}^2{)}^2={\mu }^{\{4\}}-{\sigma }^4$$

$$\displaystyle(\hat{\mu}_{\ell}^{\langle k\rangle})_{t}=\Big{(}\frac{1}{nh}\sum\limits_{i=1}^{\lfloor nt\rfloor}X_{i}^{k}\Big{)}_{t}-\Big{(}\frac{1}{nh}\sum\limits_{i=1}^{\lfloor n(t-h)\rfloor}X_{i}^{k}\Big{)}_{t}\longrightarrow\Big{(}\frac{t}{h}\mu^{\langle k\rangle}\Big{)}_{t}-\Big{(}\frac{t-h}{h}\mu^{\langle k\rangle}\Big{)}_{t}=(\mu^{\langle k\rangle})_{t}.\qed$$

$${\left[\frac{1}{\sqrt{n}}{\Gamma }^{-1/2}\cdot \left(\begin{array}{c}{\sum }_{i=1}^{\lfloor nt\rfloor }{Y}_i\\ {\sum }_{i=1}^{\lfloor nt\rfloor }{\mathcal{Y}}_i\end{array}\right)\right]}_t\stackrel{d}{⟶}{\left(\begin{array}{c}{W}_t\\ {\mathcal{W}}_t\end{array}\right)}_t,$$

Full text

Keywords: bivariate, change point detection, jcp, moving sum, multi-scale

1 Introduction

The paper contributes to the field of change point detection which provides methods for the detection of structural breaks – change points – in stochastic sequences. Change point detection finds application in many areas of research, e.g., oceanographic sciences (Killick et al., , 2010), neuroimaging (Aston and Kirch, , 2012), telecommunication (Zhang et al., , 2009), DNA sequencing (Braun et al., , 2000), econometrics (Zeileis et al., , 2010), to name a few. Change point problems are studied extensively, thus, many aspects of statistical methodology are addressed, e.g., hypothesis testing, change point estimation, model complexity, computational feasibility, data structures, practical performance, etc., see for example the textbooks Basseville and Nikiforov, (1993); Csörgő and Horváth, (1997); Brodsky and Darkhovsky, (1993); Chen and Gupta, (2000); Brodsky, (2017) or the review articles Aue and Horváth, (2013); Jandhyala et al., (2013).

We study independent univariate random variables (RVs) that have piecewise constant expectation ($\mu$) and variance ($\sigma^{2}$). The goal is to detect changes in both $\mu$ and $\sigma^{2}$. Apart from constant first and second moments within each section, the distribution is allowed to vary periodically, allowing for variability in higher moments. This periodicity can improve the mimicking of real data as compared to i.i.d. RVs, the latter of which are contained as a special case. We provide a nonparametric method for the detection of change points which may occur on multiple time scales. There is vast literature on all aspects mentioned. Regarding independence, we mention Csörgő and Horváth, (1988); Gombay and Horváth, (2002); Horváth and Hušková, (2005) ($U$-statistics), Horváth and Shao, (2007); Holmes et al., (2013) (empirical processes), and Siegmund, (1988); Frick et al., (2014); Fang et al., (2020) (likelihood ratios), the last two tackling multiple change points. For the detection of changes in $\mu$ we mention the nonparametric methods of Wolfe and Schechtman, (1984); Horváth et al., (2008); Jarušková, (2010); or Dehling et al., (2013), and in particular regarding change points on multiple time scales we refer to Spokoiny, (2009); Fryzlewicz, (2014); Matteson and James, (2014). For the detection of changes in $\sigma^{2}$ we mention the articles of Hsu, (1977); Inclán and Tiao, (1994); Chen and Gupta, (1997); Whitcher et al., (2000); Killick et al., (2013); Korkas and Fryzlewicz, (2017) and a paper for changes in scale by Gerstenberger et al., (2020). Combining both, changes in $\mu$ and $\sigma^{2}$, we refer to Pein et al., (2017) who aim at detecting changes in $\mu$ by allowing $\sigma^{2}$ to change simultaneously. Vice versa, Gao et al., (2019) or Dette et al., (2015) study changes in $\sigma^{2}$ allowing $\mu$ to vary smoothly. Further, we note that Górecki et al., (2018) and Messer et al., (2014) aim at detecting changes in $\mu$ (first moment) while a certain degree of heteroscedasticity, i.e., variability in the second moment, is allowed. In this work, the latter is extended in the sense that we detect changes in both $\mu$ and $\sigma^{2}$, while we also allow for variability in higher moments.

We propose a bivariate method to jointly quantify change in $\mu$ and $\sigma^{2}$. For that we consider two moving sum processes (MOSUM): data are pointwise restricted to adjacent windows from which first the means and second the empirical variances are compared. The first statistic is sensitive to changes in $\mu$ and the second to changes in $\sigma^{2}$. See Figure 1 for an example of the two univariate processes (B,C) and their joint consideration (D), details are explained later. Regarding univariate MOSUM we mention Steinebach and Eastwood, (1995); Antoch and Hušková, (1999); Hušková and Slabý, (2001) and Eichinger and Kirch, (2018). Importantly, the bivariate approach presented here helps overcome flawed inference as compared to separate univariate approaches. Also, it enables a straightforward interpretation of the types of changes ($\mu$, $\sigma^{2}$ or both), as well as their strengths (effect sizes), while also controlling the error (confidence) of such statements. Methodologically, we construct first a test for the null hypothesis of an absence of change points. Second, we propose an algorithm for change point detection which is run if the null hypothesis is rejected.

We mention three benefits: first, the proposed theory is highly nonparametric which enables wide-ranging applicability. Strong performance is shown under different distributional assumptions, including normal-, exponential- or gamma-distributed data. Second, multi-scale aspects are captured, meaning that the occurrence of fast as well as slow change points with different strengths of effects are tackled. Third, the method is ready to use. Theory is proven with all unknown quantities replaced by appropriate estimators. A companion R-package jcp (joint change point detection, Messer, (2020)) is provided whose graphical output facilitates interpretation.

The paper is organized as follows: in Section 2 we present basic concepts. In Section 3 we introduce the model, and in Section 4 we define the MOSUM statistics. We construct the test in Section 5 and discuss change point detection in Section 6. In Section 7 we give additional theory. In Section 8 practical aspects are considered: methodology is extended, the R-package jcp is discussed, simulation studies are performed, and a real data example is shown. Proofs and auxiliary results are given in the Appendix.

2 The idea of testing and change point detection

We show how MOSUM processes are used for testing and change point detection. Particularly, we motivate the bivariate aspect. For that, consider the three processes in Figure 1, differentiated in panels I, II, and III.

We consider $T=1000$ independent RVs that are piecewise $N(\mu,\sigma^{2})$ distributed, see red points in segments A and top D (coinciding). In panel I, there are no change points, $C=\emptyset$. In panels II and III there are three change points $C=\{250,500,750\}$. In both II and III the types of changes coincide: a change in $\mu$ at $250$ (purple cross), an increase in $\sigma^{2}$ at $500$ (blue cross), and a decrease in both $\mu$ and $\sigma^{2}$ at $750$ (purple and blue cross). II and III only differ at $250$: there is a small increase in $\mu$ in II and a prominent decrease in $\mu$ in III.

MOSUM processes are shown in segments B, C, and D. Two adjacent windows of size $h=70$ (see panel IA) are shifted through time and statistics are evaluated locally from the RVs in the windows. The first statistics $E_{h,t}$ (B, magenta) is Student’s unpooled $t$-statistic, i.e., Welch’s statistic, which compares the empirical means, see (5). It is sensitive to changes in $\mu$. The second statistic $V_{h,t}$ (C, blue) similarly compares the emprical variances and is thus sensitive to changes in $\sigma^{2}$. In segment D we see a bivariate process (black) given by the joint statistic $J_{h,t}:=(E_{h,t},V_{h,t})^{\mathsf{T}}$, where $\mathsf{T}$ denotes transposing. If the windows do not overlap a change point, then the two estimates from the left and right window typically resemble each other, resulting in a statistic close to zero. A strong deviation from zero indicates a change. We see brown rejection boundaries (B, C horizontal lines at $\approx\pm 3.59$, D circle with radius $\approx 4$. The boundaries coincide between the panels but appear different due to the scaling of axes). The idea is that the processes entirely lie within the boundaries with a predefined probability of $1-\alpha$, here $95\%$, in case null hypothesis $\mathcal{H}_{0}:C=\emptyset$ holds true, see panel I and Section 5 for details. If the boundaries are crossed at some point (II and III in B,C and D), the null hypothesis is rejected. After rejection, the set $C$ is estimated (green bars) via successive $\arg\!\max$ estimation: find the largest deviation form zero, take the argument as a change point estimate, delete the process in the $h$-neighborhood of the estimate and repeat until the remaining process lies within the boundaries (also see Figure 4). In Figure 1 panel II inference is reasonable as in B the two changes in $\mu$, in segment C the two changes in $\sigma^{2}$, and in D all three change points were detected. In III expectation change detection succeeds in B, but unfortunately variance inference in segment C fails as two changes in $\sigma^{2}$ are falsely estimated. This is caused by the prominent change in $\mu$. For an intuition, recall that the empirical variance is the mean squared deviation from the mean. Thus, when $\mu$ changes, this impacts the mean and thus the variance. This problem also appears in panel II, but it is practically negligible as the change in $\mu$ is small. Vice versa, changes in $\sigma^{2}$ practically do not impact expectation change detection, intuitively, because the mean is the first moment and not affected if only the second moment changes. Statements about robustness are subject of Section 7. One way to tackle false $\sigma^{2}$-inference is to incorporate information about $\mu$: if the empirical variances are centered correctly then $\sigma^{2}$-estimation will not systematically react falsely to changes in $\mu$. In the context of stochastic point processes this was studied in Albert et al., (2017). A second way of treating this problem is presented in this paper: the joint observation of both processes. In panel IIID the bivariate MOSUM overcomes the problem and successfully detected three change points.

The golden dartboards in segment D show the asymptotic distribution of $J_{h,c}$ at a true change point $c$. This distribution is bivariate normal, the golden cross marks the expectation and the circles the $66\%$- and $95\%$-contour lines, see Proposition 6.2. When the sliding windows run into $c$, then $(J_{h,t})_{t}$ starts an excursion approaching the center of the dartboard. The red crosses ($+$) mark $J_{h,c}$ i.e., they can be considered a realization from a dartboard. When the sliding windows trespass $c$, the process returns to zero fluctuation. Glimpse at Figure 5 for the systematics of the excursions. The type and the strength of the change affect the systematic aspect of the excursion: if there is a change only in $\mu$, the process tends to leave the circle in horizontal direction, to the right if there is an increase in $\mu$ and to the left in case of a decrease. Indeed, at time $250$ the dartboard is shifted only along the abscissa. In panel II it is found at the right side as there is an increase in $\mu$, and in III it is shifted left due to the decrease in $\mu$. Further, in panel III it lies far out as the change in $\mu$ is strong. If there is a change only in $\sigma^{2}$, the process systematically moves along the ordinate, upwards in case of an increase and downwards in case of decrease. At $500$ the dartboard is shifted only along the ordinate. In case of a change in both $\mu$ and $\sigma^{2}$, the process leaves the circle in both horizontal and vertical direction. Indeed, at $750$ the dartboard lies south-west due to the decline in both $\mu$ and $\sigma^{2}$. Vice versa, the location of a dartboard facilitates change point interpretation, see Section 6.

3 The Model

Auxiliary processes

Consider a single probability space $(\Omega,\mathcal{A},\mathbb{P})$ throughout.

Definition 3.1**.**

A sequence $\mathbf{X}=(X_{i})_{i=1,2,\ldots}$ of independent RVs of $\mathscr{L}^{4}(\Omega,\mathcal{A},\mathbb{P})$ is called an auxiliary process, if there exists a $p\in\{1,2,\ldots\}$ such that for all $m\in\{1,2,\ldots,p\}$ it holds that $X_{ip+m}\stackrel{{\scriptstyle d}}{{=}}X_{m}$ for all $i=1,2\ldots$, and also that

[TABLE]

The distribution varies with period $p$ while the first two moments are constant. This includes i.i.d. sequences. Higher moments may vary periodically. Thus, for $k=1,\ldots,4$ we abbreviate the $k$-th moment and centered moment via $\mu(m)^{\langle k\rangle}:=\mathbb{E}[X_{m}^{k}]$ and $\mu(m)^{\{k\}}:=\mathbb{E}[(X_{m}-\mu)^{k}]$. We set averaged moments within a period

[TABLE]

For $k=1,2$, (1) implies constants $\mu(m)^{\langle 1\rangle}=\mu$, $\mu(m)^{\langle 2\rangle}=\mu^{\langle 2\rangle}$, $\mu(m)^{\{1\}}=0$ and $\mu(m)^{\{2\}}=\sigma^{2}$. Further set $\nu(m)^{2}:=\mathbb{V}\!ar((X_{m}-\mu)^{2})=\mu(m)^{\{4\}}-\sigma^{4}$ and the average

[TABLE]

Note that $\sigma^{2}>0$ implies $\nu(m)^{2}>0$ and thus $\nu^{2}>0$. Further note $\mu(m)^{\{3\}}=\mathbb{C}ov(X_{m},(X_{m}-\mu)^{2})$. It is $|\rho|<1$ as $(X_{m}-\mu)^{2}$ is not linear in $X_{m}$. Call the averages $\mu$, $\sigma^{2}$, $\nu^{2}$, $\rho$, $\mu^{\{k\}}$ and $\mu^{\langle k\rangle}$ in (1), (2) and (3) the population parameters. Condition (1) is sufficient for moment estimators to appropriately estimate all population parameters – the periodicity is averaged out over large samples.

For $d=1,2,\dots$ let $N_{d}(w,\Sigma)$ denote the $d$-variate normal distribution with expectation $w\in\mathbb{R}^{d}$ and $d\times d$ covariance matrix $\Sigma$. For the mean $\hat{\mu}:=(1/n)\sum_{i=1}^{n}X_{i}$ and the empirical variance $\hat{\sigma}^{2}:=(1/n)\sum_{i=1}^{n}(X_{i}-\hat{\mu})^{2}$ we jointly obtain as $n\to\infty$

[TABLE]

while $\stackrel{{\scriptstyle d}}{{\longrightarrow}}$ denotes convergence in distribution, see Proposition 6.2. $\Sigma$ is regular as $det(\Sigma)=(1-\rho)\sigma^{2}\nu^{2}>0$.

We give an example: let $Y$ be a mix of uniform RVs with density $f_{Y}(y)=(3/(4d))\mathbbm{1}_{[\mu-d,\mu]}(y)+(1/(4d))\mathbbm{1}_{[\mu+d,\mu+2d]}(y)$ for $\mu\in\mathbb{R}$ and $d>0$. Then $\mathbb{E}[Y]=\mu$ and $\mathbb{V}\!ar(Y)=5d^{2}/6=:\sigma^{2}$, and $f_{Y}$ has positive skewness as $3>1$. Further, let $Z$ have density $f_{Z}(z)=f_{Y}(2\mu-z)$, i.e., $f_{Z}$ results from $f_{Y}$ by mirroring at $\mu$. Thus, also $\mathbb{E}[Z]=\mu$ and $\mathbb{V}\!ar(Z)=\sigma^{2}$, but skewness switches the sign. Consider independent copies of $Y$ and $Z$. Set $p$ even. Repeatedly let $p/2$ copies of $Z$ follow $p/2$ of $Y$. This is an auxiliary process with period $p$, see Figure 2 $\mathbf{X}_{1},\ldots,\mathbf{X}_{|C|+1}$.

The Model $\mathcal{M}$

A piecewise combination of auxiliary processes constitutes a process with change points: we fix $T\in\mathbb{N}\backslash\{0,1\}$ and consider a subset $C\subset\{1,\ldots,T-1\}$ of cardinality $|C|$ with ordered elements $c_{1}<c_{2}<\ldots<c_{|C|}$. Call $c_{u}$ the $u$-th change point and $C$ the set of change points. Given $C$, we consider $|C|+1$ independent auxiliary processes $\mathbf{X}_{1},\ldots,\mathbf{X}_{|C|+1}$, with $\mathbf{X}_{u}=(X_{u,i})_{i=1,2,\ldots}$ and population parameters indexed by $u$, e.g., $\mu_{u}:=\mathbb{E}[X_{u,1}]$ and $\sigma_{u}^{2}:=\mathbb{V}\!ar(X_{u,1})$, for $u=1,\ldots,|C|+1$. We assume $(\mu_{u},\sigma_{u}^{2})\not=(\mu_{u+1},\sigma_{u+1}^{2})$, for $u=1,\ldots,|C|$, i.e., at any $c_{u}$ either the expectation or the variance (or both) change, while other parameters may change too. For asymptotics we let $T$ and $c_{1},\ldots,c_{|C|}$ depend on a factor $n$, i.e., we switch to $nT$ and $nc_{1},\ldots,nc_{|C|}$. For $n=1,2,\ldots$ we define a compound process

[TABLE]

i.e., after $nc_{u}$ we enter $\mathbf{X}_{u+1}$ with $(\mu_{u+1},\sigma_{u+1}^{2})$. Given $T$, the set $\mathcal{M}$ of such processes constitutes the model. For asymptotics we let $n\to\infty$ and treat $C$, as well as all population parameters including the periods, as fixed, i.e., not depending on $n$, but unknown. By increasing $n$, the time we remain in each $\textbf{X}_{u}$ increases linearly, constituting a triangular setup. We call $n=1$ the real time scenario, see Figure 2.

For $\mathbf{X}\in\mathcal{M}$ we want to test $\mathcal{H}_{0}:C=\emptyset$. In case of rejection we aim at estimating $C$. If $C=\emptyset$ we omit indices $\mu:=\mu_{1}$, $\sigma^{2}:=\sigma_{1}^{2}$, $\nu^{2}:=\nu_{1}^{2}$ etc.

4 The moving sum processes

For $\textbf{X}\in\mathcal{M}$ we study moving sum processes in which we locally evaluate the RVs restricted to a left window $(t-h,t]$ and an adjacent right window $(t,t+h]$. A comparison of means is sensitive to changes in the expectation and a comparison of empirical variances is sensitive to changes in the variance.

Definition of $(E_{h,t}^{(n)})_{t}$ and $(V_{h,t}^{(n)})_{t}$ via local parameter estimators

Let $h\in\{2,3,\ldots,\lfloor T/2\rfloor\}$ be a window size independent of $n$, while $\lfloor\cdot\rfloor$ denotes the floor function, and let $t\in[h,T-h]$. For $n=1,2,\ldots$ set

[TABLE]

We define the estimators in (5). The subscripts $\ell$ and $r$ indicate local evaluation, while $I_{\ell}:=\{\lfloor n(t-h)\rfloor+1,\ldots,\lfloor nt\rfloor\}$ and $I_{r}:=\{\lfloor nt\rfloor+1,\ldots,\lfloor n(t+h)\rfloor\}$ are the indices associated with the windows. We strengthen the dependence on $n$, $t$ and $h$, which is inherited to the estimators below, but omitted for simplicity. For $j\in\{\ell,r\}$ and $k\in\{1,\ldots,4\}$ we define estimators for the moments in (2) and (3) via

[TABLE]

and further $\hat{\mu}_{j}:=\hat{\mu}_{j}^{\langle 1\rangle}$, $\hat{\sigma}_{j}^{2}:=\hat{\mu}_{j}^{\{2\}}$ and $\hat{\nu}_{j}^{2}:=\hat{\mu}_{j}^{\{4\}}-\hat{\sigma}_{j}^{4}$.

We study statistics in function space. For an interval $[a,b]\subset\mathbb{R}_{0}$ let $(\mathcal{D}_{\mathbb{R}^{d}}[a,b],d_{SK})$ and $(\mathcal{D}_{\mathbb{R}^{d}}[a,b],\|\cdot\|)$ denote the spaces of $\mathbb{R}^{d}$-valued càdlàg-functions on $[a,b]$ equipped with Skorokhod topology $d_{SK}$ or the supremum norm $\|\cdot\|$. Convergence w.r.t. $\|\cdot\|$ implies convergence regarding $d_{SK}$. Analogously, for $[0,\infty)$ consider $(\mathcal{D}_{\mathbb{R}^{d}}[0,\infty),d_{SK})$. For $n=1,2,\ldots$ the estimators in (6) and thus $(E_{h,t}^{(n)})_{t}$ and $(V_{h,t}^{(n)})_{t}$ constitute processes in $(\mathcal{D}_{\mathbb{R}}[h,T-h],d_{SK})$. Omitting the superscript $n$ abbreviates $n=1$, e.g., $E_{h,t}:=E_{h,t}^{(1)}$. See Figure 1 for $(E_{h,t})_{t}$ (B, magenta) and $(V_{h,t})_{t}$ (C, blue) evaluated from $\mathbb{X}$ (A, red). If $C=\emptyset$ then the estimators are functionally strongly consistent for their population parameters:

Lemma 4.1**.**

Let $\mathbf{X}\in\mathcal{M}$ with $C=\emptyset$. For $j\in\{\ell,r\}$ it holds in $(\mathcal{D}_{\mathbb{R}}[h,T-h],d_{\|\cdot\|})$ as $n\to\infty$ almost surely that $(\hat{\mu}_{j}^{\langle k\rangle})_{t}\to(\mu^{\langle k\rangle})_{t}$ and $(\hat{\mu}_{j}^{\{k\}})_{t}\to(\mu^{\{k\}})_{t}$, for $k\in\{1,\ldots,4\}$.

Consequently, $(\hat{\sigma}_{j}^{2})_{t}\to(\sigma^{2})_{t}$ and $(\hat{\nu}_{j}^{2})_{t}\to(\nu^{2})_{t}$ as $n\to\infty$ almost surely (a.s.).

The joint process $(J_{h,t}^{(n)})_{t}$

We consider the joint process $(J_{h,t}^{(n)})_{t}$ via $J_{h,t}^{(n)}:=(E_{h,t}^{(n)},V_{h,t}^{(n)})^{\mathsf{T}}$ in $(\mathcal{D}_{\mathbb{R}^{2}}[h,T-h],d_{SK})$, see Figure 1D. Weak convergence to a bivariate Gaussian process is shown under $C=\emptyset$, see Proposition 4.2 below. From this the test is constructed in Section 5. Further, convergences are extended to $C\not=\emptyset$ in Propositions 6.2 and 7.1, supporting change point detection in Section 6.

In the remainder of this section let $C=\emptyset$. In general, the components of $(J_{h,t}^{(n)})_{t}$ are correlated as both rely on X: in Figure 3A $(J_{h,t})_{t}$ primarily varies along the diagonal $x=y$. Symmetry of the RVs $\rho=\mu^{\{3\}}/(\sigma\nu)=0$ (see (3)) is necessary and sufficient for the components to be asymptotically independent. Vice versa, skewness $\rho\not=0$ results in correlated components. We abusively speak about symmetry and skewness, meaning the average $\mu^{\{3\}}=(1/p)\sum_{m=1}^{p}\mu(m)^{\{3\}}$ to be $=0$ or $\not=0$.

Skewness is captured in the correlation matrix

[TABLE]

$\Gamma$ is regular as $|\rho|<1$. The eigenvalue decomposition (7) yields $\Gamma^{1/2}$ and $\Gamma^{-1/2}$. Unit diagonals of $\Gamma$ imply that the columns of $A$ span the diagonals $x=y$ and $x=-y$ in $\mathbb{R}^{2}$. Symmetry $\rho=0$ means $\Gamma=I=\Gamma^{1/2}=\Gamma^{-1/2}$, with identity $I$. For an example of a skewed distribution consider $gamma(s,\lambda)$ with shape $s>0$ and rate $\lambda>0$, see also Figure 3A. We obtain $\sigma^{2}=s/\lambda^{2}$, $\nu^{2}=2[s^{2}+3s]/\lambda^{4}$ and $\mu^{\{3\}}=2s/\lambda^{3}$, and $\rho=(2/[s+3])^{1/2}>0$ depending only on $s$.

Proposition 4.2**.**

Let $\mathbf{X}\in\mathcal{M}$ with $C=\emptyset$. In $(\mathcal{D}_{\mathbb{R}^{2}}[h,T-h],d_{SK})$ it holds as $n\to\infty$ that $(\Gamma^{-1/2}\cdot J_{h,t}^{(n)})_{t}\stackrel{{\scriptstyle d}}{{\longrightarrow}}(\mathcal{L}_{h,t})_{t}$.

The bivariate limit process $(\mathcal{L}_{h,t})_{t}$ is given via

[TABLE]

while $(W,\mathcal{W}):=(W_{t},\mathcal{W}_{t})_{t\geq 0}$ denotes a planar Brownian motion. The double windows are preserved. $(\mathcal{L}_{h,t})_{t}$ is a continuous $2h$-dependent bivariate Gaussian process that is isotropic. It is $\mathcal{L}_{h,t}\sim N_{2}(0,I)$ for all $t$.

5 The statistical test

We test $\mathcal{H}_{0}:C=\emptyset$. A large deviation of $(J_{h,t}^{(n)})_{t}$ from zero speaks against $C=\emptyset$. We apply Proposition 4.2 to derive a rejection boundary from $(\mathcal{L}_{h,t})_{t}$ in simulations. Symmetry implies asymptotic isotropy of $(J_{h,t}^{(n)})_{t}$ and we set the boundary as a circle, see Figure 1. Generally, for possibly skewed RVs, $(J_{h,t}^{(n)})_{t}$ primarily varies along $x=y$ or $x=-y$, resulting in a rejection ellipse or square, see Figure 3.

Derivation of the test

Let $C=\emptyset$. For $(x,y)^{\mathsf{T}}\in\mathbb{R}^{2}$ let $d_{I}((x,y)^{\mathsf{T}})=(x^{2}+y^{2})^{1/2}$ denote the Euclidean and $d_{\Gamma}((x,y)^{\mathsf{T}})=[(x,y)\cdot\Gamma^{-1}\cdot(x,y)^{\mathsf{T}}]^{1/2}$ the Mahalanobis distance. Noting $d_{I}(\Gamma^{-1/2}\cdot J_{h,t}^{(n)})=d_{\Gamma}(J_{h,t}^{(n)})$, Proposition 4.2 and continuous mapping imply convergence of the maximum, as $n\to\infty$

[TABLE]

We use $M_{h}(\Gamma):=M_{h}^{(1)}(\Gamma)$ as a test statistic. We reject iff $M_{h}(\Gamma)$ exceeds the $(1-\alpha)$-quantile $Q$ of the limit distribution in (9), given $\alpha\in(0,1)$. To the best of our knowledge there is no closed formula available. We choose $Q$ as a quantile of the approximated distribution derived in Monte Carlo simulations.

In a similar setup Jarušková and Piterbarg, (2011) derived tail approximations for functionals of Brownian bridges to adjust $Q$. Here, $[d_{I}^{2}(\mathcal{L}_{h,t})]_{t}$ constitutes a $\chi^{2}$-process, which should allow to derive tail bounds as well, see Albin, (1990); Lindgren, (1980); Adler, (1990); Talagrand, (2014). We mention high accuracy of simulations.

Equivalently, instead of comparing the maximum $M_{h}(\Gamma)$ to $Q\in\mathbb{R}$, we can judge the entire process $(J_{h,t})_{t}$ w.r.t. a rejection area $\mathscr{R}\subset\mathbb{R}^{2}$. For that, define an ellipse $\mathscr{E}:=\{(x,y)^{\mathsf{T}}\in\mathbb{R}^{2}\;|\;d_{\Gamma}((x,y)^{\mathsf{T}})\leq Q\}$, see Figure 3A. $\mathcal{H}_{0}$ is rejected iff $(J_{h,t})_{t}$ enters $\mathscr{R}=\mathscr{E}^{c}$ at any time, while the superscript $c$ indicates the complement in $\mathbb{R}^{2}$. In case of symmetry $\Gamma=I$, $\mathscr{E}$ equals a circle $\mathscr{C}$ with radius $Q$, i.e., $\mathscr{C}:=\{(x,y)^{\mathsf{T}}\in\mathbb{R}^{2}\;|\;d_{I}((x,y)^{\mathsf{T}})\leq Q\},$ thus $\mathscr{R}=\mathscr{C}^{c}$, see Figure 1. For the skewed case, $\mathscr{E}$ derives from $\mathscr{C}$ by squeezing, $\mathscr{E}=\Gamma^{1/2}\cdot\mathscr{C}=\{\Gamma^{1/2}\cdot(x,y)^{\mathsf{T}}\;|\;(x,y)^{\mathsf{T}}\in\mathscr{C}\}$.

Treatment of the unknown correlation $\rho$

$M_{h}^{(n)}(\Gamma)$ depends on unknown $\rho=\mu^{\{3\}}/(\sigma\nu)$. We propose three ways to cope with $\rho$ in practice. First, the symmetry assumption $\rho=0$ yields $M_{h}^{(n)}(I)$. Second, we consistently estimate $\rho$ locally by

[TABLE]

recall (6). We define $\hat{\Gamma}$ by replacing $\rho$ with $\hat{\rho}$ in (7), and consider the statistic

[TABLE]

Convergence of $M_{h}^{(n)}(\hat{\Gamma})$ as in (9) holds true because in $(D_{\mathbb{R}}[h,T-h],d_{SK})$ it is $(\hat{\rho})_{t}\to(\rho)_{t}$ a.s., see Lemma (4.1). In terms of $(J_{h,t}^{(n)})_{t}$, the rejection ellipse $\hat{\mathscr{E}}:=\hat{\Gamma}^{1/2}\cdot\mathscr{C}$ is time dependent. Under $\mathcal{H}_{0}$ estimation is consistent and the asymptotic $\alpha$-level is kept. Third, we propose a conservative approach to avoid estimation. Define a square $\mathscr{S}:=\{(x,y)^{\mathsf{T}}\in\mathbb{R}^{2}\;|\;d_{\infty}((x,y)^{\mathsf{T}})\leq Q\}$, while $d_{\infty}((x,y)^{\mathsf{T}}):=\max(|x|,|y|)$ denotes the maximum norm. $\mathscr{S}$ has center zero and edge length equal to the diameter of $\mathscr{C}$. Because $\Gamma$ has unit diagonals, any ${\mathscr{E}}$ is trapped in $\mathscr{S}$, see Figure 3A. Thus, with $\mathscr{R}:=\mathscr{S}^{c}$ the asymptotic probability to falsely reject is less than $\alpha$. The test statistic is

[TABLE]

This approach is optimal in the sense that $\mathscr{E}$ touches each edge of $\mathscr{S}$ exactly once, thus $\mathscr{S}$ can not be shrunk. The tests are summarized as follows, see also Table 1.

Theorem 5.1**.**

Let $\mathbf{X}\in\mathcal{M}$, and for $\alpha\in(0,1)$ let $Q$ be the $(1-\alpha)$-quantile of the limit in (9). Assume $\mathcal{H}_{0}:C=\emptyset$. Then it holds $\lim_{n\to\infty}\mathbb{P}(M_{h}^{(n)}(\Gamma)>Q)=\alpha$. Further, it holds first $\lim_{n\to\infty}\mathbb{P}(M_{h}^{(n)}(I)>Q)=\alpha$ given $\rho=0$, second $\lim_{n\to\infty}\mathbb{P}(M_{h}^{(n)}(\hat{\Gamma})>Q)=\alpha$, and third $\lim_{n\to\infty}\mathbb{P}(M_{h}^{(n)}(\infty)>Q)\leq\alpha$.

Note that Proposition 6.2 below yields asymptotic power one. Further, the bivariate approach avoids double testing unlike the two univariate approaches from Section 2: Set $\alpha=5\%$. With the univariate approaches we obtain the boundaries $\approx\pm 3.59$, by simulation of the quantile of the temporal maximum of a single component of $(\mathcal{L}_{h,t})_{t}$, see Figure 1B and C. With the bivariate procedure we get $Q\approx 4.00$ (radius in D). This slight increase is the price for avoiding $\alpha$-error accumulation.

6 Change point detection

After rejection of $\mathcal{H}_{0}$ we detect change points via successive $\arg\max$ estimation, see Figure 4. Proposition 6.2 below states asymptotic normality of $J_{h,c}^{(n)}$ at $c\in C$, see dartboards in Figures 1 and 3. This facilitates inference about the effects, i.e., the strength and the type of the change, and enables interpretation in practice.

Algorithm 6.1**.**

Set $\mathscr{R}$ as either $\mathscr{C}^{c}$, $\hat{\mathscr{E}}^{c}$ or $\mathscr{S}^{c}$ as used in the test. Set $\hat{C}_{h}:=\emptyset$ and $\tau_{h}:=[h,T-h]\cap\mathbb{N}$. While $J_{h,t}\in\mathscr{R}$ for any $t\in\tau_{h}$, update $\hat{C}_{h}$ and $\tau_{h}$ as follows: among all $t\in\tau_{h}$ for which $J_{h,t}\in\mathscr{R}$, choose the candidate $\hat{c}$ that maximizes the Euclidean distance of $J_{h,t}$, add $\hat{c}$ to $\hat{C}_{h}$ and delete its $h$-neighborhood $\{\hat{c}-h+1,\ldots,\hat{c}+h\}$ from $\tau_{h}$.

Note that for any choice of $\mathscr{R}$, either $\mathscr{C}^{c}$, $\mathscr{E}^{c}$ or $\mathscr{S}^{c}$, it is always the Euclidean distance that inference is based on, see interpretation below. For single $E_{h,t}^{(n)}$ see similar procedures in Antoch and Hušková, (1999) and Eichinger and Kirch, (2018).

Asymptotic normality of $J_{h,c}^{(n)}$

Proposition 6.2**.**

Let $\mathbf{X}\in\mathcal{M}$, $h$ be a window size, and $c\in[h,T-h]$ a change point such that $(c-h,c+h]\backslash\{c\}$ is free from change points. Then it holds as $n\to\infty$

[TABLE]

with $\Delta_{h,c}^{(n)}$ and $j_{h,c}^{(n)}$ in (13) and $\Gamma_{c}$ in (14).

Near $c:=c_{u}$ the process $\mathbf{X}$ derives from $\mathbf{X}_{u}$ and $\mathbf{X}_{u+1}$. W.l.o.g. let $u=1$. Set

[TABLE]

and call $j_{h,c}^{(n)}$ the asymptotic expectation of $J_{h,c}^{(n)}$. Further, $\Delta_{h,c}^{(n)}$ replaces the true scaling of $J_{h,c}^{(n)}$ with the local estimators. Lemma 4.1 implies $\Delta_{h,c}^{(n)}\to I$ a.s. componentwise as $n\to\infty$. The asymptotic correlation matrix is

[TABLE]

Note the analogy of $\rho_{c}$ and its estimator $\hat{\rho}$ in (10). In Figures 1 and 3 the distribution of $J_{h,c}$ is approximated by $N_{2}(j_{h,c}^{(1)},\Gamma_{c})$ for the three change points. The golden dartboards describe the $66\%$- and $95\%$-contour lines with center $j_{h,c}^{(1)}$. A stronger $|\rho_{c}|$ increases the squeezing of the dartboards in Figure 3B. Symmetry $\rho_{c}=0$ implies circles in Figure 1. The size of a dartboard is unique up to squeezing as $\Gamma_{c}$ has unit diagonals: $N_{2}(j_{h,c}^{(1)},\Gamma_{c})$ implies that for $\alpha\in(0,1)$ any $\alpha$-contour ellipse is trapped optimally in a square with edge lengths $2\sqrt{q_{\alpha}}$, while $q_{\alpha}$ denotes the $\alpha$-quantile of the $\chi^{2}$-distribution with two degrees of freedom.

Change point interpretation

To interpret effects we decompose $j_{h,c}^{(n)}=\sqrt{nh}\cdot d_{I}(j_{c})\cdot(\cos(\omega),\sin(\omega))^{\mathsf{T}}$, with $j_{c}:=j_{h,c}^{(n)}/\sqrt{nh}$, and $\omega\in[0,2\pi)$ the angle between $j_{c}$ and the abscissa. Then $j_{c}$ is a signal to noise ratio that captures the effects: the Euclidean distance $d_{I}(j_{c})$ measures the strength of the change, and the angle $\omega$ the type of the change, with $\omega\in\{0,\pi\}\Leftrightarrow$ change only in expectation $(0=\textrm{increase},\pi=\textrm{decrease})$, and $\omega\in\{\pi/2,(3/2)\pi\}\Leftrightarrow$ change only in variance $(\pi/2=\textrm{increase},(3/2)\pi=\textrm{decrease})$. Generally, a dartboard is shifted towards $(cos(\omega),sin(\omega))^{\mathsf{T}}$, see Figure 1 or 3. A dartboard lies far off of zero if either $\sqrt{nh}$ or $d_{I}(j_{c})$ is large. This implies asymptotic power one, also note extensions to local changes where the parameters imply $d_{I}(j_{c})=o(1/\sqrt{nh})$, see e.g., Antoch and Hušková, (1999). Further note that the location does not depend on the symmetry of $\mathbf{X}$ which is captured in $\Gamma_{c}$. Thus, in Algorithm 6.1 it is always the Euclidean distance that $J_{h,t}^{(n)}$ is judged by, regardless of $\mathscr{R}$.

In practice, for an estimate $\hat{c}$ interpretation is based on $J_{h,\hat{c}}$, see green arrows in Figure 4. We assign the $66\%$- and $95\%$-contour lines of $N_{2}(J_{h,\hat{c}},\hat{\Gamma}_{\hat{c}})$ with $\hat{\Gamma}_{\hat{c}}$ as in (10), see green dartboards in Figure 4D. This strengthens confidence in interpretation: it is plausible that the smallest estimate refers to a change only in $\mu$ (right dartboard close to abscissa) and that the middle one indicates a change only in $\sigma^{2}$ (upper dartboard close to ordinate). A single parameter change at the largest estimate, however, is rather implausible (dartboard bottom left).

7 The joint process in case of change points

Let $\mathbf{X}\in\mathcal{M}$ with $\Gamma_{u}=I$ (symmetry). Proposition 7.1 below yields $(J_{h,t})_{t}\stackrel{{\scriptstyle d}}{{\approx}}(j_{h,t})_{t}+(\mathcal{L}_{h,t}^{*})_{t}$ approximately, with non-random $(j_{h,t})_{t}$ describing excursions from zero in the $h$-neighborhood of any $c_{u}$. $(\mathcal{L}_{h,t}^{*})_{t}$ is a Gaussian process with $\mathcal{L}_{h,t}^{*}\sim N_{2}(0,I)$.

Figure 5 extends Figure 1II and III: $J_{h,t}=(E_{h,t},V_{h,t})^{\mathsf{T}}$ fluctuates close to $j_{h,t}:=(e_{h,t},v_{h,t})^{\mathsf{T}}$ (red, green or light blue). The sensitivity of $J_{h,t}$ to the effects is explicated in $j_{h,t}$: $e_{h,t}$ is sensitive to changes in $\mu$ and robust against higher order changes, and $v_{h,t}$ reacts to changes in $\sigma^{2}$ but not to higher order changes. This supports Algorithm 6.1. Also, $v_{h,t}$ reveals the error caused by a change in $\mu$, see IIIC (red). In IIID (red) the bivariate approach overcomes the error: for any excursion, $d_{I}(j_{h,t})$ is maximal at $c_{u}$. We consider $C=\{c\}$ but note direct generalization to any $C$ with $\min_{u}|c_{u}-c_{u-1}|\geq 2h$.

Proposition 7.1**.**

Let $\mathbf{X}\in\mathcal{M}$ with $C=\{c\}$ and $h$ be a window size such that $c\in[h,T-h]$, and let $\rho_{1}=\rho_{2}=0$. In $(\mathcal{D}_{\mathbb{R}^{2}}[h,T-h],d_{SK})$ it holds as $n\to\infty$

[TABLE]

$(\mathcal{L}_{h,t}^{*})_{t}$ is a bivariate continuous $2h$-dependent process with $\mathcal{L}_{h,t}^{*}\sim N_{2}(0,I)$ for all $t$, and for $t\not\in(c-h,c+h)$ it equals $(\mathcal{L}_{h,t})_{t}$ from (8). We give $(j_{h,t}^{(n)})_{t}$ in (18), and $(\tilde{D}_{h,t})_{t}$ and $(\Delta_{h,t}^{(n)})_{t}$ in (19). The result extends Proposition 4.2. The marginal at $c$ is Proposition 6.2. We first extend Lemma 4.1.

Lemma 7.2**.**

Let $\mathbf{X}\in\mathcal{M}$ with $C=\{c\}$ and $h$ a window size such that $c\in[h,T-h]$. For $j\in\{\ell,r\}$ it holds in $(\mathcal{D}_{\mathbb{R}}[h,T-h],d_{\|\cdot\|})$ as $n\to\infty$ almost surely $(\hat{\mu}_{j})_{t}\to(\tilde{\mu}_{j})_{t}$, $(\hat{\sigma}_{j}^{2})_{t}\to(\tilde{\sigma}_{j}^{2})_{t}$ and $(\hat{\nu}_{j}^{2})_{t}\to(\tilde{\nu}_{j}^{2})_{t}$.

See (16) and (17) for the limits. They are non-random functions in $t$, in $(c-h,c+h)$ depending on parameters of both $\mathbf{X}_{1}$ and $\mathbf{X}_{2}$. The tilde marks dependence on $t$. For $j\in\{\ell,r\}$ set $(\tilde{\vartheta}_{j})_{t}$ as a placeholder for the limits, and $\vartheta_{1}$ and $\vartheta_{2}$ as the associated parameters. E.g., choose $\tilde{\mu}_{\ell}$ with $\mu_{1}$ and $\mu_{2}$, see Figure 6.

If a window does not lap $c$, then the RVs refer to a single $\mathbf{X}_{u}$ and $\tilde{\vartheta}_{j}$ equals $\vartheta_{u}$,

[TABLE]

e.g., $\tilde{\vartheta}_{r}=\vartheta_{2}$ in Figure 6, also see Lemma 4.1. If a window laps $c$, then $\tilde{\vartheta}_{j}$ depends on both $\mathbf{X}_{1}$ and $\mathbf{X}_{2}$, see left window in Figure 6. Set $\tilde{\vartheta}_{\ell}$ for $t\in[c,c+h)$

[TABLE]

and define $\tilde{\vartheta}_{r}$ for $t\in[c-h,c)$ by replacing in (17) all subscripts $\ell$ with $r$ and $t$ with $t+h$. We comment on the left window $(t-h,t]$. A proportion of $[c-(t-h)]/h$ RVs belongs to $\mathbf{X}_{1}$ and $(t-c)/h$ to $\mathbf{X}_{2}$. Then $\tilde{\mu}_{\ell}$ is a linear interpolation between $\mu_{1}$ and $\mu_{2}$. The higher orders are further affected by $\tilde{\mu}_{\ell}$ via the error $(\tilde{\mu}_{\ell}-\mu_{u})^{k}$ with $k\in\{1,2,4\}$. E.g., in $\tilde{\sigma}_{\ell}^{2}$ we find $(\tilde{\mu}_{\ell}-\mu_{u})^{2}$ deriving from $\hat{\sigma}_{\ell}^{2}=(1/nh)\sum_{i\in I_{\ell}}(X_{i}-\hat{\mu}_{\ell})^{2}$ as $\hat{\mu}_{\ell}$ relies on all RVs in $(t-h,t]$ irrespective of $c$. All errors vanish iff $\mu_{1}=\mu_{2}=\tilde{\mu}_{\ell}$.

We now define the centering $j_{h,t}^{(n)}=(e_{h,t}^{(n)},v_{h,t}^{(n)})^{\mathsf{T}}$ from $(\ref{tilde_vartheta})$ and $(\ref{tilde_left})$ via

[TABLE]

see $E_{h,t}^{(n)}$ and $V_{h,t}^{(n)}$ in (5). For $t\in(c-h,c+h)$ we discuss $e_{h,t}$ and $v_{h,t}$ w.r.t. the type of change. For $e_{h,t}$, see Figure 5B. Message 1: constant expectation implies $e_{h,t}=0$ (green). This is because $\mu_{1}=\mu_{2}$ implies $\tilde{\mu}_{\ell}=\tilde{\mu}_{r}$ see (17), regardless of higher moments. Thus, $E_{h,t}^{(n)}$ is robust against higher order changes. Message 2: a change in expectation causes $e_{h,t}$ to deviate from zero (red or light blue). Thus, $E_{h,t}^{(n)}$ is sensitive to expectation changes. $\tilde{\mu}_{j}$ linearly interpolates $\mu_{1}$ and $\mu_{2}$, thus $\tilde{\mu}_{r}-\tilde{\mu}_{\ell}$ has a hat shape peaking at $c$ and is additionally scaled by $(\tilde{\sigma}_{r}^{2}+\tilde{\sigma}_{\ell}^{2})^{1/2}$.

For $v_{h,t}$ see Figure 5C. Complexity rises as a change in expectation affects higher moments. Message 1: constant expectation results in $v_{h,t}$ behaving analogously to $e_{h,t}$: first, constant variance implies $v_{h,t}=0$, as for $\mu_{1}=\mu_{2}$ and $\sigma_{1}^{2}=\sigma_{2}^{2}$ we find $\tilde{\sigma}_{\ell}^{2}=\tilde{\sigma}_{r}^{2}$, see (17). Second, a change in variance causes a deviation from zero (green): if $\mu_{1}=\mu_{2}$, then $\tilde{\sigma}_{j}^{2}$ linearly interpolates $\sigma_{1}^{2}$ and $\sigma_{2}^{2}$, thus $\tilde{\sigma}_{r}^{2}-\tilde{\sigma}_{\ell}^{2}$ has a hat shape that is scaled with $(\tilde{\nu}_{r}^{2}+\tilde{\nu}_{\ell}^{2})^{1/2}$. Thus, given constant expectation, $V_{h,t}^{(n)}$ is sensitive to changes in variance and robust against higher order changes. Regarding $j_{h,t}$ note a linear trajectory along the ordinate, see Figure 5D (green). Message 2: a change in expectation falsely results in a systematic deviation of $v_{h,t}$ from zero even if the variance is constant (red). Thus, $V_{h,t}^{(n)}$ is also sensitive to changes in expectation.

We state the remaining functions of (15). We set $(\Delta_{h,t}^{(n)})_{t}$ and $(\tilde{D}_{h,t})_{t}$ as

[TABLE]

with $\Delta_{h,t,1}^{(n)}:=[(\tilde{\sigma}_{r}^{2}+\tilde{\sigma}_{\ell}^{2})/(\hat{\sigma}_{r}^{2}+\hat{\sigma}_{\ell}^{2})]^{1/2}$ and $\Delta_{h,t,2}^{(n)}:=[(\tilde{\nu}_{r}^{2}+\tilde{\nu}_{\ell}^{2})/(\hat{\nu}_{r}^{2}+\hat{\nu}_{\ell}^{2})]^{1/2}$ as well as $\tilde{D}_{h,t,1}:=\lim_{n\to\infty}[(\tilde{\sigma}_{r}^{2}+\tilde{\sigma}_{\ell}^{2})/(nh\mathbb{V}\!ar(\hat{\mu}_{r}-\hat{\mu}_{\ell}))]^{1/2}$ and $\tilde{D}_{h,t,2}:=\lim_{n\to\infty}[(\tilde{\nu}_{r}^{2}+\tilde{\nu}_{\ell}^{2})/(nh\mathbb{V}\!ar(\hat{\sigma}_{r}^{2}-\hat{\sigma}_{\ell}^{2}))]^{1/2}$.

$\Delta_{h,t}^{(n)}$ replaces the limits with the estimators. It is $(\Delta_{h,t}^{(n)})_{t}\to(I)_{t}$ a.s. componentwise as $n\to\infty$, see Lemma 7.2. $\Delta_{h,c}^{(n)}$ is (13). $\tilde{D}_{h,t}$ replaces the limits with the true scaling of the numerators of $J_{h,t}^{(n)}$ yielding asymptotic unit variance. $\tilde{D}_{h,t}$ captures inconsistent parameter estimation. Outside $(c-h,c+h)$ and at $c$ both windows in their entirety refer to a single $\mathbf{X}_{u}$, stating consistency $\tilde{D}_{h,t}=I$. In Figure 5IIID (red) $J_{h,t}$ is close to $j_{h,t}$ but slightly shifted, as we had to distort $j_{h,t}$ with $\tilde{D}_{h,t}$. Note that $nh\mathbb{V}\!ar(\hat{\mu}_{\ell})$ considers $c$ while $\hat{\sigma}_{\ell}^{2}$ does not: if $(t-h,t]\not\ni c$ then $nh\mathbb{V}\!ar(\hat{\mu}_{\ell})\to\sigma_{u}^{2}$ as $n\to\infty$. If $(t-h,t]\ni c$ then $nh\mathbb{V}\!ar(\hat{\mu}_{\ell})\to[(c-(t-h))/h]\sigma_{1}^{2}+[(t-c)/h]\sigma_{2}^{2}$, note the absence of the error $(\tilde{\mu}_{\ell}-\mu_{1})^{2}$ in contrast to $\tilde{\sigma}_{\ell}^{2}$ in (17).

8 Practical Performance

In Subsection 8.1 we extend methodology to multiple windows. This improves the detection of change points on multiple time scales and with different effects. In Subsection 8.2 we discuss the companion R-package jcp. In Subsection 8.3 we provide simulation studies for $\mathcal{H}_{0}:C=\emptyset$. We find the asymptotic significance level $\alpha$ of the test to be kept in many scenarios. In Subsection 8.4 we discuss simulations for $C\not=\emptyset$. We confirm reliable detection accuracy and adequate interpretation. To conclude, we present a real data example in Subsection 8.5.

8.1 Extension to multiple windows

For $\mathbf{X}\in\mathcal{M}$ let $H$ be a set of increasingly ordered windows $h_{1}<\cdots<h_{w}$, yielding multiple processes $\{(J_{h,t})_{t}\;|\;h\in H\}$. Smaller $h_{k}$ are sensitive to rapid changes, as they do no overlap adjacent change points and thus ensure an unbiased excursion of $(J_{h,t})_{t}$ which supports precise detection. Larger $h_{k}$ improve detection of small effects, as at a $c\in C$ a dartboard is shifted outwards with order $\sqrt{nh}$, see Proposition 6.2.

We extend the test: for $(\mathcal{D}_{\mathbb{R}^{2}}[h_{k},T-h_{k}],d_{SK})$, $k=1,\ldots,w$ let $(\times_{k=1}^{w}\mathcal{D}_{\mathbb{R}^{2}}[h_{k},T-h_{k}],d_{SK}^{w})$ denote the product space with distance $d_{SK}^{w}$ metrizing the product topology. Poposition 4.2 extends to

Corollary 8.1**.**

Let $\mathbf{X}\in\mathcal{M}$ with $C=\emptyset$ and $H=\{h_{1},\ldots,h_{w}\}$ be a set of windows. In $(\times_{k=1}^{w}\mathcal{D}_{\mathbb{R}^{2}}[h_{k},T-h_{k}],d_{SK}^{w})$ it holds as $n\to\infty$

[TABLE]

with $\Gamma$ in (7), and $(\mathcal{L}_{h_{k},t})_{t}$ as in (8) being evaluated on a single planar Brownian motion over all $k=1,\ldots,w$.

Under $C=\emptyset$ continuous mapping yields convergence of the global maximum

[TABLE]

for $n\to\infty$, which extends (9). We use $M(\Gamma):=M^{(1)}(\Gamma)$ as a test statistic and reject iff it exceeds the $(1-\alpha)$-quantile $Q$ of the limit distribution. The latter is approximated in simulations. We treat $\Gamma$ analogously to the case of single $h$, see Table 1: first, symmetry assumption yields $\Gamma=I$, and thus $d_{I}$. Second, consistent estimation $\hat{\Gamma}$ depends on $t$ and now also on $h$, see (10), yielding $d_{\hat{\Gamma}}$. Third, using $d_{\infty}$ avoids estimation and is conservative. For the process perspective choose $\mathscr{R}=\mathscr{R}(H)$ such that any $(J_{h,t})_{t}$ enters $\mathscr{R}$ only with probability $\alpha$, if $C=\emptyset$. Let $\mathscr{C}$ be the circle with center zero and radius $Q$. For the three cases we obtain first $\mathscr{R}=\mathscr{C}^{c}$, second $\mathscr{R}=\hat{\mathscr{E}}^{c}$ with $\hat{\mathscr{E}}=\hat{\Gamma}^{1/2}\cdot\mathscr{C}$ depending on $t$ and $h$, and third $\mathscr{R}=\mathscr{S}^{c}$ with $\mathscr{S}$ being the square with center zero and edge length $2Q$. We extend Theorem 5.1.

Theorem 8.2**.**

Let $\mathbf{X}\in\mathcal{M}$, and for $\alpha\in(0,1)$ let $Q$ be the $(1-\alpha)$-quantile of the limit in (20). Assume $\mathcal{H}_{0}:C=\emptyset$. Then it holds $\lim_{n\to\infty}\mathbb{P}(M^{(n)}(\Gamma)>Q)=\alpha$. Further, it holds first $\lim_{n\to\infty}\mathbb{P}(M^{(n)}(I)>Q)=\alpha$ given $\rho=0$, second $\lim_{n\to\infty}\mathbb{P}(M^{(n)}(\hat{\Gamma})>Q)=\alpha$, and third $\lim_{n\to\infty}\mathbb{P}(M^{(n)}(\infty)>Q)\leq\alpha$.

If $\mathcal{H}_{0}$ is rejected we aim at estimating $C$. We use a bottom-up approach similar to Messer et al., (2014). First detect candidates $\hat{C}_{h_{k}}$ for each $h_{k}$. Then merge the candidates into a final set $\hat{C}$ favoring smaller over larger windows.

Algorithm 8.3**.**

Set $\mathscr{R}(H)$ as either $\mathscr{C}^{c}$, $\hat{\mathscr{E}}^{c}$ or $\mathscr{S}^{c}$ as used in the test. For each $h_{k}\in H$ obtain $\hat{C}_{h_{k}}$ using Algorithm 6.1 w.r.t. $\mathscr{R}(H)$. Set $\hat{C}:=\hat{C}_{h_{1}}$. For increasing $k=2,3,\ldots w$ update $\hat{C}$ as follows: add any $\hat{c}\in\hat{C}_{h_{k}}$ to $\hat{C}$ that satisfies $\{\hat{c}-h_{k}+1,\ldots,\hat{c}+h_{k}\}\,\cap\,\hat{C}\,=\,\emptyset$.

A candidate from the $k$-th step is dismissed if a previously accepted estimate falls into the candidates $h_{k}$-neighborhood. See Figure 7 with $C=\{c_{1},c_{2},c_{3}\}$ and $|H|=4$ yielding three estimates $|\hat{C}|=3$. Note that $|\hat{C}_{k}|\leq 2$, i.e., no single $h_{k}$ estimates all three change points, as small $h_{k}$ lack power to detect the small effect at $c_{3}$, while larger $h_{k}$ are not sensitive to the close proximity of $c_{1}$ and $c_{2}$.

8.2 The jcp-package

The method is made available in the R-package jcp (joint change point detection) on CRAN. The package also contains a summary and a plotting routine. For the example of Figure 7 it generates Figure 8.

Panel A shows $\mathbf{X}$. $C=\emptyset$ is rejected. The bars mark $\hat{C}=\{420,500,743\}$. Panel B shows the means ($\hat{\mu}\approx 2.1,11.3,10.3,5.4$) and standard deviations ($\hat{\sigma}\approx 4.0,4.2,11.7,8.5$) calculated from all $(X_{i})_{i}$ in the estimated sections: we start at ’$+$’ (first section) and follow the arrows, i.e., point right (increase $\hat{\mu}$ at $420$), then upwards (increase in $\hat{\sigma}$ at $500$), and then south-east (decrease in both $\hat{\mu}$ and $\hat{\sigma}$ at $743$). The length of an arrow marks the effect size. Panel C facilitates interpretation, recall Section 6: the symmetry assumption yields $\mathscr{R}=\mathscr{C}^{c}$. The $(J_{h,t})_{t}$ are omitted for overview. The $h_{k}$ are color-coded (legend). The effect at $\hat{c}_{1}=420$ is strong as it is detected by the smallest $h_{1}=70$ (light green), also the dartboard lies far off. The dartboard touches the abscissa suggesting a change only in $\mu$. The effect at $\hat{c}_{2}=500$ is less strong. However, it was still $h_{1}$ that registered the change point. The dartboard lies close to the ordinate suggesting a change only in $\sigma^{2}$. The largest $h_{4}=160$ (dark red) was necessary to detect $\hat{c}_{3}=743$. The effect is the smallest among all (compare short arrow in B). The dartboard is positioned diagonally suggesting changes in both parameters.

8.3 Significance level and window choice

Significance level

We constructed an asymptotic test, where $\mathscr{R}\in\{\mathscr{C}^{c},\hat{\mathscr{E}}^{c}\}$ asymptotically keeps the $\alpha$-level, while $\mathscr{R}=\mathscr{S}^{c}$ reduces it, see Table 1 above. We fixed $h$, linked it to $n$ via $nh$, and let $n\to\infty$. In practice however, we set $n=1$. Thus, for a suitable approximation of $Q$ (resp. $\mathscr{R}$), we need the smallest window $h_{1}$ to be sufficiently large. Here we evaluate the rejection probability under $C=\emptyset$ for $n=1$ in simulations. It turns out that $h_{1}\approx 50$ typically keeps $\alpha$.

In the following let $\alpha=5\%$. Processes with $T=1000$ and $C=\emptyset$ (here i.i.d. RVs) are simulated. The relative frequency of rejections $f_{\mathscr{R}}$ ($1000$ simulations) approximates the rejection probability.

Figure 9 shows $f_{\mathscr{R}}$ as a function of the window set $H_{k}$, differentiating four distributions: $N(0,1)$ (magenta), $exp(1)$ (green), $gamma(s=1/2,\lambda=2)$ (blue) and $gamma(2,2)$ (red) with shape $s$ and rate $\lambda$. For $N(0,1)$ we use $\mathscr{R}=\mathscr{C}^{c}$, else $\mathscr{R}=\hat{\mathscr{E}}^{c}$ (solid lines) or $\mathscr{R}=\mathscr{S}^{c}$ (dotted lines). We choose $H_{k}=\{10(1+k),10(3+k),\ldots,10(11+k)\}$ for $k\in\{0,1,\ldots,15\}$, i.e., $|H_{k}|=6$ and all windows increase by $10$ when $k$ rises one unit. We see that $f_{\mathscr{R}}$ tends to the true $\alpha=5\%$ when increasing $k$, if $\mathscr{R}=\mathscr{C}^{c}$ or $\hat{\mathscr{E}}^{c}$. Also, $f_{\mathscr{R}}$ is reduced for $\mathscr{R}=\mathscr{S}^{c}$. Overall, the approximation is adequate for $h_{1}$ of about $30$ to $50$. Note that $f_{\mathscr{R}}$ is overestimated for $h_{1}$ small ($k\in\{0,1\}$): intuitively, consider $N(\mu,\sigma^{2})$ RVs, such that $E_{h,t}^{(n)}\sim t(v)$ and $t(v)$ is heavy-tailed (especially for $n$ small), while its limit is $N(0,1)$. Larger $h_{1}$ may also help decrease susceptibility to outliers.

Figure 10 shows the dependence of $f_{\mathscr{R}}$ on $\mu$ and $\sigma$. Again we see that $\alpha$ is kept for $h_{1}\approx 50$. We fix $H=\{50,75,100,125,150\}$ and consider $N(\mu,\sigma^{2})$ with $\mathscr{R}=\mathscr{C}^{c}$ (A), and $gamma(s,\lambda)$ with $\mathscr{R}=\hat{\mathscr{E}}^{c}$ (B) or $\mathscr{R}=\mathscr{S}^{c}$ (C), mentioning $s=\mu^{2}/\sigma^{2}$ and $\lambda=\mu/\sigma^{2}$. The legend in A shows the color-coding of $f_{\mathscr{R}}$: green indicates $5\%$, red an increase and blue a decrease.

For $N(\mu,\sigma^{2})$ we choose $\mu\in\{-3,-2,\ldots,3\}$ and $\sigma\in\{0.1,0.4,0.7,\ldots,2.2\}$ ($7\cdot 8=56$ combinations) and find $f_{\mathscr{R}}\approx 5\%$ for most combinations (entirely green). For $gamma(s,\lambda)$ we use $\mu,\sigma$ $\in\{0.1,0.3,0.5,\ldots,2.1\}$ ($11^{2}=121$ combinations). In B, $f_{\mathscr{R}}$ varies with $\mu$ and $\sigma$. For $\mu=\sigma$ (diagonal line, $s=1$, $exp(\lambda)$), it is $f_{\mathscr{R}}\approx 5\%$, which is decreased for $\mu<\sigma$, and increased for $\mu>\sigma$, but with $f_{\mathscr{R}}<10\%$ throughout. The figure suggests $f_{\mathscr{R}}$ to be constant for fixed $s$ and varying $\lambda$, i.e., along the lines $\mu=s\cdot\sigma$ (dotted lines, upper $s=3/4$, lower $s=3/2$), for which we mention rescaling in $J_{h,t}$ as well as the scaling property $r\cdot gamma(s,\lambda)=gamma(s,\lambda/r)$. Note that in B due to $\mathscr{R}=\hat{\mathscr{E}}^{c}$ the estimator $\hat{\Gamma}$ increases variability. The pattern in C is similar to B as the same distributions are considered, but overall it is more blue with $f_{\mathscr{R}}<3.7\%$ throughout, aligning with the conservative nature of $\mathscr{R}=\mathscr{S}^{c}$.

Window choice

We comment on the choice of $H$. For that we evaluate the dependence of $Q$ on $H$, recall (20). The radii in Figure 11 show simulations of $Q$ (also denoted $Q$, $10^{6}$ simulations).

Figure 11A shows three window sets for $T=1000$ fix: first $H_{1}=\{50\}$ yields $Q_{1}\approx 4.12$ (black), second $H_{2}=\{50,60,70,\ldots,150\}$ ($|H_{2}|=11$) gives $Q_{2}\approx 4.39$ (blue) and third $H_{3}=\{50,51,52,\ldots,150\}$ ($|H_{3}|=101$) yields $Q_{3}\approx 4.5$ (red). We mention monotonicity: $H_{A}\subset H_{B}$ implies $Q_{A}<Q_{B}$ which follows from global maximization over all $h\in H$, see (20). But we also observe a flattening in the increase, i.e., $Q_{2}-Q_{1}$ (black to blue) is larger than $Q_{3}-Q_{2}$ (blue to red) although there are much more $h$ added when going from blue to red. The reason is that in $(\ref{M})$, over all $h\in H$ the $(\mathcal{L}_{h,t})_{t}$ rely on a single Brownian motion such that $\mathcal{L}_{h_{a},t}\approx\mathcal{L}_{h_{b},t}$ for $h_{a}\approx h_{b}$ with equality if $h_{a}=h_{b}$. Thus, if $H$ is a fine grid then additional $h$ have minor impact on the global maximum. Figure 11B shows four choices of $T$ for fixed $H=\{50,60,\ldots,150\}$. $T_{1}=500$ yields $Q_{1}\approx 4.14$ (black), $T_{2}=1000$ results in $Q_{2}\approx 4.39$ (blue - same setup as in A), $T_{3}=2000$ yields $Q_{3}\approx 4.6$ (red) and $T_{4}=5000$ results in $Q_{4}\approx 4.83$ (magenta). Again, we mention monotonicity: $T_{A}<T_{B}$ implies $Q_{A}<Q_{B}$. This is because in (20) the $(\mathcal{L}_{h,t})_{t}$ are evaluated for $t\in[h,T-h]$ i.e., for larger $T$ the maxima are taken over longer intervals.

For the choice of $H$ we now argue that first, the richer $H$ the better a scenario of unknown change points and effects is exploited: for $c\in C$ we find an $h$ preferably large but with $(c-h,c+h]$ free from other change points, resulting in an unbiased excursion of $(J_{h,t})_{t}$, see Propositions 6.2 and 7.1. The order of the excursion is $\sqrt{h}$, while the competitor $Q$ is bounded (for fixed $T$ and $\alpha$) by the quantile associated with all possible windows. Second, the smallest window should be chosen large enough, e.g., $50$, for the $\alpha$-level of the test to be approximately kept.

8.4 Performance evaluation

For $C\not=\emptyset$ we evaluate the detection performance in simulations. We show precise estimation of the number and location of $C$, and appropriate interpretation of effects. First, we consider well-disposed setups using a single $h$, see Figure 12. Then we reduce effects and show improvement for multiple $H$, see Figure 13. Finally, we consider more complex processes, see Table 2. We set $\hat{\kappa}_{u}^{(d)}:=\{\hat{c}\in\hat{C}:|\hat{c}-c_{u}|\leq d\}$, describing ’correct’ estimates of $c_{u}$ at distance $\leq d$, for $d=10$ default, and $25$ or $5$. We perform $1000$ simulations with $T=1000$ and $\alpha=5\%$ throughout.

In Figure 12 and 13 we differentiate $N(\mu,\sigma^{2})$ with $\mathscr{R}=\mathscr{C}^{c}$ (left) and $gamma(s,\lambda)$ with $\mathscr{R}=\mathscr{S}^{c}$ (right). We consider piecewise i.i.d. RVs. It is $C=\{c_{1},c_{2},c_{3}\}$ with a change only in $\mu$ at $c_{1}$ (magenta cross), a change only in $\sigma^{2}$ at $c_{2}$ (blue cross), and change in both at $c_{3}$ (magenta and blue cross), see $\mathbf{X}$ in A. $1000$ simulations yield $3000$ change points. The $10$-neighborhood of any $c_{u}$ is accentuated by a gray box. $\hat{\kappa}_{u}^{(10)}$ are color-coded according to $c_{u}$ ($c_{1}=$ red, $c_{2}=$ green, $c_{3}=$ light blue, incorrect $\hat{C}\backslash\cup_{u=1}^{3}\hat{\kappa}_{u}^{(10)}=$ black).

In Figure 12 we set $H=\{100\}$, $C=\{250,500,750\}$, for $N(\mu,\sigma^{2})$ it is $\mu=2,10,10,2$ and $\sigma=4,4,16,4$, and for $gamma(s,\lambda)$ it is $\mu=0.8,2,2,4$ and $\sigma=1,1,0.1,2$. B shows $\hat{C}$: three narrow histograms around the true $c_{u}$ support precise estimation (frequencies reported in table), e.g., left $|\hat{C}|=3019$. Of thousand $c_{1}$ it is $|\hat{\kappa}_{1}^{(d)}|=998,(1000,988)$ for $d=10,(25,5)$, i.e., $|\hat{\kappa}_{1}^{(10)}|=998$ estimates yield the red, $|\hat{\kappa}_{2}^{(10)}|=948$ the green, $|\hat{\kappa}_{3}^{(10)}|=946$ the light blue histogram. It is $|\hat{C}\backslash\cup_{u=1}^{3}\hat{\kappa}_{u}^{(10)}|=127$ (black, low number hardly visible). Panel C shows $J_{h,\hat{c}}$ for all $\hat{c}\in\hat{C}$, e.g., left $998$ red points. The $J_{h,\hat{c}}$ distribute closely to the asymptotic distribution of $J_{h,c}$ for $c\in C$ (golden dartboards). Thus, interpretation of effects based on $J_{h,\hat{c}}$ is plausible: e.g., a typical red point indicates a change only in $\mu$. We see $127$ black points: those close to the brown rejection boundary represent real false positives resulting from chance. In contrast, those in the area of the dartboards refer to true $c_{u}$ but are classified incorrect due to $d=10$ (tails of the histograms in B).

Before, we showed strong performance under a well-disposed setup for a single $H=\{100\}$. We now show improvement for richer $H$. We choose closer distances of the first two change points $C=\{440,500,750\}$, see Figure 13, and also reduce effects for $N(\mu,\sigma^{2})$ to $\mu=2,10,10,6$ and $\sigma=4,4,12,10$. For $gamma(s,\lambda)$ parameters are as before. We consider three sets $H_{1}=\{100\}$ (B), $H_{2}=\{70,100,130,160\}$ (C) and $H_{3}=\{50,60,\ldots,200\}$ (D). Richer $H$ almost always improve the performance: first, $|\hat{C}|$ increases. Second, the $|\hat{\kappa}_{u}^{(d)}|$ increase. Third, B,C and D show improvement in location precision as the $\hat{c}$ lie closer to the true $c_{u}$. The small $h$ capture the rapid $c_{1}$ and $c_{2}$, and the larger $h$ improve detection of the small effect at $c_{3}$.

We further support the quality of the performance now considering five change points and different distributional assumptions, including processes with higher varying moments, see Table 2. We use $H=\{50,60,\ldots,200\}$ and set $C=\{200,260,500,720,810\}$ with $\mu=11,13,10,8,5,5$, $\sigma=1,3,3,3,4,1.3$, i.e., changes in $(\mu\;\&\;\sigma^{2})$, $\mu$, $\mu$, $(\mu\;\&\;\sigma^{2})$ and $\sigma^{2}$. We consider $N(\mu,\sigma^{2})$ with $\mathscr{R}=\mathscr{C}^{c}$, $gamma(s,\lambda)$ with $\mathscr{R}=\mathscr{E}^{c}$, $unif(a,b)$ as well as for $p=2,6,10$ the periodic processes as in Section 3 (mix of uniforms) with $\mathscr{R}=\mathscr{C}^{c}$. Indeed, the estimates $|\hat{C}|$ lie close to the true $5000$, and the $|\hat{\kappa}_{u}^{(d)}|$ often approach $1000$.

There are other prominent methods available on CRAN, e.g., mosum (Meier et al., , 2019), wbs (Baranowski and Fryzlewicz, , 2019), changepoint (Killick et al., , 2016), stepR (Pein et al., , 2020), cumSeg (Muggeo, , 2020) or FDRSeg (Li and Sieling, , 2017). But none of them captures all aspects covered by jcp, i.e., changes in both $\mu$ and $\sigma^{2}$, multiple time scales, and nonparametric methodology allowing piecewise different distributions, e.g., mosum essentially uses $E_{h,t}$ to address $\mu$ with a single $h$.

8.5 Data example

We analyze the frequency of the nucleobase uracil in the genome of the severe acute respiratory syndrome coronavirus 2 isolate Wuhan-Hu-1 (SARS-CoV-2), Wu et al., (2020, GenBank: MN908947.3). The genome contains a sequence of $29903$ bases, which we decompose into $T=996$ subsequent sections of length $30$, in each of which we compute the frequency of uracil, see Figure 14A.

We set $\alpha=0.05$, $\mathscr{R}=\mathscr{S}^{c}$ and $H=\{50,70,90,110,130\}$. $\mathcal{H}_{0}$ was rejected. Three change points were detected $\hat{C}=\{219,391,942\}$. Figure 14C indicates a decrease in $\mu$ at $391$ and $942$, and an increase in $\mu$ and also slightly in $\sigma^{2}$ at $219$. The periodicity captured in the model (see Figure 2) helps mimic the data as compared to i.i.d. sequences (see e.g. Figure 1). None of the four segments shows serial correlation. Other decomposition lengths (e.g., $20$ or $40$) yield analogous results.

9 Discussion

We proposed a method for change point detection in univariate sequences. Data are modeled as independent RVs with piecewise constant $\mu$ and $\sigma^{2}$, apart from which the distribution is allowed to vary periodically. In order to jointly detect changes in $\mu$ and $\sigma^{2}$, we developed a bivariate MOSUM approach: a comparison of means in adjacent windows is sensitive to changes in $\mu$, and a comparison of empirical variances addresses $\sigma^{2}$. Methodologically, first an asymptotic test was constructed to test the null hypothesis that no change occurred, $\mathcal{H}_{0}:C=\emptyset$. Second, an algorithm for change point detection was presented that can be run if $\mathcal{H}_{0}$ is rejected.

The method is grounded on the asymptotic behavior of the MOSUM. Under $\mathcal{H}_{0}$ it was shown to approach a zero mean Gaussian process from which the rejection boundary of the test was derived in simulations. Under the alternative it asymptotically describes a Gaussian process that systematically deviates from zero locally around a change point. The quantification of this deviation supported the sensitivity of the MOSUM process to change points. The first component was shown to be sensitive to changes in $\mu$ and robust against higher order changes. The second component was shown to be sensitive to changes in $\sigma^{2}$, but unfortunately it also reacts to changes in $\mu$. Therefore, it is the joint consideration of both components that supports adequate estimation. Inference on the strength of a change (effect size) as well as its type ($\mu$, $\sigma^{2}$, or both) is enabled with confidence.

The bivariate MOSUM was classified according to the presence or absence of symmetry (more precisely the averaged third centered moment) of the underlying distributions. Symmetry, on the one hand, was shown to result in independent asymptotic components and is thus related to Euklid’s notion of distance: under $\mathcal{H}_{0}$, symmetry results first in the isotropy of the MOSUM process and thus in a rejection area given by a circle. Second, in case of change points, we found the marginals at the change points to be asymptotically uncorrelated, depicted by round contour lines. On the other hand, a lack of symmetry of the underlying distribution was shown to result in correlated components of the MOSUM process, and is therefore related to Mahalanobis’ distance: first, under $\mathcal{H}_{0}$ the MOSUM has a preferred direction of variability along the main diagonals in $\mathbb{R}^{2}$, which results in an elliptic rejection boundary. Second, the marginals at the change points have correlated components, resulting in elliptic contour lines. We presented three ways to treat the unknown correlation in applications: first, symmetry assumption results in vanishing correlation. Second, correlation can be estimated consistently. Third, a conservative approach that avoids to address correlation was proposed.

The method was further extended to improve the detection of change points on different scales. For that, multiple bivariate MOSUM processes were applied simultaneously. Indeed, various simulation studies revealed strong performance under different distributional assumptions, including piecewise i.i.d. sequences and processes with varying higher moments. Generally, weak distributional assumptions allow for a wide range of applications. The method is implemented in the R package jcp, which performs the test and the algorithm, summarizes the results and provides a graphical output to facilitate interpretation.

Acknowledgements

The author is very grateful for valuable comments by Götz Kersting, Brooks Ferebee, Ralph Neininger and Anja Nowak, and for helpful suggestions of three anonymous referees.

Contact information

Michael Messer

Vienna University of Technology,

Institute of Statistics and Mathematical Methods in Economics,

Wiedner Hauptstraße 8-10/105, 1040 Vienna, Austria

tel.: +43 -1 -58801 -10588, email.: michael.messer@tuwien.ac.at

Appendix

We give here all proofs and additional auxiliary results.

ad Section 4

First we state a functional strong law of large numbers (SLLN) for the auxiliary processes, see Definition 3.1.

Lemma 9.1**.**

For an auxiliary process $\textbf{X}=(X_{i})_{i=1,2,\ldots}$ of period $p$ and moments $\mu^{\langle k\rangle}=(1/p)\sum_{m=1}^{p}\mu(m)^{\langle k\rangle}$ it is for $k\in\{1,\ldots,4\}$ in $(\mathcal{D}_{\mathbb{R}}[0,T],d_{\|\cdot\|})$ a.s. as $n\to\infty$

[TABLE]

Proof: The consecutive blocks $(X_{lp+1},\ldots,X_{(l+1)p})_{l=0,1,\ldots}$ of length $p$ constitute an i.i.d. sequence. Define the mean $Z_{l}^{\langle k\rangle}:=(1/p)\sum_{m=1}^{p}X_{lp+m}^{k}$. Then, for $n\to\infty$ a.s.

[TABLE]

The SLLN applied to $(Z_{l}^{\langle k\rangle})_{l=0,1,\ldots}$ shows the first summand to tend to $\mu^{\langle k\rangle}$. The second summand is the remainder of the last block and vanishes as the moment assumptions imply $(1/n)\sum_{m=1}^{p}|X_{lp+m}|^{k}\to 0$ a.s. as $n\to\infty$. From (22) we obtain $\sup_{0\leq t\leq T}|(1/n)\sum_{i=1}^{\lfloor nt\rfloor}X_{i}^{k}-t\mu^{\langle k\rangle}|\to 0$ by discretizing time. ∎

Proof of Lemma 4.1: We aim to apply Lemma 9.1, for which we express the centered population parameters from (2) through the non-centered moments as

[TABLE]

and analogously $\mu^{\{3\}}=\mu^{\langle 3\rangle}-3\mu\mu^{\langle 2\rangle}+\mu^{3}$ and $\mu^{\{4\}}=\mu^{\langle 4\rangle}-4\mu\mu^{\langle 3\rangle}+6\mu^{2}\mu^{\langle 2\rangle}-3\mu^{4}$. Note that $\mu(m)=\mu$ constant implies $(1/p)\sum_{m=1}^{p}\mu(m)^{2}=\mu^{2}$. This is favorable, as Lemma 9.1 only states estimation of $\mu^{\langle k\rangle}=(1/p)\sum_{m=1}^{p}\mu(m)^{\langle k\rangle}=(1/p)\sum_{m=1}^{p}\mathbb{E}[X_{m}^{k}]$, but not of $(1/p)\sum_{m=1}^{p}\mu(m)^{k}=(1/p)\sum_{m=1}^{p}\mathbb{E}[X_{m}]^{k}$ if $\mu(m)$ is not constant. Analogously, the average

[TABLE]

from (3) behaves well as $\sigma(m)^{2}=\sigma^{2}$ constant yields $(1/p)\sum_{m=1}^{p}(\sigma(m)^{2})^{2}=\sigma^{4}$. If $\sigma(m)^{2}$ was not constant, $\hat{\sigma}_{j}^{4}$ would fail in estimating $(1/p)\sum_{m=1}^{p}(\sigma(m)^{2})^{2}$.

We consider $j=\ell$. Note that $\hat{\mu}_{\ell}^{\{k\}}$ from (6) can be expressed through $\hat{\mu}_{\ell}^{\langle k\rangle}$ as $\hat{\mu}_{\ell}^{\{2\}}=\hat{\mu}_{\ell}^{\langle 2\rangle}-\hat{\mu}_{\ell}^{2}$, $\hat{\mu}_{\ell}^{\{3\}}=\hat{\mu}_{\ell}^{\langle 3\rangle}-3\hat{\mu}_{\ell}\hat{\mu}_{\ell}^{\langle 2\rangle}+\hat{\mu}_{\ell}^{3}$ and $\hat{\mu}_{\ell}^{\{4\}}=\hat{\mu}_{\ell}^{\langle 4\rangle}-4\hat{\mu}_{\ell}\hat{\mu}_{\ell}^{\langle 3\rangle}+6\hat{\mu}_{\ell}^{2}\hat{\mu}_{\ell}^{\langle 2\rangle}-3\hat{\mu}_{\ell}^{4}$. A comparison with (23) reveals that it suffices to have $(\hat{\mu}_{\ell}^{\langle k\rangle})_{t}\to(\mu^{\langle k\rangle})_{t}$ a.s. in $(\mathcal{D}_{\mathbb{R}}[h,T-h],d_{\|\cdot\|})$, which we obtain by Lemma 9.1 as for $n\to\infty$

[TABLE]

In order to prepare the proof of Proposition 4.2, we state a bivariate functional central limit theorem (CLT) for the first and second moment of an auxiliary process.

Lemma 9.2**.**

Let $\textbf{X}=(X_{i})_{i=1,2,\ldots}$ be an auxiliary process with population parameters $\mu$, $\sigma^{2}$, $\nu^{2}$ and $\rho$. Rescale variables as $Y_{i}:=(X_{i}-\mu)/\sigma$ and $\mathcal{Y}_{i}:=[(X_{i}-\mu)^{2}-\sigma^{2}]/\nu$ and let $\Gamma$ be as in (7). Then it holds in $(\mathcal{D}_{\mathbb{R}^{2}}[0,\infty),d_{SK})$ as $n\to\infty$

[TABLE]

while $(W_{t},\mathcal{W}_{t})^{\mathsf{T}}$ constitutes a planar Brownian motion.

Proof: The blocks $(X_{lp+1},\ldots,X_{(l+1)p})_{l=0,1,\ldots}$ are i.i.d. Set $Z_{l}:=(1/\sqrt{p})\sum_{m=1}^{p}Y_{lp+m}$ and $\mathcal{Z}_{l}:=(1/\sqrt{p})\sum_{m=1}^{p}\mathcal{Y}_{lp+m}$. Then $(Z_{l},\mathcal{Z}_{l})_{l=0,1,\ldots}$ is an i.i.d. sequence with components having zero mean and unit variance. Regarding (24) we rewrite

[TABLE]

using $\lfloor nt/p\rfloor=\lfloor\lfloor nt\rfloor/p\rfloor$.

The second summand vanishes as $(1/\sqrt{n})\sum_{m=1}^{p}|Y_{lp+m}|\to 0$ in probability as $n\to\infty$. We treat the second component of (24) analogously. As $\mathbb{C}ov((Z_{1},\mathcal{Z}_{1})^{\mathsf{T}})=\Gamma$, the result follows from a functional CLT w.r.t. $(Z_{l},\mathcal{Z}_{l})_{l=0,1,\ldots}$, see Kuelbs, (1973) working on general Banach space valued RVs. ∎

Proof of Proposition 4.2: We define a continuous map $\varphi=(\varphi_{1},\varphi_{2})^{\mathsf{T}}$ from $(\mathcal{D}_{\mathbb{R}^{2}}[0,T],d_{SK})$ to $(\mathcal{D}_{\mathbb{R}^{2}}[h,T-h],d_{SK})$, for both components $\kappa=1,2$ via

[TABLE]

Applying $\varphi$ on (24) yields in $(\mathcal{D}_{\mathbb{R}^{2}}[h,T-h],d_{SK})$ as $n\to\infty$

[TABLE]

The centering $\mu$ and $\sigma^{2}$ canceled. In the second component replace $\mu$ with $\hat{\mu}_{j}$ yielding

[TABLE]

Regarding the replacement we note that for $j\in\{\ell,r\}$

[TABLE]

while the second summand vanishes uniformly over $t$ a.s. as $n\to\infty$, as first Lemma 4.1 states $(\hat{\mu}_{j})_{t}\to(\mu)_{t}$ a.s., and second square-root scaling in $[(nh)^{-1/2}\sum_{i\in I_{j}}(X_{i}-\mu)]_{t}$ yields weak convergence to a Gaussian process limit. Finally, Lemma 4.1 also implies that $2\sigma^{2}$ and $2\nu^{2}$ in (27) can be replaced by $\hat{\sigma}_{r}^{2}+\hat{\sigma}_{\ell}^{2}$ and $\hat{\nu}_{r}^{2}+\hat{\nu}_{\ell}^{2}$. ∎

ad Section 5

Proof of Theorem 5.1: The result for $M_{h}^{(n)}(\Gamma)$ follows from Proposition 4.2 applying continuous mapping. From this, the statements for $M_{h}^{(n)}(I)$ and $M_{h}^{(n)}(\infty)$ follow by construction. For $M_{h}^{(n)}(\hat{\Gamma})$ we use consistency of $\hat{\rho}$, see Lemma 4.1. ∎

ad Section 6

Proof of Proposition 6.2: We write $J_{h,c}^{(n)}-\Delta_{h,c}^{(n)}\cdot j_{h,c}^{(n)}=\Delta_{h,c}^{(n)}\cdot D_{h,c}^{(n)}$ with

[TABLE]

We need to show $D_{h,c}^{(n)}\stackrel{{\scriptstyle d}}{{\longrightarrow}}N_{2}(0,\Gamma_{c})$ as $n\to\infty$, as $\Delta_{h,c}^{(n)}\to I$ a.s. by components, see Lemma 4.1. Set $\hat{\theta}_{h,c}^{(n)}:=(\hat{\mu}_{r},\hat{\sigma}_{r}^{2},\hat{\mu}_{\ell},\hat{\sigma}_{\ell}^{2})^{\mathsf{T}}$ and $\theta:=(\mu_{2},\sigma_{2}^{2},\mu_{1},\sigma_{1}^{2})^{\mathsf{T}}$. As $n\to\infty$

[TABLE]

using the Lindeberg-Feller CLT in each component while periodicity of the RVs implies the Lindeberg condition. Cramér-Wold device yields joint convergence. Set

[TABLE]

Then, the continuous mapping theorem yields as $n\to\infty$

[TABLE]

and further

[TABLE]

since $\Gamma_{c}=D\Sigma^{*}D^{\mathsf{T}}$, while $D$ denotes the diagonal matrix of the latter display. ∎

ad Section 7

Proof of Lemma 7.2: W.l.o.g. let $j=\ell$. Lemma 4.1 implies pointwise a.s. convergence to the population parameter if $(t-h,h]\not\ni c$. Else, a proportion of $[c-(t-h)]/h$ RVs belongs to $\mathbf{X}_{1}$, and $(t-c)/h$ to $\mathbf{X}_{2}$. Regarding $\hat{\mu}_{\ell}$ we obtain a.s. as $n\to\infty$

[TABLE]

Uniform convergence follows from Lemma 9.1. Analogously, for $\hat{\sigma}_{\ell}^{2}$ and $\hat{\nu}_{\ell}^{2}$ it remains to consider $(t-h,h]\ni c$. We replace $X_{u,i}$ in (29) and discuss the left subinterval $(t-h,c]$. For $\hat{\sigma}_{\ell}^{2}$ we find a.s. as $n\to\infty$

[TABLE]

For that, decompose $(X_{1,i}-\hat{\mu}_{\ell})^{2}=(X_{1,i}-\mu_{1})^{2}-2(X_{1,i}-\mu_{1})(\hat{\mu}_{\ell}-\mu_{1})+(\hat{\mu}_{\ell}-\mu_{1})^{2}$. The first summand yields $\sigma_{1}^{2}$, the second vanishes, and the third summand yields $(\tilde{\mu}_{\ell}-\mu_{1})^{2}$. Use here $\hat{\mu}_{\ell}\to\tilde{\mu}_{\ell}$ a.s. For $\hat{\nu}_{\ell}^{2}=\hat{\mu}_{\ell}^{\{4\}}-\hat{\sigma}_{\ell}^{4}$ it remains to comment on $\hat{\mu}_{\ell}^{\{4\}}$. We find a.s. as $n\to\infty$

[TABLE]

For this we decompose $(X_{1,i}-\hat{\mu}_{\ell})^{4}=(X_{1,i}-\mu_{1})^{4}-4(X_{1,i}-\mu_{1})^{3}(\hat{\mu}_{\ell}-\mu_{1})+6(X_{1,i}-\mu_{1})^{2}(\hat{\mu}_{\ell}-\mu_{1})^{2}-4(X_{1,i}-\mu_{1})(\hat{\mu}_{\ell}-\mu_{1})^{3}+(\hat{\mu}_{\ell}-\mu_{1})^{4},$ using $\hat{\mu}_{\ell}\to\tilde{\mu}_{\ell}$ and $\hat{\sigma}_{\ell}^{2}\to\tilde{\sigma}_{\ell}^{2}$ a.s. The fourth summand vanishes. ∎

Proof of Proposition 7.1: $\mathbf{X}$ relies on independent $\mathbf{X}_{1}$ and $\mathbf{X}_{2}$. Let $u\in\{1,2\}$. $\mathbf{X}_{u}$ has population parameters $\mu_{u}$, $\sigma_{u}^{2}$ and $\nu_{u}^{2}$. We standardize $Y_{u,i}:=(X_{u,i}-\mu_{u})/\sigma_{u}$ and $\mathcal{Y}_{u,i}:=[(X_{u,i}-\mu_{u})^{2}-\sigma_{u}^{2}]/\nu_{u}$. The functional CLT from Lemma 9.2 states that for $u\in\{1,2\}$ it holds in $(\mathcal{D}_{\mathbb{R}^{2}}[0,T],d_{SK})$ as $n\to\infty$

[TABLE]

while $(W_{u,t},\mathcal{W}_{u,t})_{t}$ is a planar Brownian motion. It is $\Gamma_{u}=I$ as $\rho_{u}=0$. The independence of $\mathbf{X}_{1}$ and $\mathbf{X}_{2}$ inherits to $(W_{1,t},\mathcal{W}_{1,t})_{t}$ and $(W_{2,t},\mathcal{W}_{2,t})_{t}$. It also implies joint convergence (30) over $u=1$ and $2$. We replace $\sigma_{u}$ and $\nu_{u}$ with $[h(\tilde{\sigma}_{r}^{2}+\tilde{\sigma}_{\ell}^{2})]^{1/2}$ and $[h(\tilde{\nu}_{r}^{2}+\tilde{\nu}_{\ell}^{2})]^{1/2}$, see (16) and (17), yielding in $(\mathcal{D}_{\mathbb{R}^{2}}[h,T-h],d_{SK})$ as $n\to\infty$

[TABLE]

To switch to the moving sum perspective, we define a continuous map $\varphi:(\mathcal{D}_{\mathbb{R}^{2}}[h,T-h]\times\mathcal{D}_{\mathbb{R}^{2}}[h,T-h],d_{SK}\otimes d_{SK})\to(\mathcal{D}_{\mathbb{R}}[h,T-h]\times\mathcal{D}_{\mathbb{R}}[h,T-h],d_{SK}\otimes d_{SK}).$ For that we write an element from the domain of $\varphi$ as $\iota:=[((f_{1},f_{2})(t))_{t},((g_{1},g_{2})(t))_{t}]$, with $f_{\kappa}$ associated with the $\kappa$-th component of the first process $(u=1)$, and $g_{\kappa}$ with the $\kappa$-th component of the second process $(u=2)$, for $\kappa=1,2$. Define $\varphi=(\varphi_{1},\varphi_{2})$ componentwise for $\kappa=1,2$ identical via

[TABLE]

We apply $\varphi_{\kappa}$ on (31). Continuous mapping preserves convergence. The first ($t<c-h$) and the fourth ($t\geq c+h$) summand in $\varphi_{\kappa}$ refer to a single $\mathbf{X}_{u}$ for which we obtain $(J_{h,t}^{(n)})_{t}\stackrel{{\scriptstyle d}}{{\longrightarrow}}(\mathcal{L}_{h,t})_{t}$ as in Proposition 4.2. We need to discuss $t\in(c-h,c+h]$. We focus on the case that $c$ lies in the right window $t\in[c-h,c)$, and consider the first component $\kappa=1$. Application of $\varphi_{1}$ on the left hand side of (31) yields

[TABLE]

For the right hand side of (31) we obtain

[TABLE]

which has zero expectation. Replacing $W_{1,t}$ and $W_{2,t}$ with a single Brownian motion $(W_{t})_{t}$ yields continuity. We replace the scaling: first, multiplication of (32) with the first entry of $\Delta_{h,t}^{(n)}$ preserves the limit as $(\Delta_{h,t}^{(n)})_{t}\to(I)_{t}$ a.s. Second, multiplication of both (32) and (33) with the first entry of $(\tilde{D}_{h,t})_{t}$ scales (33) to unit variance. Set the latter as the first component of $(\mathcal{L}_{h,t}^{*})_{t}$. In total, the first component of (15) holds true. Analogously, derive the second component. It is $\mathcal{L}_{h,t}^{*}\sim N_{2}(0,I)$. ∎

ad Section 8

Proof of Corollary 8.1: Set $\varphi^{h}:=\varphi$ as in (26) and proceed analogously to the proof of Proposition 4.2 by applying the joint map $(\varphi^{h_{1}},\ldots,\varphi^{h_{w}})$ on (24). ∎

Proof of Theorem 8.2: We argue analogously to the proof of Theorem 5.1, applying Corollary 8.1. ∎