Detection of collective and point anomalies at the presence of trend and seasonality

TL;DR

This paper introduces a new method for detecting both collective and point anomalies in time series data that accounts for trends and seasonality, with proven statistical accuracy and practical effectiveness.

Contribution

It presents a novel approach that simultaneously detects various anomaly types in complex time series, overcoming limitations of existing methods.

Findings

Accurately decomposes time series into anomaly, trend, seasonality, and residual components.

Consistently estimates the number and locations of anomalies with minimal error.

Demonstrates strong detection performance in simulations and real energy price data.

Abstract

Detecting anomalies in time series data is a challenging task with broad relevance in many applications. Existing methods work effectively only under idealized conditions, typically focusing on point anomalies or assuming a constant baseline. Our approach overcomes these limitations by detecting both collective and point anomalies, while allowing for polynomial trends and seasonal patterns. We establish statistical theory demonstrating that our method accurately decomposes the time series into anomaly, trend, seasonality, and a remainder component. We further show that it estimates the number of anomalies consistently and their locations with minimal error. Simulation studies confirm its strong detection performance with finite samples, and an application to energy price data illustrates its practical utility. An R package is available on request.

Tables

Table 1. Table 1 : Average anomaly detection accuracy of STAD compared to three oracle version of it. We vary the length and magnitude of a single anomaly. STAD performs nearly as good as the oracles unless the anomaly is longer than a quarter of the total time series.

Length (%)	Magnitude	STAD	Oracle Trend	Oracle Season	Oracle
5 (0.1%)	3.24	0.903	0.907	0.909	0.918
28 (0.56%)	1.42	0.905	0.907	0.919	0.924
51 (1.02%)	1.25	0.916	0.921	0.923	0.926
108 (2.16%)	1.24	0.921	0.920	0.926	0.923
278 (5.56%)	1.21	0.917	0.915	0.925	0.924
768 (15.36%)	1.20	0.929	0.934	0.941	0.938
1178 (23.56%)	1.19	0.909	0.917	0.920	0.926
1278 (25.56%)	1.19	0.903	0.916	0.905	0.920
1378 (27.56%)	1.19	0.887	0.920	0.886	0.926
1499 (29.98%)	1.19	0.770	0.914	0.761	0.912
1678 (33.56%)	1.18	0.565	0.920	0.566	0.926

Table 2. Table 2 : We consider a single anomaly of length 768 and vary the seasonal period P P . STAD performs almost as good as the oracles unless less than five seasonal cycles P / n P/n are present.

Seasonal Period ($P$)	STAD	Oracle Trend	Oracle Season	Oracle
50	0.936	0.937	0.941	0.938
100	0.933	0.930	0.941	0.938
250	0.933	0.934	0.941	0.938
500	0.931	0.932	0.941	0.938
800	0.924	0.911	0.941	0.938
1000	0.890	0.928	0.941	0.938
1500	0.023	0.008	0.941	0.938

Table 3. Table 3 : Average anomaly detection accuracy as a function of J J . Results improve with increasing J J . Results support our default of J = 20 J=20 as the performance increase only marginally afterwards.

$J$	STAD	Oracle Trend	Oracle Season	Oracle
5	0.701	0.920	0.715	0.926
10	0.832	0.920	0.839	0.926
20	0.887	0.920	0.886	0.926
30	0.896	0.920	0.908	0.926
40	0.909	0.920	0.915	0.926
50	0.904	0.920	0.910	0.926

Table 4. Table 4 : Average detection accuracy in a setting with multiple anomalies. We consider different anomaly lengths ( μ L \mu_{L} ) and magnitudes ( a , b ) (a,b) . Results are consistent with the single anomaly scenario in Section  4.1 . STAD performs very well unless anomalies are long, in this case μ L = 800 \mu_{L}=800 .

${\mu }_L$	$(a,b)$	STAD	Oracle Trend	Oracle Season	Oracle
30	(1, 3)	0.734	0.739	0.746	0.754
100	(1, 3)	0.845	0.845	0.847	0.852
400	(1, 3)	0.850	0.854	0.849	0.858
800	(1, 3)	0.646	0.872	0.647	0.871
30	(2, 5)	0.904	0.905	0.904	0.905
100	(2, 5)	0.909	0.909	0.909	0.909
400	(2, 5)	0.896	0.907	0.905	0.910
800	(2, 5)	0.655	0.894	0.676	0.911

Equations

$$Y=A+B+\epsilon ,$$

$$Y=A+S+T+\epsilon .$$

$$Y=A+T+S+\epsilon ,$$

$$a_i=a_i^{\mathrm{coll}}+a_i^{\mathrm{point}}.$$

$$a_i^{\mathrm{coll}}=\{\begin{array}{ll}{\mu }_k,& \text{if }{b}_k<i\le e_k,\\ 0,& \text{otherwise},\end{array}$$

$$T=X\beta ,$$

$$x_{i,j}={\left(\frac{i}{n}\right)}^{j-1},$$

$$s_i={\vartheta }_{(i\mathrm{mod}P)},$$

$$\sum _{i=1}^n{s}_i=0,$$

$${\epsilon }_i\sim \mathcal{N}(0,{\sigma }_0^2).$$

$${\hat{\sigma }}_0=\frac{1}{\sqrt{2}}\cdot \frac{\text{IQR}(\Delta y)}{{\text{IQR}}_{\Phi }},$$

$$Y=A+R+W+\epsilon .$$

$$C_{i,j}(Z):-\sum _{l=i}^j{z}_l^2,1\le i\le j\le n.$$

$$C(Z,\tilde{\mathcal{A}},\tilde{\mu }):-\sum _{k=0}^{\tilde{K}}{C}_{{\tilde{e}}_k+1,{\tilde{b}}_{k+1}}(Z)+\sum _{k=1}^{\tilde{K}}{C}_{{\tilde{b}}_k+1,{\tilde{e}}_k}(Z-{\tilde{\mu }}_k)+\tilde{K}{\lambda }_{\mathrm{coll}}{\hat{\sigma }}_0^2\log (n),$$

$$\begin{split}C_{O}(Z,\tilde{\mathcal{A}},\tilde{\mu})\coloneq&\sum_{\scriptscriptstyle i\notin\tilde{\mathcal{A}}\cup\tilde{\mathcal{O}}}z_{i}^{2}+\sum_{\scriptscriptstyle k=1}^{\scriptscriptstyle\tilde{K}}\Big{[}C_{\tilde{b}_{k}+1,\tilde{e}_{k}}(Z-\tilde{\mu}_{k})+\lambda_{\operatorname{coll}}\hat{\sigma}_{0}^{2}\log(n)\Big{]}+\sum_{\scriptscriptstyle i\in\tilde{\mathcal{O}}}\big{[}\lambda_{\operatorname{point}}\hat{\sigma}_{0}^{2}\log(n)\big{]},\end{split}$$

$$(\hat{\mathcal{A}},\hat{\mu }):-\underset{(\tilde{\mathcal{A}},\tilde{\mu })}{\mathrm{argmin}}C(Z,\tilde{\mathcal{A}},\tilde{\mu }),$$

$${\hat{a}}_i=\{\begin{array}{ll}{\hat{\mu }}_k& \text{if }\exists k:{\hat{b}}_k<i\le {\hat{e}}_k,\\ y_i& \text{if }i\in \hat{\mathcal{O}},\\ 0& \text{otherwise.}\end{array}$$

$$Y=T+A+W+\epsilon .$$

$$\hat{\beta }(U):-\underset{\beta \in {\mathbb{R}}^{Q+1}}{\mathrm{argmin}}\sum _{i\in U}\rho \left[\frac{z_i-{\sum }_{q=0}^Q{\beta }_q{x}_{i,q}}{\widehat{{\sigma }_0}}\right],$$

$${\rho }_c(x)=\{\begin{array}{ll}\frac{c^2}{6}\left[1-{\left(1-{\left(\frac{x}{c}\right)}^2\right)}^3\right],& \text{if }|x|\le c,\\ \frac{c^2}{6},& \text{if }|x|>c.\end{array}$$

$$V=B\times 2(Q+1)$$

$$\tilde{u}_{k}\coloneq\left\{\Big{\lfloor}(k-1)\frac{n}{V}\Big{\rceil}+1,\ldots,\Big{\lfloor}k\frac{n}{V}\Big{\rceil}\right\},\quad k=1,\ldots,V.$$

$$U:-\bigcup _{k=0}^Q{\tilde{u}}_{(2k+l)B+m_k},$$

$$\mathcal{U}\coloneq\big{\{}U_{\tilde{B},\tilde{Q},j}\subset\{1,\ldots,n\}:\tilde{B}\in\mathbf{B},\,\tilde{Q}\in\mathbf{Q},\,j=1,\ldots,J\big{\}},$$

$$\hat{\mathcal{T}}=\big{\{}\hat{T}_{U}:U\in\mathcal{U}\big{\}},$$

$$\hat{T}=\underset{T\in \hat{\mathcal{T}}}{\mathrm{argmin}}\underset{(\tilde{\mathcal{A}},\tilde{\mu })}{\min }C(Y-T,\tilde{\mathcal{A}},\tilde{\mu })+\tilde{Q}{\hat{\sigma }}_0\log (n).$$

$$y_{D,i}=y_{i+D}-y_i,\text{for}i=1,\dots ,n-D,$$

$$T_D=X_D{\mathbf{\beta }}_D,$$

$$(X_D{)}_{i,q}:-X_{i+D,q+1}-X_{i,q+1}={\left(\frac{i+D}{n}\right)}^q-{\left(\frac{i}{n}\right)}^q,$$

$${\hat{\beta }}_D=({\hat{\beta }}_{D,1},\dots ,{\hat{\beta }}_{D,Q})$$

Full text

1 Introduction

The ability to accurately and efficiently estimate the location of anomalies in time series data has become increasingly important in various domains, including energy (Serrano-Guerrero et al., 2020; Zhang et al., 2021), network monitoring (Wu et al., 2021) and economics (Zamanzadeh Darban et al., 2024).

Historically, there has been extensive research on anomalies within both the statistical and machine learning literature. For example, much work within statistics has been on developing inference methods that can cope with the effect of outlying observations.

Notable contributions include the foundational work of Huber (1981) and Rousseeuw and Yohai (1984) who introduced M- and S-estimators to provide robust estimates of the mean and variance, effectively mitigating the impact of outliers. In addition, several robust time series models have been developed to provide modelling approaches which are robust to anomalous observations. See, for example, work developing robust modelling for ARMA (Muler and Yohai, 2009), ARCH (Muler and Yohai, 2002), and GARCH (Muler and Yohai, 2008) settings. A recent overview was given by Maronna et al. (2019).

In parallel, the machine learning community has contributed a diverse set of methods tailored for anomaly detection in time series data. These methods range from tree-based approaches such as Isolation Forests (Liu et al., 2008) to deep learning-based techniques. For instance, Sakurada and Yairi (2014) apply auto-encoders to detect anomalies by reconstructing input data and measuring the reconstruction error. We refer the interested reader to reviews given by Chandola et al. (2009); Chalapathy and Chawla (2019); Zhao et al. (2023) and the references therein.

Traditionally, anomaly detection focuses on identifying individual data points or entire sequences that deviate from expected patterns. More recent research, however, has shifted toward integrating changepoint detection with anomaly detection to capture collective anomalies – groups of consecutive data points that appear anomalous only when considered together. A prominent example is the work of Fisch et al. (2020), who introduced the CAPA (Collective And Point Anomalies) algorithm, capable of detecting both collective and point anomalies effectively. It operates under the assumption (1) that the data sequence follows a constant baseline with a normally distributed noise, treating any deviations as anomalies:

[TABLE]

where $Y=(y_{1},\ldots,y_{n})$ is the observed time series, $A=(a_{1},\ldots,a_{n})$ represents point and collective anomalies, $B$ is the constant baseline, and $\varepsilon$ is the normally distributed noise. By employing a penalized cost method, CAPA efficiently identifies both the number and locations of point and collective anomalies, balancing model complexity with accuracy. However, real-world time series data often exhibit complex structures, including trends, seasonal patterns, and noise. Recognizing these patterns is vital for accurate analysis across various domains such as finance Liu et al. (2008); Kumar and Thenmozhi (2006), energy Hong and Fan (2016); Amjady and Keynia (2010), and healthcare Clifford et al. (2006); Liu et al. (2015). An example of such data is displayed in Figure 1, consisting of a seasonal component and a quadratic trend:

[TABLE]

Here the anomaly component $A$ represents unexpected deviations from normal behavior, the long-term deterministic trend $T$ captures the overall systematic direction of the data over time, seasonality $S$ models recurring fluctuations, and noise $\varepsilon$ accounts for random, unexplained variability. A full description of the model is given in Section 2.1.

Distinguishing between genuine anomalies and normal variations within such data presents a significant challenge, particularly when trends and seasonal effects are pronounced. In such situations, CAPA and other variants face limitations. The presence of strong trend and seasonal components can obscure or mask anomalies, leading to false positives, as seen in Figure 2.

Several decomposition methods have been introduced to separate time series data into trend, seasonal, and residual components. Examples include Seasonal and Trend decomposition using Loess (STL) (Cleveland et al., 1990), its version robust to outliers RobustSTL (Wen et al., 2019), and Prophet (Taylor and Letham, 2018). They are designed for a model without (collective) anomalies, i.e. model (2) without the anomaly component $A$ or $A$ limited to point anomalies, but allow for a more flexible trend. Decomposition may improve anomaly detection by isolating trend and seasonal components, allowing for a more focused analysis of residuals where anomalies may be more easily identified. Wen et al. (2019) demonstrate that such decomposition enhances the precision of anomaly detection by revealing hidden irregularities in the residuals. Despite the utility of decomposition methods, they struggle with time series that contain collective anomalies, which can be inadvertently incorporated into the estimated trend or seasonal components, complicating their identification and separation from regular patterns. This issue can misrepresent the anomalies and diminish the effectiveness of the decomposition approach, as shown in Figure 3. We remark that RobustSTL was too slow to be applied to this data example with $5\,000$ observations, but examples with less observations demonstrated that its decomposition has similar issues.

Consequently, current methods are not well-suited for handling such complex anomaly structures, highlighting the need for more robust approaches. To overcome the limitations of existing methods, we propose Seasonal Trend Anomaly Detection, (STAD), a novel approach extending existing anomaly detection methods by allowing a polynomial trend and seasonality as part of the typical pattern. At the same time, STAD decomposes the data into four components: anomaly, trend, seasonality, and residuals. By convention the baseline is included in the trend component. As discussed before none of the components can be estimated reliably by existing robust estimators when the other three components are present.

The method we propose is built upon a lag-differenced data, where the lag is a multiplicative of the seasonal length to remove the seasonal component. Standard robust methods can effectively estimate the polynomial trend at presence of point anomalies but struggle with collective anomalies, since they are non-symmetric and contain potentially a large number of anomalous points. We aim to overcome this issue by finding a subsample of the data that is free of collective anomalies and large enough to estimate the trend reliably. To this end, we propose a systematic sampling strategy that gives us a collection of subsamples. Under weak assumptions at least one of them will be free of collective anomalies and allow us to estimate the trend reliably. Once the trend has been removed, the seasonality can be estimated using robust estimators. Once the seasonality is estimated, the trend can be estimated without that differencing is required. Therefore, to improve accuracy we re-estimate the trend and seasonal component once before estimating the anomalies. Section 2 gives full details of the proposed methodology.

In Section 3 we use statistical theory to confirm the excellent detection properties of STAD. We show that it successfully decomposes the time series by providing asymptotic L2 bounds for the error of its estimation of the anomaly, trend and seasonal component. Furthermore, we show consistent estimation of the number of anomalies and that estimation of start and end points have a small asymptotic error. These results substantially improve on the results in (Fisch et al., 2020) for CAPA. Not only that we allow for the presence of a polynomial trend and seasonality with minor assumptions on them, but we also overcame the limitation that anomalies are at most of order $\sqrt{n}$ while we mirror the condition on the size of the anomalies despite allowing for a much more complex setting.

The finite sample performance is studied via simulation studies in Section 4. STAD estimates anomalies almost as good as its oracle version that knows the trend and seasonal component unless anomalies are extremely long (more than a quarter of the total time series) or the seasonal component has very few repetitions (less than five). In Section 5 we use STAD to detect anomalies in the price of energy data. A brief illustration of the effectiveness of our method is given in Figures 4 and 5, in which we have applied STAD to the data shown Figure 1. Different to CAPA and STL, confer Figures 2 and 3, it does not include false positives and decomposes the data successfully into anomaly, trend, seasonal, and residual components.

2 Methodology

We first detail our model, see Section 2.1, and then STAD. In Section 2.2 we explain how we estimate the anomaly component once trend and seasonality have been removed. In Section 2.3 we propose a robust method to estimate the trend at presence of collective anomalies, but under the assumption that the seasonality has been removed, either by estimation or as detailed in Section 2.4 by differencing the data. Similarly, Section 2.5 details the estimation of the seasonal component at presence of collective anomalies, but when the trend has been removed. Finally, STAD, a combination of the previous steps, is proposed in Section 2.6. All estimators are summarized in Algorithms 1, 2, 3, and 4 in Appendix A.

2.1 Model

As motivated in the introduction, we assume that the observation vector $Y=(y_{1},\ldots,y_{n})$ is given by the composition

[TABLE]

with anomalies $A=(a_{1},\ldots,a_{n})$, trend $T=(t_{1},\ldots,t_{n})$, seasonality $S=(s_{1},\ldots,s_{n})$, and noise $\varepsilon=(\varepsilon_{1},\ldots,\varepsilon_{n})$. The four components are defined as follows:

The anomalies split further into collective and point anomalies

[TABLE]

The collective anomalies are defined as

[TABLE]

with $K\in\mathbb{N}$ being the number, $0\leq b_{1}<e_{1}\leq\cdots\leq b_{K}<e_{K}\leq n$ being the start and end points, and $\mu_{1},\ldots,\mu_{K}$ being the magnitudes of the collective anomalies. Furthermore, let $\mathcal{A}=\bigcup_{k=1}^{K}(b_{k},e_{k}]$ be the union of all anomalous intervals. Let $O=\{o_{1},\ldots,o_{m}\}$ be the index set of the point anomalies, i.e. $a_{i}^{\operatorname{point}}\neq 0$ if and only if $i\in O$.

The trend is given by a polynomial of degree $Q$, i.e.

[TABLE]

with regression matrix $X=(x_{i,j})_{i=1,\ldots,n;\,j=1,\ldots,Q+1}$ having entries

[TABLE]

and coefficients $\beta=(\beta_{0},\ldots,\beta_{Q})\in\mathbb{R}^{Q+1}$.

We assume that the seasonality has a known period $P$, and is given by

[TABLE]

with $\vartheta_{0},\ldots,\vartheta_{P-1}\in\mathbb{R}$. For identifiability, we further assume

[TABLE]

i.e. the global constant is assumed to be specified by the trend, without a baseline in the anomaly component or a contribution by the seasonality.

Finally, we assume the noise $\varepsilon_{1},\ldots,\varepsilon_{n}$ to be i.i.d. centered Gaussian distributed with variance $\sigma_{0}^{2}$, i.e.

[TABLE]

The global standard deviation $\sigma_{0}$ can be pre-estimated using robust difference based estimators (Frick et al., 2014; Dette et al., 1998). More precisely, we use the IQR-based estimator from Pein et al. (2017):

[TABLE]

where $\Delta y=(y_{i+1}-y_{i})_{i=1}^{n-1}$ are the first order differences, $\text{IQR}(\Delta y)$ is the interquartile range of $\Delta y$, and $\text{IQR}_{\Phi}$ is the interquartile range of a standard normal distribution. The factor $1/\sqrt{2}$ corrects for the variance inflation introduced by differencing.

2.2 Anomaly Detection

We now present our methodology for estimating the collective and point anomalies. To this end, we assume that trend and seasonality have been removed, up to small estimation errors $R,W\in\mathbb{R}^{n}$, i.e. our method is based on the model

[TABLE]

To infer the number and locations of anomalous segments, we use a penalized cost approach analogous to CAPA (Fisch et al., 2022). We employ a change-in-mean version and assume that the trend estimation has removed any baseline, i.e. we apply CAPA with baseline $\mu_{0}=0$. More precisely, for any vector of observations $Z\coloneq(z_{1},z_{2},\ldots,z_{n})$ we define the costs for a constant segment from $i$ to $j$ as

[TABLE]

Let $\tilde{\mathcal{A}}\coloneq\bigcup_{k=1}^{\tilde{K}}(\tilde{b}_{k},\tilde{e}_{k}]$ be a candidate set of anomalies. Here, $\tilde{b}_{k}$ and $\tilde{e}_{k}$ denote the start and end of the $k$-th anomalous segment, respectively. Let further $\tilde{\mu}=(\tilde{\mu}_{1},\ldots,\tilde{\mu}_{\tilde{K}})$ be the candidate vector for the anomaly magnitudes. We then define their penalised costs as

[TABLE]

By convention, we set $\tilde{e}_{0}=0$ and $\tilde{b}_{\tilde{K}+1}=n$ and for any constant $a\in\mathbb{R}$ we write $Z-a$ for the vector $(z_{1}-a,\ldots,z_{n}-a)$. The first two terms in (12) quantify the deviation of observations from the baseline [math] and the anomaly magnitude respectively. The final term, $\tilde{K}\lambda_{\operatorname{coll}}\hat{\sigma}_{0}^{2}\log(n)$ is a SIC-type penalty. The cost function can readily be extended to allow for point anomalies. Specifically, let $\tilde{\mathcal{O}}\subset\{1,\ldots,n\}$ be a candidate set of point anomalies. Then, (12) is replaced by

[TABLE]

The reformulation of the first term in (12) follows directly from (11). The second term is unchanged and for the final term note that the unpenalised cost term is equal to [math] for a point anomaly, since the estimated function value is equal to the observation value. Similar to the collective anomaly penalty in (12), we add a SIC-type penalty $\lambda_{\operatorname{point}}\hat{\sigma}_{0}^{2}\log(n)$. By default, we recommend $\lambda_{\operatorname{coll}}=4$ for collective anomalies and $\lambda_{\operatorname{point}}=3$ for point anomalies. The reasoning behind this choice is discussed in Section 4.4.

The final estimate of the anomalies is obtained by minimizing the penalized cost

[TABLE]

and the final estimate of the anomaly component is $\hat{A}=(\hat{a}_{1},\ldots,\hat{a}_{n})$ with

[TABLE]

The penalized cost is additive over segments and can therefore be minimized via dynamic programming, leveraging the shared zero mean across non-anomalous segments (Fisch et al., 2022). To calculate (14), we employ the Generalized Functional Pruning Optimal Partitioning (GFPOP) software package (Hocking et al., 2022).

2.3 Trend Estimation

We now present our approach for polynomial trend estimation in the presence of anomalies, in particular collective anomalies. In other words, we assume that the seasonality has been removed up to a small estimation error $W\in\mathbb{R}^{n}$, that is, we model our data by

[TABLE]

Our procedure is based on finding an index set $U\subset\{1,\ldots,n\}$, such that $U\cap\mathcal{A}=\emptyset$, and $U$ is large enough and well structured to permit a reliable estimation of the trend. Specifically, we propose in Section 2.3.2 a sampling procedure that obtains a collection of index sets $\mathcal{U}$. For each $U\in\mathcal{U}$, the trend is then estimated by a robust estimator using only the observations in $U$, see Section 2.3.1. To select our final estimate we apply our anomaly detection method to the residuals of the trend estimate and choose the one that minimises the costs (12), full details are given in Section 2.3.3.

2.3.1 M-estimator

Given an index set $U\subset\{1,\ldots,n\}$, which may include both anomalous and non-anomalous observations, we estimate $\beta$ by an M-estimator using solely the data in $U$. We aim to obtain a good trend estimate if $U\cap\mathcal{A}=\emptyset$. We use the following wording in the following.

Definition 1*.*

We call an index set $U\subset\{1,\ldots,n\}$ anomalous if and only if at least one of its points is part of a collective anomaly, i.e. $U\cap\mathcal{A}\neq\emptyset$.

If $U$ is free of any anomalies, then least squares regression is best; however, since point anomalies may be present, a robust trend estimator is preferable. Following extensive experimentation, we adopt the M-estimator (Huber, 1964) with Tukey’s biweight loss. For any observation vector $Z=(z_{1},\ldots,z_{n})$ and any index set $U$, we define an M-estimator for the polynomial trend $T$ in (16) by

[TABLE]

where $\rho(\cdot)$ is the robust loss function, $x_{i,q}$ are the polynomial regressors (5) and $\hat{\sigma}_{0}$ is the estimated scale parameter (9). For the loss function $\rho$, we choose Tukey’s biweight loss function (Beaton and Tukey, 1974),

[TABLE]

Unlike for instance Huber’s loss function which only reduces the contribution of outliers to the loss, the biweight function does not increase for outliers beyond the threshold $c$; set to 4.685 by default. M-estimator (Huber, 1964) with Tukey’s biweight loss is robust against substantial outlier contamination while maintaining statistical efficiency (Özlem Gürünlü Alma, 2011).

The non-convex optimization challenge coming from using Tukey’s biweight loss is alleviated by our sampling approach as we rerun the algorithm multiple times for different observations. We employ the widely used iteratively reweighted least squares (IRLS) algorithm (Beaton and Tukey, 1974) to compute (17), details are provided in Appendix C.

2.3.2 Sampling Procedure

We aim at sampling a collection of index sets $\mathcal{U}$ that contains at least one index set $U$ that is non-anomalous, i.e. $U\cap\mathcal{A}=\emptyset$, and secondly that is large enough and well structured to allow a good estimate of $\beta$.

We use the following structured partitioning approach to sample $U$. First, we divide the index set $\{1,\ldots,n\}$ into $2(Q+1)$ disjoint segments $u_{1},\ldots,u_{2(Q+1)}$ of roughly equal length (be made more precise later). We will sample a sub-segment from either each odd or from each even indexed segment. More precisely, let $l\sim\operatorname{Bern}(\frac{1}{2})$ be a Bernoulli-distributed random variable. We then sample sub-segments from $u_{l},u_{l+2},\ldots,u_{l+2Q}$ (see below). This constructions ensures that the sampled sub-segments are well separated from each other, which is a key requirement to estimate the $Q+1$ parameters in $\beta$ well.

We now detail how we sample from $u_{k}$. We further partition each segment $u_{k}$ into $B$ sub-segments of roughly equal length. In total, we obtain

[TABLE]

sub-segments

[TABLE]

We then sample one sub-segment from each selected segment. More precisely, we obtain the index set

[TABLE]

where l\sim\operatorname{Bern}\big{(}\frac{1}{2}\big{)} is a Bernoulli-distributed random variable and m_{k}\sim\text{Uniform}\big{(}\{1,\ldots,B\}\big{)}. All random variables are independent of each other.

In summary, suppose appropriately chosen tuning parameters (see Section 2.3.3), the proposed sampling procedure obtains a collection of index sets $\mathcal{U}$ that ensures the following: With high probability at least one of the $U\in\mathcal{U}$ is non-anomalous, see Lemma 3. At the same time each $U\in\mathcal{U}$ is constructed such that if it is non-anomalous, then the trend is estimated well (with high probability), see Lemma 9. A key step was Lemma 4 in which we have shown that the determinant of $X_{U}^{\top}X_{U}$ is lower bounded. This lemma relies heavily on the fact that the $\tilde{u}_{(2k+l)B+m_{k}}$’s are well separated from each other.

2.3.3 Trend estimate selection

To obtain a good estimation of the trend $T$, we repeat the sampling procedure multiple times. We also have to select the tuning parameter $B$ and polynomial degree $Q$.

To this end, let $\mathbf{B}=\{B_{1},\ldots,B_{b}\}$ and $\mathbf{Q}=\{Q_{1},\ldots,Q_{q}\}$ be the candidate sets for $B$ and $Q$, respectively. For each combination $\tilde{B}\in\mathbf{B}$ and $\tilde{Q}\in\mathbf{Q}$ we repeat the sampling procedure in Section 2.3.2 $J$ times. We obtain the collection of index sets

[TABLE]

and trend estimates

[TABLE]

with $\hat{T}_{U}\coloneq X\hat{\beta}(U)$ and $\hat{\beta}(\cdot)$ as defined in (17).

We then select the final trend estimate $\hat{T}$ by selecting the one that minimizes the penalized cost of the anomaly estimator in (12) applied to $Y-\hat{T}_{U}$ plus a penalty for the degree of the trend, i.e.

[TABLE]

The tuning parameter $J$ is chosen to be sufficiently large such that with high probability at least one of the $J$ sampled segments is non-anomalous, confer Lemmas 2 and 3. By default we suggest $J=20$ based on the simulations in Section 4.3. The optimal choice of the tuning parameter $B$ is contingent on anomaly characteristics, e.g. number and length of collective anomalies. If $B$ is too small and collective anomalies are long, then all samples in $\mathcal{U}$ may intersect with collective anomalies. Conversely, if $B$ is too large, then the sample $U_{j,B,Q}$ is too small to allow a good estimation of the trend. These points are further illustrated in Appendix B. Finally, the correct polynomial degree $Q$ should be contained in $\mathbf{Q}$, though a too large degree may still allow good estimates. Based on empirical studies, in the absence of further information we recommend $\mathbf{B}=\{1,3,5\}$ and $\mathbf{Q}=\{0,1,2,3\}$ as they provide a good trade-off between statistical efficiency and computational demands.

2.4 Initial trend estimation

To estimate the trend without interference from the seasonality, we apply differencing to $Y$ with lag $D=C\cdot P$, where $C\in\mathbb{N}$. More precisely, we use the differenced series $Y_{D}$, with entries

[TABLE]

which has no seasonal component, confer Lemma 1. The differenced sequence has trend

[TABLE]

where the regression matrix $X_{D}$ has entries

[TABLE]

for $i=1,\ldots,n-D,\ q=1,\ldots,Q$, and parameter vector $\boldsymbol{\beta}_{D}=(\beta_{1},\ldots,\beta_{Q})$. Lemma 1 confirms this trend and details the effects of differencing on the noise and anomaly components, respectively. We then estimate $\boldsymbol{\beta}_{D}$ and hence all entries of $\boldsymbol{\beta}$ except $\beta_{0}$ by applying the trend estimation method from Section 2.3 to $Y_{D}$. We denote the resulting parameter estimate by

[TABLE]

and we obtain the initial trend $\hat{T}_{0}$ as

[TABLE]

The selection of an appropriate multiplier $C$ is crucial, as an inappropriate choice can lead to non-detectable anomalies. Theoretical justification is provided in Lemma 1. By default, $C$ is chosen to be an integer such that $D=C\cdot P\approx 0.1n$ and if $P\geq n/10$, then $C=1$.

2.5 Seasonal Component Estimation

We now introduce a robust estimator for the seasonal component that allows for the presence of anomalies. To this end, we assume that the trend has been removed, up to a small estimation error $R\in\mathbb{R}^{n}$. Under this assumption, the data follows the model:

[TABLE]

We also allow the estimator to take into account an estimation of the anomalies $\hat{\mathcal{A}}$, obtained by (23). Let $\hat{\mathcal{A}}=\emptyset$ for the initial seasonality estimation $\hat{S}_{0}$.

For $p=0,\ldots,P-1$, to estimate $\vartheta_{p}$ we use the index set

[TABLE]

and Tukey’s M-estimator

[TABLE]

with $\hat{\beta}$ being the M-estimator from (17) and $\beta\in\mathbb{R}$. For robustness, we initialize the estimator with the median of $Y_{V_{p}}$ rather than using Ordinary Least Squares (OLS). The final seasonal estimate $\hat{S}=(\hat{s}_{1},\ldots,\hat{s}_{n})$ is then defined as

[TABLE]

This ensures that $\hat{S}$ satisfies (7).

2.5.1 Smoothing

If the underlying seasonality is assumed to be smooth, we further improve the estimate $\hat{S}$ by applying spline smoothing. This enhances accuracy and interpretability. We choose spline smoothing due to its statistical robustness, computational efficiency, and flexibility in capturing complex seasonal patterns (De Boor, 1978), but other smoothing approaches would be suitable as well. More precisely, we apply cubic smoothing splines to $(\hat{\vartheta}_{0},\ldots,\hat{\vartheta}_{P-1})$ with the smoothing parameter selected by Generalized Cross-Validation. To this end, we use the R function smooth.spline.

2.6 STAD

We are now ready to define our main method STAD. It is the following combination of the previous estimators.

Step 0: Scale parameter estimation

We estimate $\sigma_{0}$ by $\hat{\sigma}_{0}$ as defined in (9).

Step 1: Initial Estimation

Obtain an initial trend estimation $\hat{T}_{0}$ by applying the proposed robust trend estimator in (23) to a differenced sequence of $Y$, which has removed the seasonality. Subsequently, obtain an initial estimate of the seasonality $\hat{S}_{0}$ by applying the robust seasonality estimator in (27) to the residuals $Y-\hat{T}_{0}$.

Step 2: Trend Estimation

Obtain a refined trend estimation $\hat{T}$ by applying the robust trend estimator in (23) to the residuals $Y-\hat{S}_{0}$. This process simultaneously estimates anomalous intervals, denoted by $\hat{\mathcal{A}}_{0}$.

Step 3: Seasonality Estimation

Obtain a refined seasonality estimation $\hat{S}$ by applying the robust seasonality estimator in (27) to the residuals $Y-\hat{T}$, but excluding the estimated anomalous intervals $\hat{\mathcal{A}}_{0}$.

Step 4: Anomaly Detection

Finally, we estimate the anomalous intervals $\mathcal{A}$ and the anomaly component $A$ by our modified version of the CAPA algorithm, detailed in Section 2.2, applied to the residuals $Y-\hat{S}-\hat{T}$.

3 Theory

With our methodology now established, we discuss theoretical properties of the STAD approach in detecting anomalies and decomposing the time series into anomalies, trend, seasonality, plus a remainder. Specifically in Theorem 1 we establish asymptotic L2-bounds on the error of each component’s estimate individually. Theorem 1 also shows that STAD estimates the number of anomalies consistently and their start and end points are close to the true ones. More precisely, we show that $\hat{\mathcal{A}}\in\mathcal{B}$, where

[TABLE]

Here, the errors $\kappa_{k}$ are defined as

[TABLE]

with $C_{1}>0$ being a large enough constant such that the last step in the proof of Lemma 15 holds. The maximum with $4^{Q+1}$ ensures that we have locally enough observations. For Theorem 1 to hold we require the following assumptions.

Assumption 1*.*

There exists a $\tilde{\delta}>0$ such that

[TABLE]

Furthermore, we require that

[TABLE]

and $D=CP$ is such that

[TABLE]

Equation (30) in Assumption 1 is a standard assumption on the length and size of an collective anomaly to be detectable. It mirrors Assumption 2 in (Fisch et al., 2020) for the CAPA algorithm in the absence of trend and seasonality. Note that the result $\hat{\mathcal{A}}\in\mathcal{B}$ is identical, up to the definition of $\kappa_{k}$ as we consider a change-in-mean instead of a change-in-mean-and-variance problem. The second requirement (31) ensures that anomalies are separated and hence individually detectable after differencing and is not required in the absence of either trend or seasonality.

Assumption 2*.*

The total length of collective anomalies is bounded. More precisely, there exists a $\gamma<1$ such that

[TABLE]

Assumption 2 ensures that we have locally enough non-anomalous observations to estimate trend and seasonality well. Remarkably, the total length of anomalies can be of order $n$. This shows that we have overcome the substantial limitation of CAPA given by Assumption 2(b) in (Fisch et al., 2020) that anomalies can be at most of order $\sqrt{n}$ to allow a good estimate of the baseline. The $4(Q+1)$-term is required for a good trend estimation, while the $8$ is required for a good seasonality estimation.

Assumption 3*.*

1. (i)

There exists a large enough $B\in\mathbf{B}$ such that (39), (48), and (50) are satisfied. 2. (ii)

We have that $Q\in\mathbf{Q}$.

Assumption 3 is on the tuning parameters $Q$ and $B$ to ensure that we can estimate the trend well. It requires the true order of the polynomial $Q$ to be in the set $\mathbf{Q}$. This could be weakened to only a $\tilde{Q}\geq Q$ being required. Secondly, the tuning parameter $B$ has to be large enough. We remark that the precise conditions work well for the asymptotic theory, but are not of practical relevance. We refer to Section 2.3.3 and Appendix B for a discussion of its choice.

Assumption 4*.*

There exits a constant $C_{2}>0$ such that

[TABLE]

and

[TABLE]

Assumption 4 is only used to bound the variance estimate $\hat{\sigma}_{0}^{2}$, see Proposition 1. The assumption could be relaxed by allowing a small proportion of points to fail conditions (33) and (34).

Under these assumptions, Theorem 1 shows that we have successfully solved the structural challenge to decompose a time series into anomalies, trend, and seasonality and to estimate collective anomalies when the typical distribution is not stationary. Our theory even accounts for the initial variance estimate. However, as, to the best of our knowledge, finite sample bounds for robust estimators are not available, we assume that we have no point anomalies and use the least squares estimator instead of Tukey’s estimator, i.e. the loss function $\rho$ is given by the $L_{2}$-norm. To estimate the seasonality we stop after the initial median, i.e. (27) is replaced by the median of the observations indexed by $V_{p}$. Using robust estimators instead will have slightly worse statistical efficiency, but better robustness to potential anomalies in $U$, either small parts of collective anomalies or point anomalies.

Theorem 1**.**

We assume the model from Section 2.1. Let $\hat{A},\hat{T},\hat{S}$, and $\hat{\mathcal{A}}$ be the STAD estimators as defined in Section 2.6. Suppose Assumptions 1, 2, 3, and 4 hold. Let $J\to\infty$, as $n\to\infty$. Then,

[TABLE]

Furthermore, there exists a constant $C_{3}>0$ such that

[TABLE]

as $n\to\infty$.

Proof.

See Appendix D. ∎

4 Simulation study

We now assess STAD’s finite sample performance by simulation studies. We compare STAD’s performance against oracle versions of itself that know the trend (Oracle-Trend), the seasonality (Oracle-Seasonality), and both components (Full Oracle). For all methods we use the default choices for tuning parameters as detailed in Section 2. A comparison with other methods is omitted, as most existing anomaly methods focus on point anomalies only and to the best of our knowledge no existing methods achieves acceptable results for detecting collective anomalies at the presence of significant trends and seasonality such that a systematic comparison is not meaningful, c.f. the illustrations in Section 1.

Our simulation studies are based on synthetic time series of length $n=5000$, generated from the model in Section 2.1. We chose the trend component to be $t_{i}=2(i/n)^{2}-2(i/n)$, a sine-wave seasonality $s_{i}=2\sin(2\pi i/P)$, where $P$ is the period length, and the noise component is $\epsilon_{i}\sim\mathcal{N}(0,1)$. Simulations (not displayed) showed that results are stable in the choice of the concrete form of the trend and seasonality, so for brevity we do not vary their form in the following. Instead we consider various settings for the anomalies and also vary in one study the period length $P$, which is otherwise chosen as $P=250$.

We found that the anomaly length, the number of seasonal repetitions, and of course the signal to noise ratio have the most dominant effect on the performance of all contrasted methods. We demonstrate the influence of the first two points in the most simplest setting with a single anomaly. In Sections 4.1 and 4.2 we vary the length of the single anomaly and the period length, respectively. Finally, in Section 4.5 we consider a more evolved setting with multiple anomalies which considers different signal to noise ratios but also different anomaly lengths. In between we discuss tuning parameter choices. Most tuning parameters have been theoretically discussed in Section 2. So, in Sections 4.3 and 4.4 we focus on the number of repetitions $J$ and the penalty, respectively.

All simulations are repeated $1\;000$ times. We then evaluate the methods based on their average accuracy in detecting the collective anomalies. We define an estimation as positive if the estimated number of anomalies is equal to the true number of anomalies and if start and end points differ by at most $\lfloor\log(n)\rfloor$. We omit other metrics for brevity and since we found for instance that the average mean square error for the anomaly, trend, and seasonal component estimation shows qualitatively the same result as the accuracy. Also changing the $\lfloor\log(n)\rfloor$ tolerance level for the locations of anomalies allows for the same comparisons.

4.1 Impact of anomaly length

We consider a single collective anomaly and vary its length. The starting position is drawn uniformly from all valid positions. To create a challenging setting, we also adjust the magnitude such that Full Oracle has slightly more than 90% accuracy.

Table 1 confirms that STAD performs very well and detects the anomaly almost as good as the oracles. Only when the anomaly is longer than a quarter of the total length of the time series its performance degenerates significantly. In these setting it fails to estimate the trend well, as evidenced by the performance of Oracle Season which performs more or the same as STAD. This is because both methods struggle to find a non-anomalous sample $U$ as even for large $B$ too many sub-segments are anomalous.

4.2 Effect of the seasonal length

In this simulation study we vary he seasonal period $P$. To this end, we reconsider a representative case from Table 1 with an anomaly length of 768, but results do not depend much on this choice as long as the anomaly is not too long.

Table 2 confirms STAD’s excellent performance, matching the oracles, across a wide range of seasonal repeats. However, when there are less than five repeats of the season, c.f. $P=1500$, both STAD and Oracle-Trend show significant degradation in performance, while oracles that know the seasonality still perform well. In such settings a very small number of points is available to estimate each value in the seasonal component and hence the anomaly effects the estimation too much.

4.3 Choice of J

We now investigate choice of the tuning parameter $J$, the number of sampling iterations used to estimate the trend. We use an anomaly length of 1378 as this was a challenging scenario in our previous experiments, c.f. Table 1. We vary $J\in\{5,10,20,30,40,50\}$.

Table 3 shows that unsurprisingly the accuracy of STAD improves with increasing $J$. Since the accuracy improves only marginally after $20$, we recommend $J=20$ as a default choice. Choosing a larger value may improve results but at the expense of larger computational demands.

4.4 Effect of the penalty

We now investigate the sensitivity of STAD to the collective anomaly penalty $\lambda_{\operatorname{coll}}$. We contrast its performance in three scenarios. Firstly, an easy scenario with a short anomaly of length 278 and moderate jump size 1.21 with many repeats of the seasonality ($P=250$). Secondly, a challenging scenario with a long anomaly of length 1378 and magnitude $1.19$ and a low-frequency seasonality ($P=1000$). Thirdly, a setting with a high signal-to-noise ratio, using the same short anomaly and seasonality as in the first case but with a significantly larger jump size of $4.2$.

Figure 6 shows how the penalty parameter $\lambda_{\operatorname{coll}}$ affects detection accuracy. As expected, in all scenarios, accuracy first rises sharply, then enters a wide, stable plateau of near-optimal performance, and finally declines as the penalty becomes too large. Based on these results, by default we recommend $\lambda_{\operatorname{coll}}=4$ for collective anomalies and $\lambda_{\operatorname{point}}=3$ for point anomalies. Note that this choice is slightly larger than the usual $\operatorname{BIC}$ choice that would be $3$ and $2$, respectively. This is to account for the additional errors resulting from the trend and seasonality estimation. However, if the anomaly is of high magnitude then the penalty can be chosen much larger. We will make use of this finding in our application in Section 5.

4.5 Multiple collective anomalies

Finally, we consider a setting with multiple anomalies. Our setting is similar to the one in (Fisch et al., 2022), but we restrict ourself to changes-in-mean only. The number of anomalies is drawn randomly from a Poisson distribution with intensity $2.5$. The length $L$ of each anomaly is drawn from a Poisson distribution with mean $\mu_{L}\in\{30,100,400,800\}$ and the magnitude of each anomaly is sampled from a uniform distribution taking values between $a$ and $b$. We consider $(a,b)=(1,3)$ (weak signals) and $(a,b)=(2,5)$ (strong signals). All random variables are independent of each other.

Again, Table 4 shows that STAD performs very well, almost matching the performance of the oracles unless anomalies are very long, i.e. $\mu_{L}=800$. This matches the findings in Section 4.1. Unsurprisingly, the performance of all methods improves when the signal to noise ratio is higher, i.e. $(a,b)=(2,5)$.

5 Application

The utility of STAD on real data is illustrated by applying it to historical electricity export prices from the Octopus Agile export tariff for the East Midlands (EM) distribution network region in the United Kingdom111Energy data: https://agile.octopushome.net/historical-data. The time series represents the electricity export price, $Y_{t}$, in pence per kilowatt-hour (p/kWh). We exported the period from April 17, 2025 (00:00 UTC) to May 26, 2025 (23:30 UTC), but we excluded weekends to maintain a daily seasonal pattern and hence to have enough repetitions of the seasonality. The recording frequency is 30 minutes, yielding to a total of $N=1344$ data points. We chose this data because of its real-world dynamics, specifically its daily seasonality and a long-term trend.

We set $P=48$, $C=1$, $J=50$, $Q\in\{1,\dots,10\}$, $B\in\{1,\dots,8\}$, $\lambda_{\operatorname{coll}}=32$, and $\lambda_{\operatorname{point}}=24$. Multiple of these values are increased in comparison to our default values as the underlying trend and seasonal patterns are rather complicated and the real data is also unlikely to perfectly fit our model. Furthermore, anomalies of interest seem to have a large magnitude, so the simulations in Section 4.4 suggests that a larger penalty is possible. For comparison we apply CAPA to the dataset with $\lambda_{\operatorname{coll}}=200$, and $\lambda_{\operatorname{point}}=150$.

Figure 7 demonstrates that CAPA, even with an extremely large penalty, detects anomalies because of the underlying seasonality. In comparison STAD identifies three significant collective anomalies which appear to the naked eye as abnormal behavior. To validate our findings, we compared the price data to renewable energy generation from the UK National Grid222UK Generation Data: National Grid ESO Data Portal, https://data.nationalgrideso.com/.

The obtained anomaly, trend, and seasonality in Figure 8 appears to be reasonable. We also see that at least two anomalies matches times with extremely large energy generation using wind.

Appendix A Pseudocodes

Appendix B Illustration of the tuning parameter $B$

In Scenario 1 we consider a simple quadratic trend and a long anomaly. If we choose $B=1$, then regardless of how large $J$ is or whether $Q$ is correctly specified, the resulting estimate will be poorly as inevitably the second sub-segment overlaps with the anomaly. If we choose $B=3$ instead, then some of the samples miss the anomaly and provide a good estimate.

In Scenario 2 we consider a more complicated cubic trend and a short anomaly. $B=3$ shows once again good results. If we choose $B$ very large, e.g. $B=15$, then even when the sample $U$ is non-anomalous it does not lead to a good estimate, since the available number of observations in $U$ is too small and the estimate is extremely noisy.

Appendix C The Iteratively Reweighted Least Squares

For completeness, we detail the Iteratively Reweighted Least Squares (IRLS) algorithm that we use to compute the $M$-estimator.

Start at $t=0$ with an initial estimate $\hat{\beta}^{(0)}$. We use ordinary least squares for the tend estimator and the median for the seasonal estimate. 2. 2.

Compute the residuals

[TABLE]

and the corresponding weights

[TABLE]

where $\psi$ denotes the influence function, i.e. the derivative of the loss function $\rho$. 3. 3.

Update the estimate $\hat{\beta}$ using weighted least squares regression

[TABLE]

where $W^{(t)}=\text{diag}(w_{1}^{(t)},w_{2}^{(t)},\dots,w_{n}^{(t)})$. 4. 4.

Increase $t$ by $1$ and repeat steps 2 and 3 until the estimated coefficients converge, i.e.

[TABLE]

or a maximum of $500$ iterations is reached.

Appendix D Proofs

We will show in Section D.2 that $C_{4}\sigma_{0}\leq\hat{\sigma}_{0}\leq C_{5}\sigma_{0}$ for constants $0<C_{4}<1$ and $C_{5}>1$ with probability converging to one. In Section D.3 we will study the effects of taking differences of lag $D$. In Sections D.4, D.5, and D.6 we will show that we estimate the trend, the seasonality, and the anomalies well given achievable assumptions. Finally, Theorem 1 will be proven in Section D.7.

D.1 Additional notations

Throughout the proof, we make use of the following notations. For any observation vector $Z=(z_{1},\ldots,z_{n})$ and index set $U\coloneq\{i_{1},\ldots,i_{m}\}\subset\{1,\ldots,n\}$ we use the definitions $Z_{U}\coloneq(z_{i_{1}},\ldots,z_{i_{m}})$ and $\overline{Z}_{U}\coloneq\frac{1}{m}\sum_{iinU}z_{i}$. For the regression matrix $X$ we denote by $X_{U}$ the matrix restricted to the rows in $U$, i.e. $X_{U}=(x_{i,j})_{i\in U,j=1,\ldots,Q+1}$. For indices $1\leq i\leq j\leq n$ we also use $i:j$ for the index set $\{i,\ldots,j\}$, e.g. we write $Z_{i:j}$ instead of $Z_{\{i,\ldots,j\}}$. Furthermore, we denote by $\hat{\theta}_{i:j}$ the best polynomial fit of degree $Q$ to $Y_{i:j}$, i.e.

[TABLE]

D.2 Bound on the variance estimate

The following proposition bounds the variance estimate.

Proposition 1**.**

We assume the model from Section 2.1. Let $\hat{\sigma}_{0}^{2}$ be the variance estimator as defined in (9). Suppose Assumption 4 holds. Then, there exists constants $0<C_{4}\leq 1$ and $C_{5}\geq 1$ such that

[TABLE]

Proof.

Let $\Delta y$ be the first order difference of the observations as defined before (9), let $\operatorname{IQR}$ be the interquartile range, let $\operatorname{IQR}_{\Phi}$ be the interquartile range of a standard Gaussian distribution, let $q_{0.75}$ and $q_{0.25}$ be the 75% and 25% quantile, respectively. Then, for a constant $\tilde{C}_{1}>0$ to be chosen later, the definition of the interquartile range and a union bound yield

[TABLE]

In the following we focus on bounding the first probability, but the bounding the second probability follows similar steps. It follows from the definition of a quantile for the first inequality, from Assumption 4 for the second inequality, and from the fact that there are at most $2K$ different $i$’s for which $a_{i}-a_{i-1}\neq 0$, since this can only occur at the beginning or end of an anomaly, for the third inequality that

[TABLE]

We choose $\tilde{C}_{1}$ such that $\frac{\tilde{C}_{1}-2C_{2}}{\sqrt{2}}$ is the $80\%$ quantile of a standard Gaussian distribution. Then, the left hand side in the probability above is binomial distributed with size $n-1$ and success probability

[TABLE]

We use the following bound from (Feller and Morse, 1958, (3.5)),

[TABLE]

if $r\geq np$, where $S_{n}$ is Binomial distributed with size $n$ and success probability $p$. It follows that

[TABLE]

as $n\to\infty$. This completes the proof with $C_{5}=2\tilde{C}_{1}\sqrt{2}\operatorname{IQR}_{\Phi}$, as the second probability can be bounded following similar arguments. Also the proof of $C_{4}\sigma_{0}\leq\hat{\sigma}_{0}$ follows similar arguments. ∎

D.3 Differenced sequence

We detail in Lemma 1 the properties of the differenced sequence $Y_{D}$ and show that it satisfies the following assumptions, which are required for a good estimation of the trend.

Assumption 5*.*

There exists a $\tilde{\delta}>0$ such that

[TABLE]

and

[TABLE]

Assumption 6*.*

The total length of collective anomalies is bounded. More precisely, there exists a $\gamma<1$ such that

[TABLE]

Lemma 1**.**

We assume the model from Section 2.1. Suppose Assumptions 1, 2, and 3 hold. Then,

[TABLE]

where

[TABLE]

with $\beta_{D}=(\beta_{1},\ldots,\beta_{Q})$ and $(X_{D})_{i,q}=\left(\frac{i+D}{n}\right)^{q}-\left(\frac{i}{n}\right)^{q}$, is a polynomial of degree $Q-1$. Furthermore, $A_{D}=(a_{D,1},\ldots,a_{D,n-D})$ contains $2K$ anomalies that satisfy Assumptions 3, 5, and 6. Finally,

[TABLE]

is a multivariate Gaussian distributed with expectation $\mathbb{E}\left[\varepsilon_{D,i}\right]=0$ and covariance

[TABLE]

Proof of Lemma 1.

We obtain

[TABLE]

In the following we look at each component separately. We start with the seasonal component. Since $D=cP$, we have that

[TABLE]

This shows that $Y_{D}$ has no seasonal component.

Next, we look at the trend component. We obtain

[TABLE]

It remains to show that $\beta_{0}$ disappears and that it is a polynomial of degree $Q-1$. The difference of powers identity yields

[TABLE]

We conclude that $T_{D}$ is a polynomial of degree $Q-1$, with regression matrix $X_{D}$, and parameter vector $\beta_{D}=(\beta_{1},\ldots,\beta_{Q})$.

Next we consider the anomaly component. We focus for now on the contribution of the $k$-th component. The final anomaly term will be the sum over all components. Assumption 1 yields $D>e_{k}-b_{k}$ and hence $e_{k}-D<b_{k}$. We distinguish a few cases.

[TABLE]

We obtain two anomalies of magnitude $\theta_{k}$ and of length $e_{k}-b_{k}$. Hence, Assumption 5 is satisfied if all anomalies are well separated. It follows from Assumption 1 that

[TABLE]

and

[TABLE]

This shows that the anomaly component contains $2K$ anomalies that satisfy Assumption 5. Assumption 6 follows from Assumption 2, since the total length of anomalies is $2\sum_{k=1}^{K}e_{k}-b_{k}$. Assumption 3 was directly formulated for the differenced sequence.

Finally, we consider the noise component. We have that $\varepsilon_{D}$ is multivariate Gaussian distributed as the difference of two multivariate Gaussian distributions. We further have that

[TABLE]

and

[TABLE]

∎

D.4 Trend estimation

In this section we show that we can estimate the trend well. We first state the setting and then the proposition.

Setting 1*.*

We assume that $Y=A+T+W+\varepsilon$, with trend $T=X\beta$ and deterministic remainder term $W=(w_{1},\ldots,w_{{m_{n}}})$, where either ${m_{n}}=n$ or ${m_{n}}=n-D$. The remainder $W=(w_{1},\ldots,w_{{m_{n}}})$ is a periodic vector with period length $P$ for which there exist a constant $C_{6}>0$ such that

[TABLE]

Let $C_{4}\sigma_{0}\leq\hat{\sigma}_{0}\leq C_{5}\sigma_{0}$ be deterministic. Suppose that Assumptions 3, 5, and 6 hold. We further assume that there are no point anomalies, i.e. $O=\emptyset$. Finally, we assume that $Y,A,T,\varepsilon,X,\beta$ are either

1. (a)

given as defined in Section 2.1, 2. (b)

or replaced by the differenced versions $Y_{D},T_{D},A_{D},\varepsilon_{D},X_{D},\beta_{D}$ as defined in Lemma 1.

We have that the trend is a polynomial of degree either ${r_{Q}}=Q$ or ${r_{Q}}=Q-1$.

Proposition 2**.**

Given Setting 1(a) there exists a constant $C_{7}>0$ such that

[TABLE]

as $n\to\infty$.

It will take multiple lemmas to prove this proposition. In Section D.4.1 we show that with large probability at least one of samples $U_{B,Q,j}$ is non-anomalous. Combined with the results in Section D.4.2 this ensures that at least one of our trend estimates is good. In Section D.4.3 we bound noise terms and then in Section D.4.4 we show that a good trend estimate leads to lower CAPA cost then a bad trend estimate. This allows us to show that the selected trend estimate has lower error. As a side product we also show that $\hat{\mathcal{A}}\in\mathcal{B}$.

D.4.1 Anomalous free subsamples

To estimate the trend well, we require that at least one of the samples $U_{B,Q,j}$ is non-anomalous, c.f. Definition 1. Given Setting 1 we show in Lemma 2 that at most $\frac{V}{2(Q+1)}$ of the sub-segments $\tilde{u}_{1},\ldots,\tilde{u}_{V}$ are anomalous. Building on this, we show in Lemma 3 that there exists a non-anomalous sample $U\in\mathcal{U}$ (21) with probability at least $1-\left(\frac{1}{2}\right)^{J}$. We remark that the sub-segments $\tilde{u}_{1},\ldots,\tilde{u}_{V}$ and the sampling procedure in Section 2.3.2 in general explicitly depends on $B$ and $Q$. For brevity we have omitted it in the definitions of the segments and sub-segments.

Lemma 2**.**

Given Setting 1, if $n$ is large enough, then there exists $B\in\mathbf{B}$ and $Q\in\mathbf{Q}$ such that at most $\frac{V}{2(Q+1)}$ of the sub-segments $\tilde{u}_{1},\ldots,\tilde{u}_{V}$ (20) are anomalous.

Proof.

Assumption 3 ensures that the correct degree $Q\in\mathbf{Q}$ and that there exists a $B\in\mathbf{B}$ such that

[TABLE]

for a constant $0<\gamma<1$. Each sub-segment $\tilde{u}_{k}$ contains at least $\lfloor n/V\rfloor$ many points. An anomalous sub-segment is either fully part of a collective anomaly or at the boundary of at least one collective anomaly. At most

[TABLE]

many sub-segments are fully part of a collective anomaly. Furthermore, at most $2K$ many sub-segments are at the boundary of a collective anomaly. So, in total we have at most

[TABLE]

anomalous sub-segments. We will bound this number in the following. Let $n$ be large enough such that $n\geq\frac{V}{1-\gamma}$, and hence $n-V\geq\gamma n$. Then,

[TABLE]

Finally, (39) and (40) yield that the total number of anomalous sub-segments is upper bounded by

[TABLE]

∎

Lemma 3**.**

Given Setting 1, if $n$ is large enough, then there exists $B\in\mathbf{B}$ and $Q\in\mathbf{Q}$ such that

[TABLE]

Consequently,

[TABLE]

Proof.

Assumption 3 ensures that the correct degree $Q\in\mathbf{Q}$ and that there exists a $B\in\mathbf{B}$ such that (39) is satisfied. Let $N$ denote the total number of anomalous sub-segments. Let $N_{k}$ be the number of anomalous sub-segments within segment $u_{k}$. Then, Lemma 2 yields

[TABLE]

By definition, $\frac{V}{2(Q+1)}=B$ is also the number of sub-segments in each segment. Thus, the segment $u_{k}$ contains $B\geq N$ sub-segments of which $N_{k}$ are anomalous. We have that l\sim\operatorname{Bern}\big{(}\frac{1}{2}\big{)} is a Bernoulli-distributed random variable and m_{k}\sim\text{Uniform}\big{(}\{1,\ldots,B\}\big{)}. All random variables are independent of each other. We will also use that $\prod_{i=1}^{m}(1-a_{i})\geq 1-\sum_{i=1}^{m}a_{i}$ for any $a_{1},\ldots,a_{m}\in[0,1]$. We obtain

[TABLE]

This proves the first statement. For the second statement note that the samples $U\in\mathcal{U}$ are independent of each other. Therefore,

[TABLE]

This completes the proof. ∎

D.4.2 Trend estimation using a single sample

In this subsection we show that we can estimate the trend well if we have a sample $U=U_{B,Q,j}\subset\{1,\ldots,n\}$ from the sampling procedure in Section 2.3.2 for which $U\cap\mathcal{A}=\emptyset$ holds. We estimate the trend by

[TABLE]

We then show that the error $\hat{T}_{U}-T$ is small, i.e. there exists a constant $C_{8}>0$ such that

[TABLE]

holds with probability going to one. To this end, we lower bound in Lemma 4 the determinant of $X_{U}^{\top}X_{U}$. We use this to bound in Lemma 8 the trace of a matrix related to the covariance matrix, before we finally show (42) in Lemma 9.

Lemma 4**.**

Let $U\coloneq U_{B,Q,j}\subset\{1,\ldots,n\}$ be a sample from the sampling procedure in Section 2.3.2. Suppose that $X$ is given by Setting 1. Then there exists a constant $C_{9}>0$ such that

[TABLE]

We require multiple definitions and technical lemmas to show Lemma 4. Note that adding a multiple of one row to another row does not change the determinant. We apply Gauss reduction to $X_{U}^{\top}X_{U}$ to obtain an upper right diagonal matrix. The matrix after $s\in\{1,\ldots,{r_{Q}}\}$ steps is defined recursively by

[TABLE]

and

[TABLE]

The next lemma simplifies this iterative process.

Lemma 5**.**

Under the assumptions of Lemma 4 we have that

[TABLE]

where

[TABLE]

for any $u\in U$, and

[TABLE]

for any $(u_{0},u_{1})\in U^{2^{s-1}}\times U^{2^{s-1}}$.

Proof.

We use induction. We exemplify the induction for Setting 1(a), but Setting 1(b) follows the same steps. The induction beginning for $s=1$ is true by definition,

[TABLE]

Let the statement be true for $s$. Then,

[TABLE]

For the second last equality note that summands with $u_{0}=u_{1}$ are [math]. This proves the statement. ∎

By construction of $U$ (21) there exists $w_{1}<\cdots<w_{{r_{Q}}+1}\in\{1,\ldots,{m_{n}}\}$ such that $U=\bigcup_{j=1}^{{r_{Q}}+1}\{w_{j}+1,\ldots,w_{j}+\tilde{L}_{j}\}$ with $L-2\leq\tilde{L}_{j}\leq L$ for all $j$ and $L$, where there exists a constant $C_{10}>0$ such that $L=C_{10}n$. Let $W=\{w_{1},\ldots,w_{{r_{Q}}+1}\}$. The following lemma bounds the difference between using $U$ and $W$.

Lemma 6**.**

Under the assumptions of Lemma 4, for any $i,s\in\mathbb{N}$, $w\in W^{2^{s-1}}$, $l\in\{1,\ldots,L\}^{2^{s-1}}$ there exists a constant $C_{i,s}>0$ only depending on $s$ and $i$ such that

[TABLE]

Proof.

Let $X$ be given by Setting 1(a). We use induction to show the statement. For any $w\in W$ and $l\in\{1,\ldots,L\}$, the difference of powers identity yield

[TABLE]

Hence, the statements holds for $s=1$ with $C_{i,1}=i-1$.

Suppose the statement is true for $s$. Note that it follows from (44) and $m_{i}^{(1)}(u)=\left(\frac{u}{n}\right)^{i-1}\leq 1$ that here exists a constant $\tilde{C}_{1}>0$ such that \big{|}m_{s}^{(s)}(u)\big{|}\leq\tilde{C}_{1} for all $u\in U^{2^{s-1}}$. Thus, for any $w,\tilde{w}\in W^{2^{s-1}}$ and $l,\tilde{l}\in\{1,\ldots,L\}^{2^{s-1}}$,

[TABLE]

This completes the induction with C_{i,s+1}=2\big{(}\tilde{C}_{1}C_{i,s}+\tilde{C}_{1}C_{s,s}+C_{i,s}C_{s,s}\big{)}.

Let now $X$ be given by Setting 1(b). The proof follows the same steps with two modifications. The induction beginning follows from applying difference of powers identity twice. We obtain

[TABLE]

So, $C_{i,s}=(i-1)i$. The first three lines also show that

[TABLE]

Consequently, there exists a constant $\tilde{C}_{2}>0$ such that

[TABLE]

Hence, the statement follows from the same induction.

∎

We now lower bound $m_{i}^{(s)}$ for a special sequence of $w$’s. We consider $v=(v_{1},\ldots,v_{2^{s-1}})\in W^{2^{s-1}}$ of the form

[TABLE]

Lemma 7**.**

Under the assumptions of Lemma 4, we have that

[TABLE]

Proof.

Let $X$ be given by Setting 1(a). We show by induction that $m_{i}^{(s)}(v)=P_{s}(w)S_{i,s}(w)$, where

[TABLE]

The induction beginning is satisfied for $s=1$ by

[TABLE]

Suppose the statement is true for $s$. Let $w=(w_{1},\ldots,w_{s})$, $\tilde{w}=(w_{1},\ldots,w_{s-1},w_{s+1})$, and $\breve{w}=(w_{1},\ldots,w_{s+1})$. It follows from (44), the induction assumption, $S_{s,s}(w)=1$, and the fact that $v_{2^{s-1}+1:2^{s}}$ is equal to $v_{1:2^{s-1}}$ except that the last entry is $w_{s+1}$ instead of $w_{s}$ that

[TABLE]

By construction we have that

[TABLE]

Finally, the difference of power identity yields

[TABLE]

This completes the induction. Thus,

[TABLE]

Let $X$ now be given by Setting 1(b). We follow the same steps, but with

[TABLE]

and

[TABLE]

where $s$ rows appear. The induction beginning for $s=1$ follows from the difference of power identity

[TABLE]

For the induction step we focus on the last sum as all other terms follow the same pattern as before. Again, by difference of power identity,

[TABLE]

This completes the induction. Thus,

[TABLE]

∎

We are now able to proof the main Lemma.

Proof of Lemma 4.

Let $X$ be given by Setting 1(a). We show by induction that there exists constants $\tilde{C}_{Q,s}>0$ only depending on $Q$ and $s$ such that

[TABLE]

The induction begin is satisfied, since

[TABLE]

Suppose (47) holds for $s-1$. Then from Lemma 5 and the induction assumption it follows that

[TABLE]

It remains to show that there exists a constant $\tilde{C}_{1}>0$ such that

[TABLE]

In the following we first upper, then lower bound this expression.

It follows from (44) and $m_{i}^{(1)}(u)=\left(\frac{u}{n}\right)^{i-1}\leq 1$ that there exists a constant $\tilde{C}_{2}>0$ such that $m_{s}^{(s)}(u)\leq\tilde{C}_{2}$. Consequently,

[TABLE]

We now obtain a lower bound as well. Let $v$ be as defined in (46). Then, it follows from Lemmas 6 and 7 that

[TABLE]

We have that the constant $\tilde{C}_{1}$ exists if $B$ is chosen large enough such that

[TABLE]

This completes the induction to show (47). The proof of the statement follows, since $L\geq C_{10}n$ and hence

[TABLE]

with $C_{9}=C_{10}^{Q+1}\prod_{s=1}^{Q+1}\tilde{C}_{Q,s}$.

Let $X$ now be given by Setting 1(b). The proof follows the same steps. We show by induction that there exists constants $\tilde{C}_{Q,s}>0$ only depending on $Q$ and $s$ such that

[TABLE]

The induction begin is satisfied with $\tilde{C}_{Q,1}=Q$, since

[TABLE]

Suppose (49) holds for $s-1$. Then from Lemma 5 and the induction assumption it follows that

[TABLE]

It remains to show that there exists a constant $\tilde{C}_{3}>0$ such that

[TABLE]

In the following we first upper, then lower bound this expression.

We had shown in (45) that there exists a constant $\tilde{C}_{4}>0$ such that

[TABLE]

Consequently,

[TABLE]

We now obtain a lower bound as well. Let $v$ be as defined in (46). Then, it follows from Lemmas 6 and 7 that

[TABLE]

We have that the constant $\tilde{C}_{3}$ exists if $B$ is chosen large enough such that

[TABLE]

This completes the induction to show (49). The proof follows, since $L\geq C_{10}n$ and hence

[TABLE]

with $C_{9}=C_{10}^{Q}\prod_{s=1}^{Q}\tilde{C}_{Q,s}$. ∎

Lemma 8**.**

Suppose that the assumptions of Lemma 4 hold. Then there exists a constant $C_{11}>0$ such that

[TABLE]

Proof of Lemma 8.

For any matrix $A$ and any symmetric, positive definite matrix $B$, with dimensions such that the products exist, we have that $\operatorname{tr}(ABA^{\top})=\operatorname{tr}(BA^{\top}A)$. We obtain

[TABLE]

Let $X$ now be given by Setting 1(a). Then, we have that $0\leq X_{i,j}\leq 1$ for all $i=1,\ldots,n;\ j=1,\ldots,Q+1$. Consequently, $(X^{\top}X)_{i,j}\leq n$ for all $i=1,\ldots,Q+1;\ j=1,\ldots,Q+1$.

Next, we focus on the entries in the inverse matrix $(X_{U}^{\top}X_{U})^{-1}$. Since $X_{U}$ is a sub-matrix of $X$, we have once again that $(X_{U})_{i,j}\leq 1$ and $(X_{U}^{\top}X_{U})_{i,j}\leq n$. The inverse $(X_{U}^{\top}X_{U})^{-1}$ is given by the cofactor matrix of $X_{U}^{\top}X_{U}$ divided by the determinant of $X_{U}^{\top}X_{U}$. Each entry in the cofactor matrix is given by the determinant of a $Q\times Q$ sub-matrix of $X_{U}^{\top}X_{U}$. Thus, each entry is bounded by $C_{12}n^{Q}$, where $C_{12}>0$ is a constant only depending on $Q$. Together with Lemma 4, which says $\operatorname{det}\left(X_{U}^{\top}X_{U}\right)\geq C_{9}n^{Q+1}$, we obtain

[TABLE]

Consequently,

[TABLE]

for all $i,j\in\{1,\ldots,Q+1\}$. Thus,

[TABLE]

Note that $C_{11}>0$ only depends on $Q$. This completes the proof if $X$ is given by Setting 1(a).

Let $X$ now be given by Setting 1(b). Then, we have that $0\leq X_{i,j}\leq Q\frac{D}{n}$ for all $i=1,\ldots,n-D;\ j=1,\ldots,Q$. Consequently, there exists a constant $\tilde{C}_{1}>0$ such that $(X^{\top}X)_{i,j}\leq\tilde{C}_{1}\frac{D^{2}}{n}$ for all $i=1,\ldots,Q;\ j=1,\ldots,Q$.

Next, we focus on the entries in the inverse matrix $(X_{U}^{\top}X_{U})^{-1}$. Since $X_{U}$ is a sub-matrix of $X$, we have once again that $(X_{U})_{i,j}\leq Q\frac{D}{n}$ and $(X_{U}^{\top}X_{U})_{i,j}\leq\tilde{C}_{1}\frac{D^{2}}{n}$. The inverse $(X_{U}^{\top}X_{U})^{-1}$ is given by the cofactor matrix of $X_{U}^{\top}X_{U}$ divided by the determinant of $X_{U}^{\top}X_{U}$. Each entry in the cofactor matrix is given by the determinant of a $(Q-1)\times(Q-1)$ sub-matrix of $X_{U}^{\top}X_{U}$. Thus, each entry is bounded by $C_{12}\left(\frac{D^{2}}{n}\right)^{Q-1}$, where $C_{12}>0$ is a constant only depending on $Q$. Together with Lemma 4, which says $\operatorname{det}\left(X_{U}^{\top}X_{U}\right)\geq C_{9}\left(\frac{D^{2}}{n}\right)^{Q}$, we obtain

[TABLE]

Consequently,

[TABLE]

for all $i,j\in\{1,\ldots,Q\}$. Thus,

[TABLE]

Note that $C_{11}>0$ only depends on $Q$. This completes the proof. ∎

We now bound the error $\hat{T}_{U}-T$ of the trend estimate given by $T_{U}$ in (41).

Lemma 9**.**

Let $U\coloneq U_{B,Q,j}\subset\{1,\ldots,{m_{n}}\}$ be a sample from the sampling procedure in Section 2.3.2 and suppose that $U\cap\mathcal{A}=\emptyset$. Suppose that our data follows Setting 1. Then there exists a constant $C_{8}>0$ such that

[TABLE]

Proof of Lemma 9.

The least squares regression estimator is given by $\hat{T}_{U}=X(X_{U}^{\top}X_{U})^{-1}X_{U}^{\top}Y_{U}$.

Since $U\cap\mathcal{A}=\emptyset$, we have that $Y_{U}=X_{U}\beta+W+\varepsilon_{U}$. Hence, $\hat{T}_{U}-T=X(X_{U}^{\top}X_{U})^{-1}X_{U}^{\top}(W+\varepsilon_{U})$. Therefore, we have that

[TABLE]

For the first term, we obtain

[TABLE]

where the final inequality follows from (38) and Lemma 8.

We now focus on the second term. For any symmetric positive definite matrix $A$ and any random variable $X$ with $\mathbb{E}[X]=0$, it follows from Kendrick (2002, (B.1)) that \mathbb{E}\left[X^{\top}AX\right]=\operatorname{tr}\big{(}A\operatorname{Cov}[X]\big{)}. Therefore,

[TABLE]

where $\operatorname{colSum}\left(|\operatorname{Cov}\left[\varepsilon_{U}\right]|\right)$ is the column sum of the element-wise absolute values of the covariance matrix. The final inequality follows from Lemmas 1 and 8.

Next, since $\varepsilon_{U}^{\top}X_{U}(X_{U}^{\top}X_{U})^{-1}(X^{\top}X)(X_{U}^{\top}X_{U})^{-1}X_{U}^{\top}\varepsilon_{U}\geq 0$, Markov’s inequality gives us

[TABLE]

It follows that (51) holds with $C_{8}\coloneq 2C_{6}C_{11}+2$. ∎

D.4.3 Uniform Bounds on the Noise

In this section we bound relevant noise terms.

Lemma 10**.**

Suppose that our data follows Setting 1. Then, there exists a constant $C_{13}<\infty$ such that

[TABLE]

as $n\to\infty$.

Proof.

We first consider the Setting 1(a). Note that $H_{i:j}$ is a projection matrix on a $Q+1$ dimensional column space. Hence, $H_{i:j}$ is symmetric, idempotent, and has rank $\min\{Q+1,j-i+1\}$. Consequently, (Mathai and Provost, 1992, Theorem 5.1.1) implies that

[TABLE]

is chi-squared distributed with $\min\{Q+1,j-i+1\}$ degrees of freedom. Thus, (Laurent and Massart, 2000, Lemma 1) yields

[TABLE]

for any $t>0$. Let $t=C_{13}\sigma_{0}^{2}\log(n)$ for a constant $C_{13}>0$. If $n$ is large enough, then

[TABLE]

Thus, if $n$ is large enough, we obtain

[TABLE]

If $C_{13}>6$, combining this bound with a union bound gives us

[TABLE]

as $n\to\infty$.

We now consider the Setting 1(b). Lemma 1 implies that $\varepsilon_{D}=\varepsilon_{D+1:n}-\varepsilon_{1:{n-D}}$. Consequently, $(\varepsilon_{D})_{i:j}=\varepsilon_{i+D:j+D}-\varepsilon_{i:j}$. Therefore,

[TABLE]

Following the same steps as in Setting 1(a) completes the proof. ∎

Lemma 11 bounds the mixed term between noise and a polynomial fit.

Lemma 11**.**

Let \tilde{T}=\big{(}\tilde{t}_{1},\ldots,\tilde{t}_{{m_{n}}}\big{)} be a polynomial of degree ${r_{Q}}$, then for any $1\leq i\leq j\leq{m_{n}}$, we have that

[TABLE]

Proof.

Since $\tilde{T}$ is a polynomial of degree ${r_{Q}}$, there exists a $\tilde{\beta}\in\mathbb{R}^{{r_{Q}}+1}$ such that $\tilde{T}=X\tilde{\beta}$. Let $C^{2}\coloneq\tilde{\beta}^{\top}X_{i:j}^{\top}X_{i:j}\tilde{\beta}$ and $\breve{\beta}=\frac{1}{C}\tilde{\beta}$. Then,

[TABLE]

It remains to solve the optimisation problem on the right. We use Lagrangian multipliers. We have to maximise

[TABLE]

We obtain

[TABLE]

Hence, $\lambda=\sqrt{\varepsilon_{i:j}^{\top}H_{i:j}\varepsilon_{i:j}}$. Consequently,

[TABLE]

Similarly, we have that

[TABLE]

This completes the proof. ∎

D.4.4 Bounding the cost difference

To show that we select the right trend estimate using the cost of CAPA, we contrast in the following the costs of a good trend estimate with the one of a bad trend estimate. We require multiple lemmas to bound these costs.

We start by rewriting the anomaly model as a changepoint model. Let $\tau_{2k}\coloneq b_{k}$ and $\tau_{2k+1}\coloneq e_{k}$, for $k=0,\ldots,K$. Furthermore, we define the function values $\theta_{2k-1}\coloneq\mu_{k}$, for $k=1,\ldots,K$, and $\theta_{2k}\coloneq 0$. We then define the jump sizes $\Delta_{2k-1}\coloneq|\mu_{k}|$ and $\Delta_{2k}\coloneq|\mu_{k}|$ for $k=1,\ldots,K$. We also use the notation $K_{\tau}\coloneq 2K$.

In several of the following lemmas we make use of the following inequalities. Let $R=(r_{1},\ldots,r_{{m_{n}}})$ be a polynomial of degree ${r_{Q}}$ for which

[TABLE]

holds. Suppose that

[TABLE]

Lemma 12 bounds the reduction in costs by having an estimated polynomial trend instead of the true parameter.

Lemma 12**.**

Suppose that our data follows Setting 1. Suppose that (52) and (53) hold. Then, there exists a constant $C_{15}>0$ such that for any $k\in\{0,\ldots,{K_{\tau}}\}$ and any $i,j$ such that $\tau_{k}<i\leq j\leq\tau_{k+1}$,

[TABLE]

Proof.

Note that $y_{l}=\theta_{k}+t_{l}+w_{l}+\varepsilon_{l},\ l=i,\ldots,j$. We further remark that $\theta_{k}+T$ is a polynomial of degree ${r_{Q}}$, so $\hat{\theta}_{i:j}=\theta_{k}+T_{i:j}+H_{i:j}W_{i:j}+H_{i:j}\varepsilon_{i:j}$. For the last equality before the inequalities note that $H_{i:j}$ is a projection matrix, i.e. $H_{i:j}^{\top}=H_{i:j}$, $H_{i:j}^{2}=H_{i:j}$. We obtain the first inequality from Lemma 11, the Cauchy-Schwarz inequality, and $x^{\top}H_{i:j}x\leq x^{\top}x$ for any vector $x$ since the largest eigenvalue of $H_{i:j}$ is $1$ as $H_{i:j}$ is a projection matrix. The second inequality follows from $\sum_{l=i}^{j}r_{l}^{2}\leq\sum_{l=1}^{{m_{n}}}r_{l}^{2}$ and $\sum_{l=i}^{j}w_{l}^{2}\leq\sum_{l=1}^{{m_{n}}}w_{l}^{2}$. The final inequality is a direct consequence of (52), (53),  and (38). In total, we obtain

[TABLE]

where $C_{15}\coloneq 2\sqrt{C_{13}C_{14}}+\sqrt{C_{14}C_{6}}+C_{14}+C_{13}+2\sqrt{C_{13}C_{6}}+C_{6}$. ∎

Lemma 13 bounds the increase in costs from having an extra estimated changepoint.

Lemma 13**.**

Suppose that our data follows Setting 1. Suppose that (53) hold. If $n$ is large enough, then for any $k\in\{0,\ldots,K\}$ and any $a,b,c$ such that $\tau_{k}<a\leq b\leq c\leq\tau_{k+1}$, then there exits a constant $C_{16}>0$ such that

[TABLE]

Proof.

Similarly to the proof of Lemma 12, we have that

[TABLE]

where the final inequality holds for large enough $n$. ∎

Lemma 14 lower bounds the quadratic difference between a step function and any polynomial of degree $Q$.

Lemma 14**.**

For any $Q\in\mathbb{N}$ there exists a constant $C_{17}>0$ only depending on $Q$ such that for any polynomial $f$ of degree $Q$, any step function with changepoint $\tau$ and function values $\theta_{0}$ and $\theta_{1}$, and any distance $\kappa_{0}\geq\max\{\kappa_{1},4^{Q+1}\}$, we have that

[TABLE]

The statement also holds with the roles of $\kappa_{0}$ and $\kappa_{1}$ reversed.

Proof.

Without loss of generality we assume that $\theta_{0}=0$ and $\theta_{1}>0$, as subtracting a constant or multiplying everything inside the squares with $-1$ do not change the result. Furthermore, without loss of generality assume that $\kappa_{0}$ and $\kappa_{1}$ are as given in the lemma. Otherwise consider the observations in reverse order which is not changing the costs.

We now distinguish the cases $\kappa_{1}<4^{Q+1}$ and $\kappa_{1}\geq 4^{Q+1}$. If $\kappa_{1}<4^{Q+1}$, then we distinguish the cases whether

[TABLE]

holds or not. If it does not hold, then

[TABLE]

So the statement holds with $C_{17}=4^{-Q-1}$.

If $\kappa_{1}\geq 4^{Q+1}$, then we distinguish the cases whether

[TABLE]

holds or not. If it does not hold, then we continue with proof below but we consider the observations after $\tau$ instead of the ones before.

The remaining proof uses the following idea. We consider a large number of small boxes around the step function. If $f$ is fully outside of any of these boxes, then the statement holds with $C_{17}$ resembling the size of the box. If $f$ is inside all boxes for at least one time point, then by using repeatedly the mean value theorem we will show that the $(Q+1)$-th derivative of $f$ must be non-zero and hence $f$ cannot be a polynomial of degree $Q$.

We start by constructing the boxes. For any $\iota=(i_{1},\ldots,i_{m})\in\{0,1\}^{m}$, $m\in\mathbb{N}$, let $I_{\iota}\coloneq\left(L_{\iota},R_{\iota}\right]$ with

[TABLE]

For notational convenience we extend the notation to the empty set and define $I_{\emptyset}\coloneq(\tau-\kappa_{0},\tau]$. Hence, $I_{(i_{1},\ldots,i_{m})}$ is of length $\kappa_{0}4^{-m}$ and $I_{(i_{1},\ldots,i_{m-1})}$ is the union of $I_{(i_{1},\ldots,i_{m-1},0)}$, $I_{(i_{1},\ldots,i_{m-1},1)}$, and the interval in between, which is of length $2\kappa_{0}4^{-m}$.

Let

[TABLE]

If there exists an $\iota\in\{0,1\}^{Q}$ such that for all $x\in I_{\iota}$ we have that

[TABLE]

then

[TABLE]

with $C_{17}>0$ being a constant only depending on $Q$, since $\lfloor R_{\iota}\rfloor-\lceil L_{\iota}\rceil\geq\frac{\kappa_{0}}{4^{Q}}-1$ and $\kappa_{0}\geq 4^{Q+1}$ ensuring that $I_{\iota}$ contains at least one observation.

Otherwise, we show by induction that for every $m\in\{Q,\ldots,0\}$ and for every $\iota\in\{0,1\}^{m}$ there exists an $x_{\iota}\in I_{\iota}$ for which

[TABLE]

Furthermore, for every $m\in\{Q,\ldots,-1\}$ there exists a $\xi_{m}\in(x_{\mathbf{1}_{m+1}},\tau+1]$, where $\mathbf{1}_{m}=(1,\ldots,1)\in\{0,1\}^{m}$ and $\xi_{Q}=\tau+0.5$ or $\xi_{Q}=\tau+1$, depending on whether (54) or (56) holds, for which

[TABLE]

For the first statement, the induction start for $m=Q$ is given by construction of the boxes and that $f$ is inside all boxes for at least one point. Suppose that the statement is true for an $m$. Then, by the mean value theorem we have that there exists a point

[TABLE]

for which

[TABLE]

since by construction the distance between the two points is at least $2\kappa_{0}4^{-m}$. The lower bound follows from the same arguments but with a minus in front. Hence, the first statement holds for $m-1$.

For the second statement, the induction start is given by $f(\xi_{Q})\geq\frac{\theta_{1}-\theta_{0}}{2}$ since (54) or (56) holds. Suppose that the statement holds for an $m$. Then, by the mean value theorem we have that there exists a point

[TABLE]

for which

[TABLE]

since by construction the distance between the two points is at most $\kappa_{0}4^{-m}+1\leq 2\kappa_{0}4^{-m}$. This completes the induction step.

All in all, we have shown that there exists a $\xi_{-1}\in(\tau-\kappa_{0},\tau+1]$ for which, by the definition of $b$,

[TABLE]

Thus, the $(Q+1)$-th derivative is non-zero and hence $f$ cannot be a polynomial of degree $Q$. This is a contradiction. Hence $f$ must be in the boxes and the statement holds as shown before. ∎

Lemma 15 bounds the increase in costs from missing a changepoint.

Lemma 15**.**

Suppose that our data follows Setting 1. Let $\hat{T}=X\hat{\beta}$ be any polynomial of degree ${r_{Q}}$. Suppose that (52) and (53) hold. If $n$ is large enough, then there exists a constant $C_{18}>0$ such that for any $k\in\{1,\ldots,2K\}$,

[TABLE]

Proof.

Similarly to the proof of Lemma 12, we obtain for the second and third term that

[TABLE]

For the first term, we have that

[TABLE]

Note that $\theta_{k-1}+T-\hat{T}$ and $\theta_{k}+T-\hat{T}$ are polynomials of degree ${r_{Q}}$. Hence, Lemma 11, $\sqrt{a}+\sqrt{b}\leq\sqrt{2a+2b}$ for any $a,b>0$, and (53) yield

[TABLE]

Similarly, the Cauchy-Schwarz inequality and (38) yield

[TABLE]

Moreover, note that $\hat{T}-T$ is a polynomial of degree ${r_{Q}}$. We can apply Lemma 14, since $\kappa_{k}\geq 4^{Q+1}$. Thus, Lemma 14 yields

[TABLE]

Therefore, all in all,

[TABLE]

where the last three inequalities follow from the definitions of $\kappa_{k}$, $C_{1}$ and $n$ being large enough. ∎

We can now bound the costs of CAPA between a good and any other trend estimate.

Lemma 16**.**

Suppose that our data follows Setting 1. Let $\hat{T}=X\hat{\beta}$ be any polynomial of degree ${r_{Q}}$. Let $\tilde{T}$ be a polynomial of degree ${r_{Q}}$ such that (52) holds for $R=(r_{1},\ldots,r_{{m_{n}}})\coloneq\tilde{T}-T$. Suppose that (53) holds. If $n$ is large enough, then for any $\hat{\mathcal{A}}\not\in\mathcal{B}$, $\hat{\mu}$ we have that

[TABLE]

Proof.

There exists a constant $\tilde{C}_{1}>0$ such that

[TABLE]

since $\mathbf{Q}$ is a finite set and $\hat{\sigma}_{0}\leq C_{5}\sigma_{0}$. Let

[TABLE]

with $\hat{K}$ being the number of anomalies, $\hat{s}_{1},\ldots,\hat{s}_{\hat{K}}$ and $\hat{e}_{1},\ldots,\hat{e}_{\hat{K}}$ being the points before the anomalies start and the last point of the anomalies, respectively. We require $0\leq\hat{s}_{1}<\hat{e}_{1}\leq\hat{s}_{2}<\cdots<\hat{e}_{\hat{K}-1}\leq\hat{s}_{\hat{K}}<\hat{e}_{\hat{K}}\leq{m_{n}}$. Furthermore, we define $\hat{e}_{0}=0$ and $\hat{s}_{\hat{K}+1}={m_{n}}$.

Furthermore, we rewrite them as a changepoint model: We use the notation $\hat{K}_{\tau}\coloneq 2\hat{K}$. Furthermore, we define $\hat{\tau}_{2k}\coloneq\hat{e}_{k}$ and $\hat{\tau}_{2k+1}\coloneq\hat{s}_{k}$, for $k=0,\ldots,\hat{K}$. Moreover, we have function values $\hat{\theta}_{2k-1}\coloneq\hat{\mu}_{k}$, for $k=1,\ldots,\hat{K}$, and $\hat{\theta}_{2k}\coloneq 0$, for $k=0,\ldots,\hat{K}$.

Using these definitions, we have that

[TABLE]

Furthermore, using the definitions for $\mathcal{A}$, we have that

[TABLE]

We will split the costs at $\tau_{k}\pm\kappa_{k}$. Note that Assumption 5 yields that $\tau_{k}-\kappa_{k}>\tau_{k-1}$ and $\tau_{k}+\kappa_{k}<\tau_{k+1}$ for large enough $n$.

We further remark that $\hat{T}+\hat{\theta}_{k}$ is a polynomial of degree ${r_{Q}}$ and hence by definition for any $k$, for any $\hat{\tau}_{k}\leq\tau\leq\hat{\tau}_{k+1}$ we have that

[TABLE]

Furthermore, $Y-\tilde{T}=\varepsilon+A+R+W$. Hence, on the set $\{\tau_{k}+1,\tau_{k+1}\}$, we have that $Y-\tilde{T}-\theta_{k}=\varepsilon+R+W$. Hence, for any $k$, for any $\tau_{k}\leq\tau\leq\tau_{k+1}$ we have that

[TABLE]

For $k\in\{0,\ldots,K_{\tau}\}$, let $\tilde{K}_{k}$ be the number of estimated change-points between $\tau_{k}+\kappa_{k}$ and $\tau_{k+1}-\kappa_{k+1}$. That is, there exists an $l$ such that

[TABLE]

We write $\tilde{\tau}_{k,i}\coloneq\hat{\tau}_{l+i}$ for $i=1,\ldots,\tilde{K}_{k}$. Moreover, we write $\tilde{\tau}_{k,0}\coloneq\tau_{k}+\kappa_{k}$ and $\tilde{\tau}_{k,\tilde{K}_{k}+1}\coloneq\tau_{k+1}-\kappa_{k+1}$. Note that $k=0$ and $k=K_{\tau}$ are special boundary cases, for which we have $\kappa_{0}=\kappa_{2K+1}=0$.

Similarly, for $k\in\{1,\ldots,K_{\tau}\}$, let $\hat{K}_{k}$ be the number of estimated change-points between $\tau_{k}-\kappa_{k}$ and $\tau_{k}+\kappa_{k}$. That is, there exists an $l$ such that

[TABLE]

We write $\hat{\tau}_{k,i}\coloneq\hat{\tau}_{l+i}$ for $i=1,\ldots,\hat{K}_{k}$. Moreover, we write $\hat{\tau}_{k,0}\coloneq\tau_{k}-\kappa_{k}$ and $\hat{\tau}_{k,\hat{K}_{k}+1}\coloneq\tau_{k}-\kappa_{k}$.

Using these definitions, we have that

[TABLE]

where

[TABLE]

and

[TABLE]

Lemmas 17 and 18 bound $T_{1,k}$ and $T_{2,k}$, respectively.

If there exists a $k\in\{1,\ldots,K_{\tau}\}$ such that $\hat{\tau}_{1},\ldots,\hat{\tau}_{\hat{K}}\not\in\{\tau_{k}-\kappa_{k},\ldots,\tau_{k}+\kappa_{k}\}$, then $\hat{K}_{k}=0$. Note that if $\hat{K}<K$, then such a $k$ must exist. Consequently, Lemmas 17 and 18 and (57) imply that

[TABLE]

where the final inequality follows from $n$ being large enough.

Furthermore, if $\hat{K}>K$, then there exists a $k\in\{0,\ldots,K_{\tau}\}$ such that $\tilde{K}_{k}\geq 1$ or there exists a $k\in\{1,\ldots,K_{\tau}\}$ such that $\hat{K}_{k}>1$. Consequently, Lemmas 17 and 18 and (57) imply that

[TABLE]

where the final inequality follows from $n$ being large enough.

All in all, if $\hat{\mathcal{A}}\not\in\hat{\mathcal{B}}$, then C\big{(}Y-\hat{T},\hat{\mathcal{A}},\hat{\mu}\big{)}-C\big{(}Y-\tilde{T},\mathcal{A},\mu\big{)}-\max_{Q\in\mathbf{Q}}Q\,\hat{\sigma}_{0}\log(n)>0, as $\hat{\mathcal{A}}\not\in\hat{\mathcal{B}}$ implies that one of the above cases must occur. ∎

Lemma 17**.**

Let $T_{1,k}$ be as defined in (58). Under the assumptions of Lemma 16, for any $k=0,\ldots,K_{\tau}$, we have that

[TABLE]

Proof.

If $\tilde{K}_{k}=0$, then the statement follows directly from Lemma 12. If $\tilde{K}_{k}\geq 1$, then applying Lemma 13 iteratively and finally Lemma 12 yields

[TABLE]

where the last inequality follows from $n$ being large enough. ∎

Lemma 18**.**

Let $T_{2,k}$ be as defined in (59). Under the assumptions of Lemma 16, for any $k=1,\ldots,K_{\tau}$, we have that

[TABLE]

Proof.

If $\hat{K}_{k}=0$, then Lemma 15 gives the desired statement.

From now on, let $\hat{K}_{k}\geq 1$. We split the costs at $\tau_{k}$. Let $\hat{K}_{k,1}$ be the number of estimated change-points between $\hat{\tau}_{k,0}+1$ and $\tau_{k}$ and $\hat{K}_{k,2}$ be the number of estimated change-points between $\tau_{k}+1$ and $\hat{\tau}_{k,\hat{K}_{k}+1}$. Furthermore, let $\hat{\tau}_{k,1,l}\coloneq\hat{\tau}_{k,l}$, for $l=0,\ldots,\hat{K}_{k,1}$, and $\hat{\tau}_{k,2,l}\coloneq\hat{\tau}_{k,l+\hat{K}_{k,1}}$, for $l=1,\ldots,\hat{K}_{k,2}+1$. We also define $\hat{\tau}_{k,1,\hat{K}_{k,1}+1}\coloneq\tau_{k}\eqcolon\hat{\tau}_{k,2,0}$. Since, $\hat{K}_{k}=\hat{K}_{k,1}+\hat{K}_{k,2}$ and $\hat{K}_{k}\geq 1$, we have that $\hat{K}_{k,1}\geq 1$ or $\hat{K}_{k,2}\geq 1$ or both. Without loss of generality let $\hat{K}_{k,1}\geq 1$. Otherwise, we can switch the role of the two terms in the following. Thus, we obtain,

[TABLE]

with

[TABLE]

and

[TABLE]

If $\hat{K}_{k}=1$, then $\hat{K}_{k,1}=1$ and $\hat{K}_{k,2}=0$. Thus, Lemma 12 yields

[TABLE]

and

[TABLE]

Thus,

[TABLE]

If $\hat{K}_{k}>1$, then $\hat{K}_{k,1}>1$ and $\hat{K}_{k,2}>0$ (or both). Applying Lemma 13 iteratively and then the same steps as before gives us

[TABLE]

for large enough $n$. ∎

We are now ready to show that the estimated anomalies are within the set $\mathcal{B}$ with probability converging to one. This result is of interest by its own but also helps us to show that we select a good trend estimate.

Proposition 3**.**

Suppose that our data follows Setting 1. Let $J\to\infty$, as $n\to\infty$. Let $\hat{\mathcal{A}}$ be as defined in (23). Then,

[TABLE]

Proof.

Suppose that there exists an arbitrary but fixed polynomial $\tilde{T}\in\mathcal{T}$ (22) of degree ${r_{Q}}$ such that $(\tilde{T}-T)^{\top}(\tilde{T}-T)\leq\tilde{C}_{1}\sigma_{0}^{2}\log(n)$. Assume further that

[TABLE]

for a constant $\tilde{C}_{2}>0$. Let $n$ be large enough. Then, Lemma 16 implies that $\hat{\mathcal{A}}\in\mathcal{B}$, since $\mathcal{A}\in\mathcal{B}$ and $\hat{\mathcal{A}}$ is the minimiser of the costs in (23).

Consequently, it follows from a union bound, Lemmas 3, 9, and 10, the independence of $\mathcal{U}$ and $\varepsilon$, and $J\to\infty$, as $n\to\infty$, that there exists $\tilde{C}_{1},\tilde{C}_{2}>0$ such that

[TABLE]

as $n\to\infty$. ∎

Proof of Proposition 2.

Let $\tilde{T}=(\tilde{t}_{1},\ldots,\tilde{t}_{{m_{n}}})$ be as in Proposition 3, i.e. $(\tilde{T}-T)^{\top}(\tilde{T}-T)\leq\tilde{C}_{1}\sigma_{0}^{2}\log(n)$. Assume further that

[TABLE]

for a constant $\tilde{C}_{2}>0$. In the proof of Proposition 3 we have shown that both events hold simultaneously with probability converging to one as $n\to\infty$.

We now show, by contradiction, that the statement in Proposition 2 holds. If not, then there must exists $\hat{T}$, $\hat{\mathcal{A}}$, and $\hat{\mu}$ such that

[TABLE]

for a sequence $\gamma_{n}\to\infty$, as $n\to\infty$, and

[TABLE]

Furthermore, Proposition 3 yields that $\hat{\mathcal{A}}\in\mathcal{B}$. Note that $\mathcal{A}\in\mathcal{B}$ by definition.

Then Lemma 11 and the Cauchy-Schwarz inequality for the first inequality, as well as (60), (38), and Lemma 10 for the second inequality imply that

[TABLE]

where $\tilde{C}_{3}\coloneq 2\sqrt{C_{13}C_{8}}+2\sqrt{C_{6}C_{8}}+C_{8}$.

Let $\{\tilde{\tau}_{0},\tilde{\tau}_{1},\ldots,\tilde{\tau}_{\tilde{K}},\tilde{\tau}_{\tilde{K}+1}\}=\{e_{0},b_{1},\hat{b}_{1},e_{1},\hat{e}_{1},\ldots,b_{K},\hat{b}_{K},e_{K},\hat{e}_{K},b_{K+1}\}$ with $\tilde{\tau}_{0}\leq\tilde{\tau}_{1}\leq\cdots\leq\tilde{\tau}_{\tilde{K}}\leq\tilde{\tau}_{\tilde{K}+1}$ and $\tilde{K}=4K$. Then, the first inequality is required by (62) and the remaining steps follow the same arguments as before as well as (57)

[TABLE]

This can only hold if there exists a constant $\tilde{C}_{4}>0$ such that

[TABLE]

This, combined with the Cauchy-Schwarz inequality, and (61) gives us

[TABLE]

This can only hold if there exists a constant $\tilde{C}_{5}>0$ such that

[TABLE]

Since $\tilde{K}$ is finite, there must exist an $l\in\{0,\ldots,\tilde{K}\}$ and a constant $\tilde{C}_{6}>0$ such that

[TABLE]

We will show the statement by contradiction. We distinguish cases depending on which interval $[\tilde{\tau}_{l}+1,\tilde{\tau}_{l+1}]$ is.

Suppose $[\tilde{\tau}_{l}+1,\tilde{\tau}_{l+1}]$ is non-anomalous and to be estimated non-anomalous, i.e. there exists a $k\in\{1,\ldots,K+1\}$ such that $\tilde{\tau}_{l}=\max\{e_{k-1},\hat{e}_{k-1}\}$ and $\tilde{\tau}_{l+1}=\min\{b_{k},\hat{b}_{k}\}$. Then, $\hat{a}_{i}=a_{i}=0$ for all $\tilde{\tau}_{l}<i\leq\tilde{\tau}_{l+1}$. Thus,

[TABLE]

This is a contradiction to (64). 2. 2.

Suppose $[\tilde{\tau}_{l}+1,\tilde{\tau}_{l+1}]$ is either anomalous or to be estimated anomalous, but not both, i.e. there exists a $k\in\{1,\ldots,K\}$ such that either $\tilde{\tau}_{l}=\min\{b_{k},\hat{b}_{k}\}$ and $\tilde{\tau}_{l+1}>\min\{b_{k},\hat{b}_{k}\}$ or $\tilde{\tau}_{l}<\max\{e_{k},\hat{e}_{k}\}$ and $\tilde{\tau}_{l+1}=\max\{e_{k},\hat{e}_{k}\}$. We focus on the first case. The second case follows in the same way but with $\kappa_{0}$ and $\kappa_{1}$ reversed.

Let $\kappa_{1}=\tilde{\tau}_{l+1}-\tilde{\tau}_{l}$ and $\kappa_{0}=\tilde{\tau}_{l}-\tilde{\tau}_{l-1}$. Because of the definition of $\kappa_{k}$ and Assumption 5 we have that $\kappa_{0}>\max\{\kappa_{1},4^{Q+1}\}$ for $n$ large enough. We have that $a_{i}=\hat{a}_{i}=0$ for all $\tilde{\tau}_{l-1}<i\leq\tilde{\tau}_{l}$ since this interval must be non-anomalous and estimated to be non-anomalous by the order of the $\tau_{l}$’s. Secondly, there exists a $\hat{\mu}\in\mathbb{R}$ such that $a_{i}-\hat{a}_{i}=\hat{\mu}$ for $\tilde{\tau}_{l}+1\leq i\leq\tilde{\tau}_{l+1}$. Consequently, since $T-\hat{T}$ is a polynomial of degree $Q$, Lemma 14, and (64) imply

[TABLE]

This is a contradiction to (63). 3. 3.

Suppose $[\tilde{\tau}_{l}+1,\tilde{\tau}_{l+1}]$ is anomalous and to be estimated anomalous, i.e. there exists a $k\in\{1,\ldots,K\}$ such that $\tilde{\tau}_{l}=\max\{b_{k},\hat{b}_{k}\}$ and $\tilde{\tau}_{l+1}=\min\{e_{k},\hat{e}_{k}\}$.

Let $\kappa_{0}=\tilde{\tau}_{l+1}-\tilde{\tau}_{l}$ and $\kappa_{1}=\tilde{\tau}_{l+2}-\tilde{\tau}_{l+1}$. Because of the definition of $\kappa_{k}$ and Assumption 5 we have that $\kappa_{0}>\max\{\kappa_{1},4^{Q+1}\}$ for $n$ large enough. If $\kappa_{1}=0$, then consider the next interval and we may have to switch the roles of $\kappa_{0}$ and $\kappa_{1}$. If there exists a sequence $\tilde{\gamma}_{n}\to\infty$ such that $\sum_{i=\tilde{\tau}_{l+1}+1}^{\tilde{\tau}_{l+2}}(t_{i}-\hat{t}_{i})^{2}>\tilde{\gamma}_{n}\sigma_{0}^{2}\log(n)$, then it follows from the same arguments that lead to (64) that there exists a constant $\tilde{C}_{7}>0$ such that $\sum_{i=\tilde{\tau}_{l+1}+1}^{\tilde{\tau}_{l+2}}(a_{i}-\hat{a}_{i})^{2}>\tilde{C}_{7}\tilde{\gamma}_{n}\sigma_{0}^{2}\log(n)$. Consequently, considering $l-1$ and case 2. leads to a contradiction. Otherwise, there exists a constant $\tilde{C}_{8}$ such that

[TABLE]

Then, from the same arguments as before, there exists constants $\tilde{C}_{9},\tilde{C}_{10}>0$ such that

[TABLE]

This is a contradiction to (63).

In summary all cases lead to a contradiction. This completes the proof. ∎

D.5 Seasonality

We show that we can estimate the seasonality well if the following assumption and setting is assumed. Recall that we use $\tilde{\vartheta}_{p}\coloneq\operatorname{median}(Y_{V_{p}})$.

Assumption 7*.*

The total length of collective anomalies is bounded, i.e.

[TABLE]

Setting 2*.*

We assume that

[TABLE]

where $Y,A,S,\varepsilon$ are as defined in Section 2.1, $B=(b,\ldots,b)\in\mathbb{R}^{n}$ is a constant baseline with $b=0$ or $b=\beta_{0}$, and $R=(r_{1},\ldots,r_{n})$ being a deterministic vector such that there exists a constant $C_{19}>0$ such that

[TABLE]

We further assume that

[TABLE]

where $\hat{\mathcal{A}}$ is deterministic. Finally, suppose that Assumptions 5 and 7 holds.

Under Setting 2 we show the following proposition.

Proposition 4**.**

Suppose that $Y$ follows Setting 2. Let $\tilde{\vartheta}_{p}\coloneq\operatorname{median}(Y_{V_{p}})$. Then, there exists a constant $C_{20}>0$ such that

[TABLE]

as $n\to\infty$.

We require a few definitions and lemmas before we can prove Proposition 4. Let

[TABLE]

where the latter matches (26). By construction, $\tilde{V}_{0}=\{P,2P,\ldots,\tilde{n}_{0}P\}$, with $\tilde{n}_{0}\coloneq\lfloor\frac{n}{P}\rfloor$, and $\tilde{V}_{p}=\{p,p+P,\ldots,p+(\tilde{n}_{p}-1)P\}$, with $\tilde{n}_{p}=\tilde{n}_{0}$ or $\tilde{n}_{p}=\tilde{n}_{0}+1$.

Lemma 19**.**

Under the setting of Proposition 4, with probability converging to one, $\tilde{\vartheta}_{p}$, i.e. the median of the observations in $V_{p}$, is in the interval

[TABLE]

with $\gamma\coloneq\sqrt{\frac{\log(n)}{n}}\sigma_{0}$ and $\tilde{q}$ being the median of the observations

[TABLE]

We will prove this lemma after we have shown a few more technical lemmas. The following lemma bounds the cardinality of $V_{0}$ and how many observations are anomalous.

Lemma 20**.**

Under the setting of Proposition 4, if $n$ is large enough, then there exists a constant $C_{21}>0$ such that the set $V_{0}$ contains $\frac{5}{6}\frac{n}{P}-C_{21}\leq N_{0}\leq\frac{n}{P}$ many observations. Moreover, at most $\frac{1}{6}\frac{n}{P}+C_{21}$ many observations are anomalous and at least $\frac{5}{6}\frac{n}{P}-C_{21}$ many observations are non-anomalous.

Proof.

Because of (65), $V_{0}$ contains at least $\tilde{n}_{0}$ many observations minus the points that are part of an estimated anomaly. A season containing an estimated collective anomalies is either fully part of an estimated collective anomaly or at the boundary of at least one estimated collective anomaly. The definition of $\mathcal{B}$ yield that at most

[TABLE]

seasons are fully part of an estimated collective anomaly. Furthermore, at most $2K$ many seasons are at the boundary of an estimated anomaly. So, in total $V_{p}$ contains at least

[TABLE]

many observations. Let $C_{21}=2K+1$. Let $n$ be large enough such that

[TABLE]

This holds for large $n$ because of Assumption 5. Then, it follows from Assumption 7 that

[TABLE]

This shows the first statement. For the second statement note that a season containing anomalous points is either fully part of a collective anomaly or at the boundary of at least one collective anomaly. At most $\frac{\sum_{k=1}^{K}e_{k}-b_{k}}{P}$ seasons are fully part of a collective anomaly and at most $2K$ many seasons are at the boundary of an anomaly. So, in total at most

[TABLE]

many observations in $V_{0}$ are anomalous. It follows from Assumption 7 that

[TABLE]

The number of non-anomalous points follows conversely. ∎

We have two technical lemmas on how to bound the density and cumulative distribution function of a standard Gaussian random variable.

Lemma 21**.**

Let $\phi$ be the density of a standard Gaussian random variable. If $x\leq 1$, then

[TABLE]

Proof.

We obtain the Taylor series

[TABLE]

The statement follows from the fact that the terms have alternating signs and decay in absolute values, since $x\leq 1$. ∎

Lemma 22**.**

Let $\Phi$ be the cumulative distribution function of a standard Gaussian random variable. If $x\leq 1$, then

[TABLE]

Proof.

Taylor’s theorem gives us

[TABLE]

for an $\tilde{x}\in[0,x]$. The statement follows from $\frac{\tilde{x}^{3}}{6\sqrt{2\pi}}\leq\frac{x^{3}}{6\sqrt{2\pi}}\leq\frac{1}{6\sqrt{2\pi}}x$, since $\tilde{x}\leq x\leq 1$. ∎

The following lemma bounds $\tilde{q}$.

Lemma 23**.**

Under the setting of Proposition 4, if $n$ is large enough, then there exists a constant $0<C_{22}<0.752$ and a constant $0<C_{23}<1$ such that

[TABLE]

Proof.

Let $n_{0}$ be the size of $V_{0}$. As the median, $\tilde{q}>C_{22}\,\sigma_{0}$ requires

[TABLE]

We obtain from Lemma 20 and $n_{0}=\lceil\frac{n}{P}\rceil\leq\frac{n}{P}-1$ that

[TABLE]

The random variable on the left hand side in the probability is Binomial distributed with size $\frac{5}{6}\frac{n}{P}-C_{21}$ and success probability $\mathbb{P}(\varepsilon_{1}\geq C_{22}\,\sigma_{0})$. We use the following bound from (Feller and Morse, 1958, (3.5)),

[TABLE]

if $r\geq np$, where $S_{n}$ is Binomial distributed with size $n$ and success probability $p$. It follows from Lemma 22 that

[TABLE]

Thus,

[TABLE]

∎

The following lemma bounds the probability that the noise is between $\tilde{q}$ and $\tilde{q}+\gamma$.

Lemma 24**.**

Under the setting of Proposition 4, if $n$ is large enough, then there exists a constant $C_{24}>0$ such that

[TABLE]

Proof.

Let $q\coloneq 0.752\,\sigma_{0}$. Then, if $n$ is large enough, there exists a constant $0<\tilde{C}_{1}<1$ such that

[TABLE]

In the first inequality we lower bound the integral by minimal value times the length of the interval of integration. We use Lemma 21 for the second and (69) for the third inequality. We obtain

[TABLE]

It follows from Lemma 23 and a union bound that

[TABLE]

with $C_{24}\coloneq\frac{1}{\sqrt{2\pi}}\left(1-\frac{1}{2}\tilde{C}_{1}^{2}\right)-C_{23}>0$. ∎

With these lemmas we are finally able to prove Lemma 19.

Proof of Lemma 19.

For any $p\in\{0,\ldots,P-1\}$, let $n_{p}$ be the number of observations in $V_{p}$. As $\tilde{\vartheta}_{p}$ is the median of the observations in $V_{p}$,

[TABLE]

We focus on the first probability, but bounding the second probability follows similar arguments. The following relies on the fact that $V_{0}$ and $V_{p}$ are similar for most observations. For $p\neq 0$, we have that

[TABLE]

The last inequality holds, since there are at most $4K+1$ many $j^{\prime}s$ for which either $a_{p+jP}\neq a_{(j+1)P}$ or $(j+1)P\in V_{0}$, but $p+jP\in V_{p}$, since this can only happen at the beginning or end of a collective anomaly or an estimated collective anomaly or for the last observation in $V_{p}$. The last line has the same distribution as

[TABLE]

For similar reasons we have that $n_{p}\geq n_{0}-2K$. Hence, there exists a constant $\tilde{C}_{1}>0$ such that

[TABLE]

If $p=0$, we can obtain the same inequality with $\tilde{C}_{1}=0$. By the definition of $\tilde{q}$ (66),

[TABLE]

Thus,

[TABLE]

where the final inequality follows from Lemma 20.

The random variable on the left hand side in the probability is Binomial distributed with size $\frac{5}{6}\frac{n}{P}-C_{21}$ and success probability $\mathbb{P}(\tilde{q}\leq\varepsilon_{1}\leq\tilde{q}+\gamma)$. We use the following bound from (Feller and Morse, 1958, (3.6)),

[TABLE]

if $r\leq np$, where $S_{n}$ is Binomial distributed with size $n$ and success probability $p$. From this and Lemma 24 we obtain

[TABLE]

as $n\to\infty$. This completes the proof as showing that the second probability converges to [math] follows similar arguments. ∎

We now in the position to show Proposition 4.

Proof of Proposition 4.

It follows from Lemma 19, $\hat{s}_{l}=\tilde{s}_{l}-\frac{1}{P}\sum_{p=0}^{P-1}\tilde{s}_{p}$ that

[TABLE]

as $n\to\infty$, with $C_{20}\coloneq 1+C_{19}$. ∎

D.6 Anomaly estimation

We will show that under the following setting we can estimate the anomaly component well. In Proposition 5 we show that $\hat{\mathcal{A}}\in\mathcal{B}$ with probability converging to one. In Proposition 6 we will bound the L2 error of the anomaly estimate. For both proofs we will make great use of similar results from Section D.4.

Setting 3*.*

We assume that $Y=A+R+W+\varepsilon$, with $A,\varepsilon,X,\beta$ as defined in Section 2.1. We further assume that there are no point any anomalies, i.e. $O=\emptyset$. The remainder $R=(r_{1},\ldots,r_{n})$ is a deterministic polynomial of order $Q$ for which there exists a constant $C_{19}>0$ such that

[TABLE]

The remainder $W=(w_{1},\ldots,w_{n})$ is a deterministic periodic vector with period length $P$ for which there exist a constant $C_{6}>0$ such that

[TABLE]

Let $C_{4}\sigma_{0}\leq\hat{\sigma}_{0}\leq C_{5}\sigma_{0}$ be deterministic. Suppose that Assumptions 3, 5, and 6 hold.

Proposition 5**.**

Given Setting 3 we have that

[TABLE]

Proof.

Let $\tilde{T}$ be the estimate that was used to obtain the remainder $R$ and let $\mathcal{T}=\{\tilde{T}\}$. Then, the statement is a direct consequence of Proposition 3. ∎

Proposition 6**.**

Given Setting 3 there exists a constant $C_{25}>0$ such that

[TABLE]

as $n\to\infty$.

Proof.

It follows from the proof of Propositions 3 and 5 that the events $\hat{\mathcal{A}}\in\mathcal{B}$ and

[TABLE]

holds simultaneously with a probability converging to one as $n\to\infty$. In the remaining proof we show that this implies the desired statement which completes the proof.

Since $\mathcal{A}\in\mathcal{B}$, it follows from the definition of $\hat{\mathcal{A}}$ that

[TABLE]

Also note that $\hat{\mathcal{A}}\in\mathcal{B}$ implies $\hat{K}=K$.

Let $\{\tilde{\tau}_{0},\tilde{\tau}_{1},\ldots,\tilde{\tau}_{\tilde{K}},\tilde{\tau}_{\tilde{K}+1}\}=\{e_{0},b_{1},\hat{b}_{1},e_{1},\hat{e}_{1},\ldots,b_{K},\hat{b}_{K},e_{K},\hat{e}_{K},b_{K+1}\}$ with $\tilde{\tau}_{0}\leq\tilde{\tau}_{1}\leq\cdots\leq\tilde{\tau}_{\tilde{K}}\leq\tilde{\tau}_{\tilde{K}+1}$ and $\tilde{K}=4K$. Then, we obtain from (72) for the first inequality, Lemma 11 and the Cauchy-Schwarz inequality for the second inequality, as well as Lemma 10, (70) and (71) that

[TABLE]

This and $(a+b+c)^{2}\leq 3a^{2}+3b^{2}+3c^{2}$ for any $a,b,c\in\mathbb{R}$ yield

[TABLE]

with $C_{25}\coloneq 48KC_{13}+12C_{19}+12C_{6}$. This completes the proof. ∎

D.7 Proof of Theorem 1

Proof of Theorem 1.

Proposition 1 implies that there exists constants $0<C_{4}\leq 1$ and $C_{5}\geq 1$ such that

[TABLE]

Let $\tilde{E}_{1}$ denote the event in (73). It follows from Lemma 1 that if the event $\tilde{E}_{1}$ holds, then the differenced series $Y_{D}$ follows Setting 1(b).

Given Setting 1(b), Proposition 2 implies that there exists a constant $\tilde{C}_{1}>0$ such that

[TABLE]

Let $\tilde{E}_{2}$ be the event in (74). Since $\hat{T}_{0}-T$ is a polynomial of degree $Q$ the event $\tilde{E}_{2}$ implies that there exists a constant $\tilde{C}_{2}>0$ such that

[TABLE]

Given that they have the same parameter vector, it also implies that there exists a constant $\tilde{C}_{3}>0$ such that

[TABLE]

Therefore, given $\tilde{E}_{2}$ the residuals $Y-X\hat{\beta}$ satisfy Setting 2.

Given Setting 2, Proposition 4 implies that there exists a constant $\tilde{C}_{4}>0$ such that

[TABLE]

Let $\tilde{E}_{3}$ be the event in (75). Then, given $\tilde{E}_{3}$ the residuals $Y-\hat{S}_{0}$ satisfy Setting 1(a).

Given Setting 1(a), Proposition 2 implies that there exists a constant $C_{7}>0$ such that

[TABLE]

as $n\to\infty$. Moreover, Proposition 3 ensures that there exists an anomaly estimate $\tilde{\mathcal{A}}$ for which

[TABLE]

Let $\tilde{E}_{4}$ be union of the events in (76) and (77). Since $\hat{T}-T$ is a polynomial of degree $Q$ the event $\tilde{E}_{4}$ implies that there exists a constant $\tilde{C}_{5}>0$ such that

[TABLE]

Therefore, given $\tilde{E}_{4}$ the residuals $Y-\hat{T}$ and $\tilde{\mathcal{A}}$ satisfy Setting 2.

Given Setting 2, Proposition 4 implies that there exists a constant $C_{20}>0$ such that

[TABLE]

as $n\to\infty$. Let $\tilde{E}_{5}$ be the event in (78). Then, given $\tilde{E}_{4}\cup\tilde{E}_{5}$ the residuals $Y-\hat{T}-\hat{S}$ satisfy Setting 3.

Given Setting 3, Proposition 5 implies that (35) holds. Finally, given Setting 3 Proposition 6 yields that there exists a constant $C_{25}>0$ such that

[TABLE]

as $n\to\infty$.

Thus, both statements in Theorem 1 hold given the events $\tilde{E}_{1},\tilde{E}_{2},\tilde{E}_{3},\tilde{E}_{4}$, and $\tilde{E}_{5}$ with $C_{3}=\max\{C_{7},C_{20},C_{25}\}$. The proof follows from a union bound as this is a finite number of events and all of them hold with probability converging to one as $n\to\infty$. ∎