# Detecting linear trend changes and point anomalies in data sequences

Hyeyoung Maeng and Piotr Fryzlewicz
Department of Statistics, London School of Economics

###### Abstract

We propose TrendSegment, a methodology for detecting multiple change-points
corresponding to linear trend changes or point anomalies in one dimensional
data. A core ingredient of TrendSegment is a new Tail-Greedy Unbalanced
Wavelet transform: a conditionally orthonormal, bottom-up transformation of
the data through an adaptively constructed unbalanced wavelet basis, which
results in a sparse representation of the data. The bottom-up nature of this
multiscale decomposition enables the detection of point anomalies and linear
trend changes at once as the decomposition focuses on local features in its
early stages and on global features next. To reduce the computational
complexity, the proposed method merges multiple regions in a single pass over
the data. We show the consistency of the estimated number and locations of
change-points. The practicality of our approach is demonstrated through
simulations and two real data examples, involving Iceland temperature data and
sea ice extent of the Arctic and the Antarctic. Our methodology is implemented
in the R package trendsegmentR, available from CRAN.

Keywords: change-point detection; bottom-up algorithms; piecewise-linear
signal; anomaly detection; wavelets

##  1 Introduction

Multiple change-point detection is a problem of importance in many
applications; recent examples include automatic detection of change-points in
cloud data to maintain the performance and availability of an app or a website
(James et al., , 2016), climate change detection in tropical cyclone records
(Robbins et al., , 2011), detecting exoplanets from light curve data (Fisch et
al., , 2018), detecting changes in the DNA copy number (Olshen et al., , 2004;
Jeng et al., , 2012; Bardwell et al., , 2017), estimation of stationary
intervals in potentially cointegrated stock prices (Matteson et al., , 2013),
estimation of change-points in multi-subject fMRI data (Robinson et al., ,
2010) and detecting changes in vegetation trends (Jamali et al., , 2015).

This paper considers the change-point model

$X_{t}=f_{t}+\varepsilon_{t},\quad t=1,\ldots,T,$ (1)

where $f_{t}$ is a deterministic and piecewise-linear signal containing $N$
change-points, i.e. time indices at which the slope and/or the intercept in
$f_{t}$ undergoes changes. These changes occur at unknown locations
$\eta_{1},\eta_{2},\ldots,\eta_{N}$. The $\varepsilon_{t}$’s are iid random
errors following the normal distribution with mean zero and variance
$\sigma^{2}$. Both continuous and discontinuous changes in the linear trend
are permitted. A point anomaly can be viewed as a separate data segment
containing only one data point. Therefore, if $f_{\eta_{i}}$ is a point
anomaly, then the two consecutive change-points that define it, $\eta_{i-1}$
and $\eta_{i}$, are linked via $\eta_{i-1}=\eta_{i}-1$ under the definition of
a change-point specified later in (15). Our main interest is in the estimation
of $N$ and $\eta_{1},\eta_{2},\ldots,\eta_{N}$ under some assumptions that
quantify the difficulty of detecting each $\eta_{i}$; therefore, our aim is to
segment the data into sections of linearity and/or point anomalies in $f_{t}$.
In particular, a point anomaly can only be detected when it has a large enough
jump size with respect to the signal levels to its right and left, while a
change-point capturing a small size of linear trend change requires a longer
distance from its adjacent change-points to be detected. Detecting both linear
trend changes and point anomalies is an important applied problem in a variety
of fields, including climate change, as illustrated in Section 5.

The change-point detection procedure proposed in this paper is referred to as
TrendSegment; it is designed to work well in detecting not only long trend
segments and point anomalies, but also short trend segments that are not
necessarily classified as point anomalies. The engine underlying TrendSegment
is a new Tail-Greedy Unbalanced Wavelet (TGUW) transform: a conditionally
orthonormal, bottom-up transformation for univariate data sequences through an
adaptively constructed unbalanced wavelet basis, which results in a sparse
representation of the data. In this article, we show that TrendSegment offers
good performance in estimating the number and locations of change-points
across a wide range of signals containing constant and/or linear segments
and/or point anomalies. TrendSegment is also shown to be statistically
consistent and computationally efficient.

In earlier related work regarding linear trend changes, Bai and Perron,
(1998) consider the estimation of linear models with multiple structural
changes by least-squares and present Wald-type tests for the null hypothesis
of no change. Kim et al.,  (2009) and Tibshirani et al.,  (2014) consider
‘trend filtering’ with the $L_{1}$ penalty and Maidstone et al.,  (2017)
detect changes in the slope with an $L_{0}$ regularisation via a dynamic
programming algorithm. Spiriti et al.,  (2013) study two algorithms for
optimising the knot locations in least-squares and penalised splines.
Baranowski et al.,  (2016) propose a multiple change-point detection device
termed Narrowest-Over-Threshold (NOT), which focuses on the narrowest segment
among those whose contrast exceeds a pre-specified threshold. Anastasiou and
Fryzlewicz,  (2018) propose the Isolate-Detect (ID) approach which
continuously searches expanding data segments for changes.

Keogh et al.,  (2004) mention that sliding windows, top-down and bottom-up
approaches are three principal categories which most time series segmentation
algorithms can be grouped into. Keogh et al.,  (2004) apply those three
approaches to the detection of changes in linear trends (but not point
anomalies) in 10 different signals and discover that the performance of
bottom-up methods is better than that of top-down methods and sliding windows,
notably when the underlying signal has jumps, sharp cusps or large
fluctuations. Bottom-up procedures have rarely been used in change-point
detection. Matteson and James,  (2014) use an agglomerative algorithm for
hierarchical clustering in the context of change-point analysis. Keogh et al.,
(2004) merge adjacent segments of the data according to a criterion involving
the minimum residual sum of squares (RSS) from a linear fit, until the RSS
falls under a certain threshold; but the lack of precise recipes for the
choice of this threshold parameter causes the performance of this method to be
somewhat unstable, as we report in Section 4.

As illustrated later in this paper, many existing change-point detection
methods for the piecewise-linear model fail in signals that include frequent
change-points or abrupt local features. The TGUW transform, which underlies
TrendSegment, is able to handle scenarios involving possibly frequent change-
points. It constructs, in a bottom-up way, an adaptive wavelet basis by
consecutively merging neighbouring segments of the data starting from the
finest level (throughout the paper, we refer to a wavelet basis as adaptive if
it is constructed in a data-driven way). This enables it to identify local
features at an early stage, before it proceeds to focus on more global
features corresponding to longer data segments.

Fryzlewicz,  (2018) introduces the Tail-Greedy Unbalanced Haar (TGUH)
transform, a bottom-up, agglomerative, data-adaptive transformation of
univariate sequences that facilitates change-point detection in the piecewise-
constant sequence model. The current paper extends this idea to adaptive
wavelets other than adaptive Haar, which enables change-point detection in the
piecewise-linear model (and, in principle, to higher-order piecewise
polynomials, but we do not pursue this in the current work). We emphasise that
this extension from TGUH to TGUW is both conceptually and technically non-
trivial, due to the fact that it is not a priori clear how to construct a
suitable wavelet basis in TGUW for wavelets other than adaptive Haar; this is
due to the non-uniqueness of the local orthonormal matrix transformation for
performing each merge in TGUW, which does not occur in TGUH. We solve this
issue by imposing certain guiding principles in the way the merges are
performed, which enables the detection of changes in the linear trend and
point anomalies. The TGUW transform is fast and its computational cost is the
same as that of TGUH. Important properties of the TGUW transform include
orthonormality conditional on the merging order, nonlinearity and “tail-
greediness”, and will be investigated in Section 2. The TGUW transform is the
first step of our proposed TrendSegment procedure, which involves four steps.

The detection of point anomalies has been widely studied in both time series
and machine learning literature and the reader is referred to Chandola et al.,
(2009) for an extensive review. Our framework is different from a model
recently studied by Fisch et al.,  (2018) in that our focus is on linear trend
changes and point anomalies, while they do not focus on trends but only on
point and collective anomalies with respect to a constant baseline
distribution.

The remainder of the article is organised as follows. Section 2 gives a full
description of the TrendSegment procedure and the relevant theoretical results
are presented in Section 3. The supporting simulation studies are described in
Section 4 and our methodology is illustrated in Section 5 through climate
datasets. The proofs of our main theoretical results are in Appendix A and
Section A of the supplementary materials, and the geometric interpretation of
the TGUW transformation can be found in Section D of the supplementary
materials. The TrendSegment procedure is implemented in the R package
trendsegmentR.

##  2 Methodology

###  2.1 Summary of TrendSegment

The TrendSegment procedure for estimating the number and the locations of
change-points includes four steps. We give a broad picture first and outline
details in later sections.

  1. 1.

TGUW transformation. Perform the TGUW transform; a bottom-up unbalanced
adaptive wavelet transformation of the input data $X_{1},\ldots,X_{T}$ by
recursively applying local conditionally orthonormal transformations. This
produces a data-adaptive multiscale decomposition of the data with $T-2$
detail-type coefficients and 2 smooth coefficients. The resulting
conditionally orthonormal transform of the data hopes to encode most of the
energy of the signal in only a few detail-type coefficients arising at coarse
levels. This sparse representation of the data justifies thresholding in the
next step.

  2. 2.

Thresholding. Set to zero those detail coefficients whose magnitude is smaller
than a pre-specified threshold as long as all the non-zero detail coefficients
are connected to each other in the tree structure. This step performs
“pruning” as a way of deciding the significance of the sparse representation
obtained in step 1.

  3. 3.

Inverse TGUW transformation. Obtain an initial estimate of $f_{t}$ by carrying
out the inverse TGUW transformation of the thresholded coefficient tree. The
resulting estimator can be shown to be $l_{2}$-consistent, but not yet
consistent for $N$ or $\eta_{1},\ldots,\eta_{N}$.

  4. 4.

Post-processing. Post-process the estimate from step 3 by removing some
change-points perceived to be spurious, which enables us to achieve estimation
consistency for $N$ and $\eta_{1},\ldots,\eta_{N}$.

We devote the following four sections to describing each step above in order.

###  2.2 TGUW transformation

####  2.2.1 Key principles of the TGUW transform

In the initial stage, the data are considered smooth coefficients and the TGUW
transform iteratively updates the sequence of smooth coefficients by merging
the adjacent sections of the data which are the most likely to belong to the
same segment. The merging is done by performing an adaptively constructed
orthonormal transformation to the chosen triplet of the smooth coefficients
and in doing so, a data-adaptive unbalanced wavelet basis is established. The
TGUW transform is completed after $T-2$ such orthonormal transformations and
each merge is performed under the following principles.

  1. 1.

In each merge, three adjacent smooth coefficients are selected and the
orthonormal transformation converts those three values into one detail and two
(updated) smooth coefficients. The size of the detail coefficient gives
information about the strength of the local linearity and the two updated
smooth coefficients are associated with the estimated parameters (intercept
and slope) of the local linear regression performed on the raw observations
corresponding to the initially chosen three smooth coefficients.

  2. 2.

“Two together” rule. The two smooth coefficients returned by the orthonormal
transformation are paired in the sense that both contain information about one
local linear regression fit. Thus, we require that any such pair of smooth
coefficients cannot be separated in any subsequent merges. We refer to this
recipe as the “two together” rule.

  3. 3.

To decide which triplet of smooth coefficients should be merged next, we
compare the corresponding detail coefficients as their magnitude represents
the strength of the corresponding local linear trend; the smaller the
(absolute) size of the detail, the smaller the local deviation from linearity.
Smooth coefficients corresponding to the smallest detail coefficients have
priority in merging.

As merging continues under the “two together” rule, all mergings can be
classified into one of three forms, Type 1: merging three initial smooth
coefficients, Type 2: merging one initial and a paired smooth coefficient and
Type 3: merging two sets of (paired) smooth coefficients (this is composed of
two merges of triplets; more details are given later).

Table 1: Notation. See Section 2.2.3 for formulae for the terms listed. $X_{p}$ |  $p^{\textnormal{th}}$ element of the observation vector $\boldsymbol{X}=\\{X_{1},X_{2},\ldots,X_{T}\\}^{\text{T}}$.
---|---
$s^{0}_{p,p}$ |  $p^{\textnormal{th}}$ initial smooth coefficient of the vector $\boldsymbol{s}^{0}$ where $\boldsymbol{X}=\boldsymbol{s}^{0}$.
$d_{p,q,r}$ |  detail coefficient obtained from $\\{X_{p},\ldots,X_{r}\\}$ (merges of Types 1 or 2).
$s^{1}_{p,r},s^{2}_{p,r}$ |  smooth coefficients obtained from $\\{X_{p},\ldots,X_{r}\\}$, paired under the “two together” rule.
$d^{1}_{p,q,r},d^{2}_{p,q,r}$ |  paired detail coefficients obtained by merging two adjacent subintervals, $\\{X_{p},\ldots,X_{q}\\}$ and $\\{X_{q+1},\ldots,X_{r}\\}$, where $r>q+2$ and $q>p+1$ (merge of Type 3).
$\boldsymbol{s}$ |  data sequence vector containing the (recursively updated) smooth and detail coefficients from the initial input $\boldsymbol{s}^{0}$.

####  2.2.2 Example

We now provide a simple example of the TGUW transformation; the accompanying
illustration is in Figure 1. The notation for this example and for the general
algorithm introduced later is in Table 1. This example shows single merges at
each pass through the data. We will later generalise it to multiple passes
through the data, which will speed up computation (this device is referred to
as “tail-greediness”). We refer to $j^{\text{th}}$ pass through the data as
scale $j$. Assume that we have the initial input
$\boldsymbol{s}^{0}=(X_{1},X_{2},\ldots,X_{8})$, so that the complete TGUW
transform consists of 6 merges. We show 6 example merges one by one under the
rules introduced in Section 2.2.1. This example demonstrates all three
possible types of merges.

Scale $j=1$. From the initial input $\boldsymbol{s}^{0}=(X_{1},\ldots,X_{8})$,
we consider 6 triplets $(X_{1},X_{2},X_{3})$, $(X_{2},X_{3},X_{4})$,
$(X_{3},X_{4},X_{5})$, $(X_{4},X_{5},X_{6})$, $(X_{5},X_{6},X_{7})$,
$(X_{6},X_{7},X_{8})$ and compute the size of the detail for each triplet,
where the formula can be found in (2). Suppose that $(X_{2},X_{3},X_{4})$
gives the smallest size of detail, $|d_{2,3,4}|$, then merge
$(X_{2},X_{3},X_{4})$ through the orthogonal transformation formulated in (4)
and update the data sequence into
$\boldsymbol{s}=(X_{1},s^{1}_{2,4},s^{2}_{2,4},d_{2,3,4},X_{5},X_{6},X_{7},X_{8})$.
We categorise this transformation into Type 1 (merging three initial smooth
coefficients).

Figure 1: Construction of tree for the example in Section 2.2.2; each diagram
shows all merges performed up to the given scale.

Scale $j=2$. From now on, the “two together” rule is applied. Ignoring any
detail coefficients in $\boldsymbol{s}$, the possible triplets for next
merging are $(X_{1},s^{1}_{2,4},s^{2}_{2,4})$,
$(s^{1}_{2,4},s^{2}_{2,4},X_{5})$, $(X_{5},X_{6},X_{7})$,
$(X_{6},X_{7},X_{8})$. We note that $(s^{2}_{2,4},X_{5},X_{6})$ cannot be
considered as a candidate for next merging under the two “together rule” as
this triplet contains only one (not both) of the paired smooth coefficients
returned by the previous merging. Assume that $(X_{5},X_{6},X_{7})$ gives the
smallest size of detail coefficient $|d_{5,6,7}|$ among the four candidates,
then we merge them through the orthogonal transformation formulated in (4) and
now update the sequence into
$\boldsymbol{s}=(X_{1},s^{1}_{2,4},s^{2}_{2,4},d_{2,3,4},s^{1}_{5,7},s^{2}_{5,7},d_{5,6,7},X_{8})$.
This transformation is also Type 1.

Scale $j=3$. We now compare four candidates for merging,
$(X_{1},s^{1}_{2,4},s^{2}_{2,4})$, $(s^{1}_{2,4},s^{2}_{2,4},s^{1}_{5,7})$,
$(s^{2}_{2,4},s^{1}_{5,7},s^{2}_{5,7})$ and $(s^{1}_{5,7},s^{2}_{5,7},X_{8})$.
The two triplets in middle, $(s^{1}_{2,4},s^{2}_{2,4},s^{1}_{5,7})$ and
$(s^{2}_{2,4},s^{1}_{5,7},s^{2}_{5,7})$, are paired together as they contain
two sets of paired smooth coefficients, $(s^{1}_{2,4},s^{2}_{2,4})$ and
$(s^{1}_{5,7},s^{2}_{5,7})$, and if we were to treat these two triplets
separately, we would be violating the “two together” rule. The summary detail
coefficient for this pair of triplets is obtained as
$d_{2,4,7}=\max(|d^{1}_{2,4,7}|,|d^{2}_{2,4,7}|)$, which is compared with
those of the other triplets. Now suppose that
$(X_{1},s^{1}_{2,4},s^{2}_{2,4})$ has the smallest size of detail; we merge
this triplet and update the data sequence into
$\boldsymbol{s}=(s^{1}_{1,4},s^{2}_{1,4},d_{1,1,4},d_{2,3,4},$
$s^{1}_{5,7},s^{2}_{5,7},d_{5,6,7},X_{8})$. This transformation is of Type 2.

Scale $j=4$. We now have two pairs of paired coefficients:
$(s^{1}_{1,4},s^{2}_{1,4})$ and $(s^{1}_{5,7},s^{2}_{5,7})$. Therefore, with
the “two together” rule in mind, the only possible options for merging are: to
merge the two pairs into $(s^{1}_{1,4},s^{2}_{1,4},s^{1}_{5,7},s^{2}_{5,7})$,
or to merge $(s^{1}_{5,7},s^{2}_{5,7})$ with $X_{8}$. Suppose that the first
merging is preferred. The merge of $(s^{1}_{1,4},s^{2}_{1,4})$ and
$(s^{1}_{5,7},s^{2}_{5,7})$ into
$(s^{1}_{1,4},s^{2}_{1,4},s^{1}_{5,7},s^{2}_{5,7})$ is of Type 3 and is
performed in two stages as follows. In the first stage, we merge
$(s^{1}_{1,4},s^{2}_{1,4},s^{1}_{5,7})$ and then update the sequence
temporarily as
$\boldsymbol{s}=(s^{1^{\prime}}_{1,7},s^{2^{\prime}}_{1,7},d_{1,1,4},d_{2,3,4},d^{1}_{1,4,7},s^{2}_{5,7},d_{5,6,7},X_{8})$.
In the second stage, we merge
$(s^{1^{\prime}}_{1,7},s^{2^{\prime}}_{1,7},s^{2}_{5,7})$, which gives the
updated sequence
$\boldsymbol{s}=(s^{1}_{1,7},s^{2}_{1,7},d_{1,1,4},d_{2,3,4},d^{1}_{1,4,7},d^{2}_{1,4,7},d_{5,6,7},$
$X_{8})$. As a summary detail coefficients for this merge, we use
$d_{1,4,7}=\max(|d^{1}_{1,4,7}|,|d^{2}_{1,4,7}|)$.

Scale $j=5$. The only available triplet is now
$(s^{1}_{1,7},s^{2}_{1,7},X_{8})$, thus we perform this Type 2 merge and
update the data sequence into
$\boldsymbol{s}=(s^{1}_{1,8},s^{2}_{1,8},d_{1,1,4},d_{2,3,4},d^{1}_{1,4,7},d^{2}_{1,4,7},d_{5,6,7},d_{1,7,8})$.
The transformation is completed with the updated data sequence which contains
$T-2=6$ detail and $2$ smooth coefficients.

####  2.2.3 TGUW transformation: general algorithm

In this Section, we formulate in generality the TGUW transformation
illustrated in Section 2.2.2. One of the important principles is “tail-
greediness” (Fryzlewicz, , 2018) which enables us to reduce the computational
complexity by performing multiple merges over non-overlapping regions in a
single pass over the data. More specifically, it allows us to perform up to
$\max\\{2,\lceil\rho\alpha_{j}\rceil\\}$ merges at each scale $j$, where
$\alpha_{j}$ is the number of smooth coefficients in the data sequence
$\boldsymbol{s}$ and $\rho\in(0,1)$ (the lower bound of 2 is essential to
permit a Type 3 transformation, which consists of two merges).

Sometimes, we will be referring to a detail coefficient $d_{p,q,r}^{\cdot}$ as
$d_{p,q,r}^{(j,k)}$ or $d^{(j,k)}$, where $j=1,\ldots,J$ is the scale of the
transform (i.e. the consecutive pass through the data) at which
$d_{p,q,r}^{\cdot}$ was computed, $k=1,\ldots,K(j)$ is the location index of
$d_{p,q,r}^{\cdot}$ within all scale $j$ coefficients, and $d_{p,q,r}^{\cdot}$
is $d_{p,q,r}^{1}$ or $d_{p,q,r}^{2}$ or $d_{p,q,r}$, depending on the type of
merge. We now describe the TGUW algorithm.

  1. 1.

At each scale $j$, find the set of triplets that are candidates for merging
under the “two together” rule and compute the corresponding detail
coefficients. Regardless of the type of merge, a detail coefficient
$d_{p,q,r}^{\cdot}$ is, in general, obtained as

$d_{p,q,r}^{\cdot}=a\boldsymbol{s}_{p:r}^{1}+b\boldsymbol{s}_{p:r}^{2}+c\boldsymbol{s}_{p:r}^{3},$
(2)

where $p\leq q<r$, $\boldsymbol{s}_{p:r}^{k}$ is the $k^{\text{th}}$ smooth
coefficient of the subvector $\boldsymbol{s}_{p:r}$ with a length of $r-p+1$
and the constants $a,b,c$ are the elements of the detail filter
$\boldsymbol{h}=(a,b,c)^{\text{T}}$. We note that $(a,b,c)$ also depends on
$(p,q,r)$, but this is not reflected in the notation, for simplicity. The
detail filter is a weight vector used in computing the weighted sum of a
triplet of smooth coefficients which should satisfy the condition that the
detail coefficient is zero if and only if the corresponding raw observations
over the merged regions have a perfect linear trend. If $(X_{p},\ldots,X_{r})$
are the raw observations associated with the triplet of the smooth
coefficients
$(\boldsymbol{s}_{p:r}^{1},\boldsymbol{s}_{p:r}^{2},\boldsymbol{s}_{p:r}^{3})$
under consideration, then the detail filter $\boldsymbol{h}$ is obtained in
such a way as to produce zero detail coefficient only when
$(X_{p},\ldots,X_{r})$ has a perfect linear trend, as the detail coefficient
itself represents the extent of non-linearity in the corresponding region of
data. This implies that the smaller the size of the detail coefficient, the
closer the alignment of the corresponding data section with linearity.
Specifically, the detail filter $\boldsymbol{h}=(a,b,c)^{\text{T}}$ is
established by solving the following equations,

$\displaystyle\begin{split}&a\boldsymbol{w}^{c,1}_{p:r}+b\boldsymbol{w}^{c,2}_{p:r}+c\boldsymbol{w}^{c,3}_{p:r}=0,\\\
&a\boldsymbol{w}^{l,1}_{p:r}+b\boldsymbol{w}^{l,2}_{p:r}+c\boldsymbol{w}^{l,3}_{p:r}=0,\\\
&a^{2}+b^{2}+c^{2}=1,\end{split}$ (3)

where $\boldsymbol{w}^{\cdot,k}_{p:r}$ is $k^{\text{th}}$ nonzero element of
the subvector $\boldsymbol{w}^{\cdot}_{p:r}$ with a length of $r-p+1$, and
$\boldsymbol{w}^{c}$ and $\boldsymbol{w}^{l}$ are weight vectors of constancy
and linearity, respectively, in which the initial inputs have a form of
${\boldsymbol{w}^{c}_{0}}=(1,1,\ldots,1)^{\text{T}},{\boldsymbol{w}^{l}_{0}}=(1,2,\ldots,T)^{\text{T}}$.
The last condition in (3) is to preserve the orthonormality of the transform.
The solution to (3) is unique up to multiplication by $-1$.

  2. 2.

Summarise all $d_{p,q,r}^{\cdot}$ constructed in step 1 to a (equal length or
shorter) sequence of $d_{p,q,r}$ by finding a summary detail coefficient
$d_{p,q,r}=\max(|d^{1}_{p,q,r}|,|d^{2}_{p,q,r}|)$ for any pair of detail
coefficients constructed by type 3 merges.

  3. 3.

Sort the size of the summarised detail coefficients $|d_{p,q,r}|$ obtained in
step 2 in non-decreasing order.

  4. 4.

Extract the (non-summarised) detail coefficient(s) $|d_{p,q,r}^{\cdot}|$
corresponding to the smallest (summarised) detail coefficient $|d_{p,q,r}|$
where both $|d^{1}_{p,q,r}|$ and $|d^{2}_{p,q,r}|$ should be extracted only if
$d_{p,q,r}=\max(|d^{1}_{p,q,r}|,|d^{2}_{p,q,r}|)$. Repeat the extraction until
$\max\\{2,\lceil\rho\alpha_{j}\rceil\\}$ (or all possible, whichever is the
smaller number) detail coefficients have been obtained, as long as the region
of the data corresponding to each detail coefficient extracted does not
overlap with the regions corresponding to the detail coefficients already
drawn.

  5. 5.

For each $|d_{p,q,r}^{\cdot}|$ extracted in step 4, merge the corresponding
smooth coefficients by updating the corresponding triplet in $\boldsymbol{s}$,
$\boldsymbol{w}^{c}$ and $\boldsymbol{w}^{l}$ through the orthonormal
transform as follows,

$\displaystyle\begin{pmatrix}s^{1}_{p,r}\\\ s^{2}_{p,r}\\\
d^{\cdot}_{p,q,r}\end{pmatrix}=\begin{pmatrix}&&\boldsymbol{\ell}_{1}^{\text{T}}&&\\\
&&\boldsymbol{\ell}_{2}^{\text{T}}&&\\\ &&\boldsymbol{h}^{\text{T}}&&\\\
\end{pmatrix}\begin{pmatrix}\boldsymbol{s}_{p:r}^{1}\\\
\boldsymbol{s}_{p:r}^{2}\\\
\boldsymbol{s}_{p:r}^{3}\end{pmatrix}=\Lambda\begin{pmatrix}\boldsymbol{s}_{p:r}^{1}\\\
\boldsymbol{s}_{p:r}^{2}\\\ \boldsymbol{s}_{p:r}^{3}\end{pmatrix},$ (4)
$\displaystyle\begin{pmatrix}{w}^{c,1}_{p,r}\\\ {w}^{c,2}_{p,r}\\\
0\end{pmatrix}=\Lambda\begin{pmatrix}\boldsymbol{w}^{c,1}_{p:r}\\\
\boldsymbol{w}^{c,2}_{p:r}\\\
\boldsymbol{w}^{c,3}_{p:r}\end{pmatrix},\quad\begin{pmatrix}{w}^{l,1}_{p,r}\\\
{w}^{l,2}_{p,r}\\\
0\end{pmatrix}=\Lambda\begin{pmatrix}\boldsymbol{w}^{l,1}_{p:r}\\\
\boldsymbol{w}^{l,2}_{p:r}\\\ \boldsymbol{w}^{l,3}_{p:r}\end{pmatrix}.$ (5)

The key step is finding the $3\times 3$ orthonormal matrix, $\Lambda$, which
is composed of one detail and two low-pass filter vectors in its rows. Firstly
the detail filter $\boldsymbol{h}^{\text{T}}$ is determined to satisfy the
conditions in (3), and then the two low-pass filters
($\boldsymbol{\ell}_{1}^{\text{T}},\boldsymbol{\ell}_{2}^{\text{T}}$) are
obtained by satisfying the orthonormality of $\Lambda$. There is no uniqueness
in the choice of
($\boldsymbol{\ell}_{1}^{\text{T}},\boldsymbol{\ell}_{2}^{\text{T}}$), but
this has no effect on the transformation itself. The details of this mechanism
can be found in Section D of the supplementary materials.

  6. 6.

Go to step 1 and repeat at new scale $j=j+1$ as long as we have at least three
smooth coefficients in the updated data sequence $\boldsymbol{s}$.

More specifically, the detail coefficient in (2) is formulated for each type
of merging introduced in Section 2.2.2 as follows.
Type 1: merging three initial smooth coefficients
$(s^{0}_{p,p},s^{0}_{p+1,p+1},s^{0}_{p+2,p+2})$,

$d_{p,p+1,p+2}=a_{p,p+1,p+2}s^{0}_{p,p}+b_{p,p+1,p+2}s^{0}_{p+1,p+1}+c_{p,p+1,p+2}s^{0}_{p+2,p+2}.$
(6)

Type 2: merging one initial and a paired smooth coefficient
$(s^{0}_{p,p},s^{1}_{p+1,r},s^{2}_{p+1,r})$,

$d_{p,p,r}=a_{p,p,r}s^{0}_{p,p}+b_{p,p,r}s^{1}_{p+1,r}+c_{p,p,r}s^{2}_{p+1,r},\quad\text{where}\quad
p+2<r,$ (7)

similarly, when merging a paired smooth coefficient and one initial,
$(s^{1}_{p,r-1},s^{2}_{p,r-1},s^{0}_{r,r})$,

$d_{p,r-1,r}=a_{p,r-1,r}s^{1}_{p,r-1}+b_{p,r-1,r}s^{2}_{p,r-1}+c_{p,r-1,r}s^{0}_{r,r},\quad\text{where}\quad
p+2<r.$ (8)

Type 3: merging two sets of (paired) smooth coefficients,
$(s^{1}_{p,q},s^{2}_{p,q})$ and $(s^{1}_{q+1,r},s^{2}_{q+1,r})$,

$\begin{aligned}
d^{1}_{p,q,r}=a^{1}_{p,q,r}s^{1}_{p,q}+b^{1}_{p,q,r}s^{2}_{p,q}+c^{1}_{p,q,r}s^{1}_{q+1,r}\\\
d^{2}_{p,q,r}=a^{2}_{p,q,r}s^{01}_{p,r}+b^{2}_{p,q,r}s^{02}_{p,r}+c^{2}_{p,q,r}s^{2}_{q+1,r}\end{aligned}\qquad\Longrightarrow\qquad\begin{aligned}
d_{p,q,r}=\max(|d^{1}_{p,q,r}|,|d^{2}_{p,q,r}|),\\\ \end{aligned}$ (9)

where $q>p+1$ and $r>q+2$. Importantly, the two consecutive merges in (9) are
achieved by visiting the same two adjacent data regions twice. In this case,
after the first detail coefficient, $d^{1}_{p,q,r}$, has been obtained, we
instantly update the corresponding triplets $\boldsymbol{s}$,
$\boldsymbol{w}^{c}$ and $\boldsymbol{w}^{l}$ via an orthonormal transform as
defined in (4) and (5). Therefore, the second detail filter,
$(a^{2}_{p,q,r},b^{2}_{p,q,r},c^{2}_{p,q,r})$, is constructed with the updated
$\boldsymbol{w}^{c}$ and $\boldsymbol{w}^{l}$ in a way that satisfies the
conditions (3).

The TGUW transform eventually converts the input data sequence
$\boldsymbol{X}$ of length $T$ into the sequence containing 2 smooth and $T-2$
detail coefficients through $T-2$ orthonormal transforms. The detail
coefficients $d^{(j,k)}$ can be regarded as scalar products between
$\boldsymbol{X}$ and a particular unbalanced wavelet basis $\psi^{(j,k)}$,
where the formal representation is given as $\\{d^{(j,k)}=\langle
X,\psi^{(j,k)}\rangle,_{j=1,\ldots,J,k=1,\;\ldots,K(j)}\\}$ for detail
coefficients and $s^{1}_{1,T}=\langle X,\psi^{(0,1)}\rangle$,
$s^{2}_{1,T}=\langle X,\psi^{(0,2)}\rangle$ for the two smooth coefficients.
The set $\\{\psi^{(j,k)}\\}$ is an orthonormal unbalanced wavelet basis for
$\mathbb{R}^{T}$. Some additional properties of the TGUW transform such as
sparse representation and computational complexity are discussed in Section
2.6.

###  2.3 Thresholding

Because at each stage, the TGUW transform constructs the smallest possible
detail coefficients, but it is at the same time orthonormal and so preserves
the $l_{2}$ energy of the input data, the variability (= deviation from
linearity) of the signal tends to be mainly encoded in only a few detail
coefficients computed at the later stages of the transform. The resulting
sparsity of representation of the input data in the domain of TGUW
coefficients justifies thresholding as a way of deciding the significance of
each detail coefficient (which measures the local deviation from linearity).

We propose to threshold the TGUW detail coefficients under two important
rules, which should simultaneously be satisfied; we refer to these as the
“connected” rule and the “two together” rule. The “connected” rule prunes the
branches of the TGUW detail coefficients if and only if the detail coefficient
itself and all of its children coefficients fall below a certain threshold in
absolute value. For instance, referring to the example of Section 2.2.2, if
both $d_{1,1,4}$ and $d_{1,7,8}$ were to survive the initial thresholding, the
“connected” rule would mean we also had to keep $d_{1,4,7}^{1}$ and
$d_{1,4,7}^{2}$, which are the children of $d_{1,7,8}$ and the parents of
$d_{1,1,4}$ in the TGUW coefficient tree.

The “two together” rule in thresholding is similar to the one in the TGUW
transformation except it targets pairs of detail rather than smooth
coefficients, and only applies to pairs of detail coefficients arising from
Type 3 merges. One such pair in the example of Section 2.2.2 is
$(d_{1,4,7}^{1},d_{1,4,7}^{2})$. The “two together” rule means that both such
detail coefficients should be kept if at least one survives the initial
thresholding. This is a natural requirement as a pair of Type 3 detail
coefficients effectively corresponds to a single merge of two adjacent
regions.

Through the thresholding, we wish to estimate the underlying signal $f$ in (1)
by estimating $\mu^{(j,k)}=\langle f,\psi^{(j,k)}\rangle$ where $\psi^{(j,k)}$
is an orthonormal unbalanced wavelet basis constructed in the TGUW transform
from the data. Throughout the thresholding procedure, the “connected” and “two
together” rules are applied in this order. We firstly threshold and apply the
“connected” rule, which gives us $\hat{\mu}_{0}^{(j,k)}$, the initial
estimator of $\mu^{(j,k)}$, as

$\hat{\mu}_{0}^{(j,k)}=d^{(j,k)}_{p,q,r}\cdot\mathbb{I}\;\Big{\\{}\,\exists(j^{\prime},k^{\prime})\in\mathcal{C}_{j,k}\quad\big{|}d^{(j^{\prime},k^{\prime})}_{p^{\prime},q^{\prime},r^{\prime}}\big{|}>\lambda\,\Big{\\}},$
(10)

where $\mathbb{I}$ is an indicator function and

$\mathcal{C}_{j,k}=\\{(j^{\prime},k^{\prime}),j^{\prime}=1,\ldots,j,k^{\prime}=1,\ldots,K(j^{\prime}):d^{(j^{\prime},k^{\prime})}_{p^{\prime},q^{\prime},r^{\prime}}\text{
is such that }[p^{\prime},r^{\prime}]\subseteq[p,r]\\}.$ (11)

Now the “two together” rule is applied to the initial estimators
$\hat{\mu}_{0}^{(j,k)}$ to obtain the final estimators $\hat{\mu}^{(j,k)}$. We
firstly note that two detail coefficients, $d^{(j,k)}_{p,q,r}$ and
$d^{(j^{\prime},k+1)}_{p^{\prime},q^{\prime},r^{\prime}}$ are called “paired”
when they are formed by Type 3 mergings and when
$(j,p,q,r)=(j^{\prime},p^{\prime},q^{\prime},r^{\prime})$. The “two together”
rule is formulated as below,

$\displaystyle\hat{\mu}^{(j,k)}=$ $\displaystyle\hat{\mu}_{0}^{(j,k)},$ if
$d_{p,q,r}^{(j,k)}$ is not paired,  $\displaystyle\hat{\mu}^{(j,k)}=$
$\displaystyle\hat{\mu}_{0}^{(j,k)},$ if $d_{p,q,r}^{(j,k)}$ is paired with
$d_{p,q,r}^{(j,k^{\prime})}$ and both $\hat{\mu}_{0}^{(j,k)}$ and
$\hat{\mu}_{0}^{(j,k^{\prime})}$ are zero or nonzero,
$\displaystyle\hat{\mu}^{(j,k)}=$ $\displaystyle d^{(j,k)},$ if
$d_{p,q,r}^{(j,k)}$ is paired with $d_{p,q,r}^{(j,k^{\prime})}$ and
$\hat{\mu}_{0}^{(j,k^{\prime})}\neq 0$ and $\hat{\mu}_{0}^{(j,k)}=0$. (12)

It is important to note that the application of the two rules ensures that
$\tilde{f}$ is a piecewise-linear function composed of best linear fits (in
the least-squares sense) for each estimated interval of linearity. As an
aside, we note that the number of survived detail coefficients does not
necessarily equal the number of change-points in $\tilde{f}$ as a pair of
detail coefficients arising from a Type 3 merge are associated with a single
change-point.

###  2.4 Inverse TGUW transformation

The estimator $\tilde{f}$ of the true signal $f$ in (1) is obtained by
inverting (= transposing) the orthonormal transformations in (4) in reverse
order to that in which they were originally performed. This inverse TGUW
transformation is referred to as $\text{TGUW}^{-1}$, and thus

$\tilde{f}=\text{TGUW}^{-1}\big{\\{}\;\hat{\mu}^{(j,k)},j=1,\ldots,J,\,k=1,\ldots,K(j)\;\|\;s^{1}_{1,T},s^{2}_{1,T}\;\big{\\}},$
(13)

where $\|$ denotes vector concatenation.

###  2.5 Post processing for consistency of change-point detection

As will be specified in Theorem 1 of Section 3, the piecewise-linear estimator
$\tilde{f}$ in (13) possibly overestimates the number of change-points. To
remove the spurious estimated change-points and to achieve the consistency of
the number and the locations of the estimated change-points, we borrow the
post-processing framework of Fryzlewicz,  (2018). Lin et al.,  (2017) show
that we can usually post-process $l_{2}$-consistent estimators in this way as
a fast enough $l_{2}$ error rate implies that each true change-point has an
estimator nearby. The post-processing methodology includes two stages, i)
execution of three steps, TGUW transform, thresholding and inverse TGUW
transform, again to the estimator $\tilde{f}$ in (13) and ii) examination of
regions containing only one estimated change-point to check for its
significance.

##### Stage 1.

We transform the estimated function $\tilde{f}$ in (13) with change-points
$(\tilde{\eta}_{1},\tilde{\eta}_{2},\ldots,\tilde{\eta}_{\tilde{N}})$ into a
new estimator $\tilde{\tilde{f}}$ with corresponding change-points
$(\tilde{\tilde{\eta}}_{1},\tilde{\tilde{\eta}}_{2},\ldots,\tilde{\tilde{\eta}}_{\tilde{\tilde{N}}})$.
Using $\tilde{f}$ in (13) as an input data sequence $\boldsymbol{s}$, we
perform the TGUW transform as presented in Section 2.2.3, but in a greedy
rather than tail-greedy way such that only one detail coefficient $d^{(j,1)}$
is produced at each scale $j$, and thus $K(j)=1$ for all $j$. We repeat to
produce detail coefficients until the first detail coefficient such that
$|d^{(j,1)}|>\lambda$ is obtained where $\lambda$ is the parameter used in the
thresholding procedure described in Section 2.3. Once the condition,
$|d^{(j,1)}|>\lambda$, is satisfied, stop merging and relabel the surviving
change-points as
$(\tilde{\tilde{\eta}}_{1},\tilde{\tilde{\eta}}_{2},\ldots,\tilde{\tilde{\eta}}_{\tilde{\tilde{N}}})$
and construct the new estimator $\tilde{\tilde{f}}$ as

$\tilde{\tilde{f}}_{t}=\hat{\theta}_{i,1}+\hat{\theta}_{i,2}\;t\quad\text{for}\quad
t\in\big{[}\tilde{\tilde{\eta}}_{i-1}+1,\tilde{\tilde{\eta}}_{i}\big{]},\quad
i=1,\ldots,\tilde{\tilde{N}},$ (14)

where $\tilde{\tilde{\eta}}_{0}=1$,
$\tilde{\tilde{\eta}}_{\tilde{\tilde{N}}+1}=T+1$ and
($\hat{\theta}_{i,1},\hat{\theta}_{i,2}$) are the OLS intercept and slope
coefficients, respectively, for the corresponding pairs
$\\{(t,X_{t}),\;t\in\big{[}\tilde{\tilde{\eta}}_{i-1}+1,\tilde{\tilde{\eta}}_{i}\big{]}\\}$.
The exception is when the region under consideration only contains a single
data point $X_{t_{0}}$ (a situation we refer to as a point anomaly throughout
the paper), in which case fitting a linear regression is impossible, so we
then set $\tilde{\tilde{f}}_{t_{0}}=X_{t_{0}}$.

##### Stage 2.

From the estimator $\tilde{\tilde{f}}_{t}$ in Stage 1, we obtain the final
estimator $\hat{f}$ by pruning the change-points
$(\tilde{\tilde{\eta}}_{1},\tilde{\tilde{\eta}}_{2},\ldots,\tilde{\tilde{\eta}}_{\tilde{\tilde{N}}})$
in $\tilde{\tilde{f}}_{t}$. For each $i=1,\ldots,\tilde{\tilde{N}}$, compute
the corresponding detail coefficient $d_{p_{i},q_{i},r_{i}}$ as described in
(7)-(9), where
$p_{i}=\Big{\lfloor}\frac{\tilde{\tilde{\eta}}_{i-1}+\tilde{\tilde{\eta}}_{i}}{2}\Big{\rfloor}+1$,
$q_{i}=\tilde{\tilde{\eta}}_{i}$ and
$r_{i}=\Big{\lceil}\frac{\tilde{\tilde{\eta}}_{i}+\tilde{\tilde{\eta}}_{i+1}}{2}\Big{\rceil}$.
Now prune by finding the minimiser
$i_{0}=\operatorname*{arg\,min}_{i}{|d_{p_{i},q_{i},r_{i}}|}$ and removing
$\tilde{\tilde{\eta}}_{i_{0}}$ and setting
$\tilde{\tilde{N}}:=\tilde{\tilde{N}}-1$ if
$|d_{p_{i_{0}},q_{i_{0}},r_{i_{0}}}|\leq\lambda$ where $\lambda$ is same as in
Section 2.3. Then relabel the change-points with the subscripts
$i=1,\ldots,\tilde{\tilde{N}}$ under the convention
$\tilde{\tilde{\eta}}_{0}=0$, $\tilde{\tilde{\eta}}_{\tilde{\tilde{N}}+1}=T$.
Repeat the pruning while we can find $i_{0}$ which satisfies the condition
$\big{|}d_{p_{i_{0}},q_{i_{0}},r_{i_{0}}}\big{|}<\lambda$. Otherwise, stop,
set $\hat{N}$ as the number of detected change-points and reconstruct the
change-points $\hat{\eta}_{i}$ in increasing order for $i=0,\ldots,\hat{N}+1$
where $\hat{\eta}_{0}=1$ and $\hat{\eta}_{\hat{N}+1}=T+1$. The estimated
function $\hat{f}$ is obtained by simple linear regression for each region
determined by the final change-points
$\hat{\eta}_{1},\ldots,\hat{\eta}_{\hat{N}}$ as in (14), with the exception
for point anomalies as described in Stage 1 above.

Through these two stages of post processing, the estimation of the number and
the locations of change-points becomes consistent, and further details can be
found in Section 3.

###  2.6 Extra discussion of TGUW transformation

##### Sparse representation.

The TGUW transform is nonlinear, but it is also linear and orthonormal
conditional on the order in which the merges are performed. The orthonormality
of the unbalanced wavelet basis, $\\{\psi^{(j,k)}\\}$, implies Parseval’s
identity,
$\sum_{t=1}^{T}X_{t}^{2}=\sum_{j=1}^{J}\sum_{k=1}^{K(j)}(d^{(j,k)})^{2}+(s^{1}_{1,T})^{2}+(s^{2}_{1,T})^{2}$
where $d^{(j,k)}=\langle X,\psi^{(j,k)}\rangle$, $s^{1}_{1,T}=\langle
X,\psi^{(0,1)}\rangle$ and $s^{2}_{1,T}=\langle X,\psi^{(0,2)}\rangle$.
Furthermore, the filters $(\psi^{(0,1)},\psi^{(0,2)})$ corresponding to the
two smooth coefficients $s^{1}_{1,T}$ and $s^{2}_{1,T}$ form an orthonormal
basis of the subspace
$\\{(x_{1},x_{2},\ldots,x_{T})\;|\;x_{1}-x_{2}=x_{2}-x_{3}=\cdots=x_{T-1}-x_{T}\\}$
of $\mathbb{R}^{T}$; see Section D of the supplementary materials for further
details. This implies that
$\sum_{t=1}^{T}X_{t}^{2}-(s^{1}_{1,T})^{2}-(s^{2}_{1,T})^{2}=\sum_{t=1}^{T}(X_{t}-\hat{X}_{t})^{2}$,
where $\hat{\boldsymbol{X}}=s^{1}_{1,T}\psi^{(0,1)}+s^{2}_{1,T}\psi^{(0,2)}$
is the best linear regression fit to $\boldsymbol{X}$ achieved by minimising
the sum of squared errors. This, combined with the Parseval’s identity above,
implies
$\sum_{t=1}^{T}(X_{t}-\hat{X}_{t})^{2}=\sum_{j=1}^{J}\sum_{k=1}^{K(j)}(d^{(j,k)})^{2}$.

By construction, the detail coefficients $|d^{(j,k)}|$ obtained in the initial
stages of the TGUW transform tend to be small in magnitude. Therefore, the
above Parseval’s identity implies that a large portion of
$\sum_{t=1}^{T}(X_{t}-\hat{X}_{t})^{2}$ is explained by only a few large
$|d^{(j,k)}|$’s arising in the later stages of the transform; in this sense,
the TGUW transform provides sparsity of signal representation.

##### Computational complexity.

Assume that $\alpha_{j}$ smooth coefficients are available in the data
sequence $\boldsymbol{s}$ at scale $j$. We allow the algorithm to merge up to
$\big{\lceil}\rho\alpha_{j}\big{\rceil}$ many triplets (unless their
corresponding data regions overlap) where $\rho\in(0,1)$ is a constant. This
gives us at most $(1-\rho)^{j}T$ smooth coefficients remaining in
$\boldsymbol{s}$ after $j$ scales. Solving for $(1-\rho)^{j}T\leq 2$ gives the
largest number of scales $J$ as
$\Big{\lceil}\log(T)/\log\big{(}(1-\rho)^{-1}\big{)}+\log(2)/\log(1-\rho)\Big{\rceil}$,
at which point the TGUW transform terminates with two smooth coefficients
remaining. Considering that the most expensive step at each scale is sorting
which takes $O(T\log(T))$ operations, the computational complexity of the TGUW
transformation is $O(T\log^{2}(T))$.

##  3 Theoretical results

We study the $l_{2}$ consistency of $\tilde{f}$ and $\tilde{\tilde{f}}$, and
the change-point detection consistency of $\hat{f}$, where the estimators are
defined in Section 2. The $l_{2}$ risk of an estimator $\tilde{f}$ is defined
as
$\big{\|}\tilde{f}-f\big{\|}_{T}^{2}=T^{-1}\sum_{i=1}^{T}(\tilde{f}_{i}-f_{i})^{2}$,
where $f$ is the underlying signal as in (1). We note the true change-points
$\\{\eta_{i},i=1,\ldots,N\\}$ are such that,

$\displaystyle\begin{split}&f_{t}=\theta_{\ell,1}+\theta_{\ell,2}\;t\;\text{
for }\;t\in[\eta_{\ell-1}+1,\eta_{\ell}],\;\ell=1,\ldots,N+1\\\ &\text{ where
}\;f_{\eta_{\ell}}+\theta_{\ell,2}\neq f_{\eta_{\ell}+1}\;\text{ for
}\;\ell=1,\ldots,N.\end{split}$ (15)

This definition permits both continuous and discontinuous changes and if
$f_{\eta_{i}}$ is a point anomaly, there exist two consecutive change-points
at $\eta_{i}-1$ and $\eta_{i}$ where $\eta_{i-1}=\eta_{i}-1$. We firstly
investigate the $l_{2}$ behaviour of $\tilde{f}$. The proofs of Theorems 1-3
can be found in Appendix A.

###### Theorem 1

$X_{t}$ follows model (1) with $\sigma=1$ and $\tilde{f}$ is the estimator in
(13). If the threshold $\lambda=C_{1}\\{2\log(T)\\}^{1/2}$ with a constant
$C_{1}$ large enough, then we have

$\mathbb{P}\;\bigg{(}\;\|\tilde{f}-f\|_{T}^{2}\;\leq\;C_{1}^{2}\;\frac{1}{T}\;\log(T)\;\Big{\\{}4+8N\;\lceil\;\log(T)/\log(1-\rho)^{-1}\;\rceil\;\Big{\\}}\;\bigg{)}\;\rightarrow\;1,$
(16)

as $T\rightarrow\infty$ and the piecewise-linear estimator $\tilde{f}$
contains $\tilde{N}\leq CN\log(T)$ change-points where $C$ is a constant.

Thus, $\tilde{f}$ is $l_{2}$ consistent under the strong sparsity assumption
i.e. if $N$ is finite. The crucial mechanism of $l_{2}$ consistency is the
“tail-greediness” which allows up to $K(j)\geq 1$ smooth coefficients to be
removed at each scale $j$. In other words, consistency is generally
unachievable if we proceed in a greedy (as opposed to tail-greedy) way, i.e.
if we only merge one triplet at each scale of the TGUW transformation.

We now move onto the estimator $\tilde{\tilde{f}}$ obtained in the first stage
of post-processing.

###### Theorem 2

$X_{t}$ follows model (1) with $\sigma=1$ and $\tilde{\tilde{f}}$ is the
estimator in (14). Let the threshold $\lambda$ be as in Theorem 1 and let the
number of true change-points, $N$, be finite. Then we have
$\big{\|}\tilde{\tilde{f}}-f\big{\|}_{T}^{2}\;=\;O\big{(}NT^{-1}\log^{2}(T)\big{)}$
with probability approaching $1$ as $T\rightarrow\infty$ and there exist at
most two estimated change-points between each pair of true change-points
$(\eta_{i},\eta_{i+1})$ for $i=0,\ldots,N$, where $\eta_{0}=0$ and
$\eta_{N+1}=T$. Therefore $\tilde{\tilde{N}}\leq 2(N+1)$.

We see that $\tilde{\tilde{f}}$ is $l_{2}$ consistent, but inconsistent for
the number of change-points. Now we investigate the final estimators,
$\hat{f}$ and $\hat{N}$.

###### Theorem 3

$X_{t}$ follows model (1) with $\sigma=1$ and ($\hat{f}$, $\hat{N}$) are the
estimators obtained in Section 2.5. Let the threshold $\lambda$ be as in
Theorem 1 and suppose that the number of true change-points, $N$, be finite.
Let
$\Delta_{T}=\min_{i=1,\ldots,N}\Big{\\{}\Big{(}\underaccent{\bar}{f}_{T}^{i}\Big{)}^{2/3}\cdot\delta_{T}^{i}\Big{\\}}$
where
$\underaccent{\bar}{f}_{T}^{i}=\min\Big{(}|f_{\eta_{i+1}}-2f_{\eta_{i}}+f_{\eta_{i-1}}|,|f_{\eta_{i+2}}-2f_{\eta_{i+1}}+f_{\eta_{i}}|\Big{)}$
and
$\delta_{T}^{i}=\min\Big{(}|\eta_{i}-\eta_{i-1}|,|\eta_{i+1}-\eta_{i}|\Big{)}$.
Assume that $T^{1/3}R_{T}^{1/3}=o\Big{(}\Delta_{T}\Big{)}$ where
$\big{\|}\tilde{\tilde{f}}-f\big{\|}_{T}^{2}=O_{p}(R_{T})$ is as in Theorem 2.
Then we have

$\mathbb{P}\;\bigg{(}\hat{N}=N,\quad\max_{i=1,\ldots,N}\bigg{\\{}|\hat{\eta}_{i}-\eta_{i}|\cdot\Big{(}\underaccent{\bar}{f}_{T}^{i}\Big{)}^{2/3}\bigg{\\}}\leq
CT^{1/3}R_{T}^{1/3}\bigg{)}\;\rightarrow\;1,$ (17)

as $T\rightarrow\infty$ where $C$ is a constant.

Our theory indicates that in the case in which
$\min_{i}\underaccent{\bar}{f}_{T}^{i}$ is bounded away from zero, the
consistent estimation of the number and locations of change-point is achieved
by assuming $T^{1/3}R_{T}^{1/3}=o(\delta_{T})$ where
$\delta_{T}=\min_{i=1,\ldots,N+1}|\eta_{i}-\eta_{i-1}|$. In addition, when
point anomalies exist in the set of true change-points, a point anomaly
$\eta_{k}$ and its neighbouring change-point $\eta_{k-1}=\eta_{k}-1$ can be
detected exactly at their true locations only if the corresponding
$\underaccent{\bar}{f}_{T}^{i}$s satisfy the condition
$\min\Big{(}\underaccent{\bar}{f}_{T}^{k},\underaccent{\bar}{f}_{T}^{k-1}\Big{)}\gtrsim\log(T)$.

##  4 Simulation study

###  4.1 Parameter choice and setting

_Post-processing._ In what follows, we disable Stages 1 and 2 of post-
processing by default: our empirical experience is that Stage 1 rarely makes a
difference in practice but comes with an additional computational cost, and
Stage 2 occasionally over-prunes change-point estimates.

_Choice of threshold $\lambda$._ Motivated by Theorem 1, we use the threshold
of the form $\lambda=C\sigma(2\log T)^{1/2}$ and estimate $\sigma$ using the
Median Absolute Deviation (MAD) estimator (Hampel, , 1974) defined as
$\hat{\sigma}=\text{Median}(|X_{1}-2X_{2}+X_{3}|,\ldots,|X_{T-2}-2X_{T-1}+X_{T}|)/(\Phi^{-1}(3/4)\sqrt{6})$
where $\Phi^{-1}$ is the quantile function of the Gaussian distribution. We
use $C=1.3$ as a default as it empirically led to the best performance over
the range $C\in[1,1.4]$.

_Choice of the “tail-greediness” parameter._ $\rho\in(0,1)$ is a constant
which controls the greediness level of the TGUW transformation in the sense
that it decides how many merges are performed in a single pass over the data.
A large $\rho$ can reduce the computational cost but it makes the procedure
less adaptive, whereas a small $\rho$ gives the opposite effect. Based on our
empirical experience, the best performance is achieved in the range
$\rho\in(0,0.05]$ and we use $\rho=0.04$ as a default in the simulation study
and data analyses.

We consider i.i.d. Gaussian noise and simulate data from model (1) using 8
signals, (M1) wave1, (M2) wave2, (M3) mix1, (M4) mix2, (M5) mix3, (M6)
lin.sgmts, (M7) teeth and (M8) lin, shown in Figure 2. (M1) is continuous at
change-points, while (M2) has discontinuities. (M3) has a mix of continuous
and discontinuous change-points and contains both constant and linear
segments, whereas (M4) is of the same type but also contains two point
anomalies. In addition, (M5) has two particularly short segments. (M6)
contains isolated spike-type short segments. (M7) is piecewise-constant, and
(M8) is a linear signal without change-points. The signals and R code for all
simulations can be downloaded from our GitHub repository (Maeng and
Fryzlewicz, , 2019) and the simulation results under dependent or heavy-tailed
errors can be found in Section B of the supplementary materials.

(a) (M1) wave1

(b) (M2) wave2

(c) (M3) mix1

(d) (M4) mix2

(e) (M5) mix3

(f) (M6) lin.sgmts

(g) (M7) teeth

(h) (M8) lin

Figure 2: Examples of data with their underlying signals studied in Section 4.
(a)-(h) data series $X_{t}$ (light grey) and true signal $f_{t}$ (black).

###  4.2 Competing methods and estimators

We perform the TrendSegment procedure based on the parameter choice in Section
4.1 and compare the performance with that of the following competitors:
Narrowest-Over-Threshold detection (NOT, Baranowski et al.,  (2016))
implemented in the R package not from CRAN, Isolate-Detect (ID, Anastasiou and
Fryzlewicz,  (2018)) available in the R package IDetect, trend filtering (TF,
Kim et al.,  (2009)) available from https://github.com/glmgen/genlasso,
Continuous-piecewise-linear Pruned Optimal Partitioning (CPOP, Maidstone et
al.,  (2017)) available from
https://www.maths.lancs.ac.uk/~fearnhea/Publications.html and a bottom-up
algorithm based on the residual sum of squares (RSS) from a linear fit (BUP,
Keogh et al.,  (2004)). The TrendSegment methodology is implemented in the R
package trendsegmentR.

As BUP requires a pre-specified number of change-points (or a well-chosen
stopping criterion which can vary depending on the data), we include it in the
simulation study (with the stopping criterion optimised for the best
performance using the knowledge of the truth) but not in data applications. We
do not include the methods of Spiriti et al.,  (2013) and Bai and Perron,
(2003) implemented in the R packages freeknotsplines and strucchange as we
have found them to be particularly slow. For instance, the minimum segment
size in strucchange can be adjusted to be small as long as it is greater than
or equal to 3 for detecting linear trend changes. This cannot capture point
anomalies but is suitable for detecting very short segments (e.g in (M6)
lin.sgmts). However, this setting is accompanied by extremely heavy
computation: with this minimum segment size constraint in place, a single
signal simulated from (M6) took us over three hours to process on a standard
PC.

Out of the competing methods tested, ID, TF and CPOP are in principle able to
classify two consecutive time point as change-points, and therefore they are
able to detect point anomalies. NOT and BUP are not designed to detect point
anomalies as their minimum distance between two consecutive change-points is
restricted to be at least two. For NOT, we use the contrast function for not
necessarily continuous piecewise-linear signals. Regarding the tuning
parameters for the competing methods, we follow the recommendation of each
respective paper or the corresponding R package.

Table 2: Distribution of $\hat{N}-N$ for models (M1)-(M4) and all methods listed in Section 4.1 and 4.2 over 100 simulation runs. Also the average MSE (Mean Squared Error) of the estimated signal $\hat{f}_{t}$ defined in Section 4.3, the average Hausdorff distance $d_{H}$ given by (18) and the average computational time in seconds using an Intel Core i5 2.9 GHz CPU with 8 GB of RAM, all over 100 simulations. Bold: methods within 10% of the highest empirical frequency of $\hat{N}-N=0$ or within 10% of the lowest empirical average $d_{H}(\times 10^{2})$. |  |  |  |  | $\hat{N}-N$ |  |  |  |  |  |
---|---|---|---|---|---|---|---|---|---|---|---
Model | Method |  $\leq$-3 | -2 | -1 | 0 | 1 | 2 |  $\geq$3 | MSE | $d_{H}(\times 10^{2})$ | time
(M1) | TS | 0 | 0 | 0 | 99 | 1 | 0 | 0 | 0.044 | 2.79 | 1.12
NOT | 0 | 0 | 0 | 99 | 1 | 0 | 0 | 0.034 | 2.09 | 0.29
ID | 0 | 0 | 0 | 99 | 1 | 0 | 0 | 0.029 | 1.45 | 0.22
TF | 0 | 0 | 0 | 0 | 0 | 0 | 100 | 0.016 | 4.29 | 36.30
CPOP | 0 | 0 | 0 | 99 | 1 | 0 | 0 | 0.014 | 0.78 | 8.55
BUP | 0 | 1 | 18 | 81 | 0 | 0 | 0 | 0.069 | 3.88 | 2.62
(M2) | TS | 0 | 0 | 2 | 98 | 0 | 0 | 0 | 0.109 | 1.90 | 1.06
NOT | 0 | 0 | 2 | 98 | 0 | 0 | 0 | 0.092 | 1.56 | 0.35
ID | 0 | 0 | 0 | 94 | 6 | 0 | 0 | 0.089 | 1.44 | 0.23
TF | 0 | 0 | 0 | 0 | 0 | 0 | 100 | 0.065 | 2.31 | 31.34
CPOP | 0 | 0 | 0 | 93 | 7 | 0 | 0 | 0.065 | 1.15 | 2.09
BUP | 100 | 0 | 0 | 0 | 0 | 0 | 0 | 0.752 | 4.69 | 2.21
(M3) | TS | 0 | 0 | 1 | 97 | 2 | 0 | 0 | 0.032 | 3.23 | 1.47
NOT | 0 | 0 | 0 | 100 | 0 | 0 | 0 | 0.020 | 2.35 | 0.36
ID | 0 | 0 | 1 | 94 | 5 | 0 | 0 | 0.047 | 2.37 | 0.33
TF | 0 | 0 | 0 | 0 | 0 | 0 | 100 | 0.023 | 5.87 | 45.31
CPOP | 0 | 0 | 0 | 61 | 32 | 6 | 1 | 0.024 | 2.34 | 21.11
BUP | 0 | 0 | 0 | 3 | 18 | 47 | 32 | 0.041 | 5.41 | 3.50
(M4) | TS | 0 | 0 | 5 | 76 | 18 | 1 | 0 | 0.030 | 1.81 | 1.48
NOT | 0 | 100 | 0 | 0 | 0 | 0 | 0 | 0.066 | 2.10 | 0.33
ID | 0 | 11 | 52 | 35 | 2 | 0 | 0 | 0.163 | 1.83 | 0.30
TF | 0 | 0 | 0 | 0 | 0 | 0 | 100 | 0.080 | 6.10 | 44.78
CPOP | 0 | 0 | 2 | 22 | 45 | 27 | 4 | 0.025 | 1.60 | 7.79
BUP | 0 | 0 | 8 | 31 | 45 | 13 | 3 | 0.092 | 5.30 | 3.62

Table 3: Distribution of $\hat{N}-N$ for models (M5)-(M8) and all methods listed in Section 4.1 and 4.2 over 100 simulation runs. Also the average MSE (Mean Squared Error) of the estimated signal $\hat{f}_{t}$ defined in Section 4.3, the average Hausdorff distance $d_{H}$ given by (18) and the average computational time in seconds using an Intel Core i5 2.9 GHz CPU with 8 GB of RAM, all over 100 simulations. Bold: methods within 10% of the highest empirical frequency of $\hat{N}-N=0$ or within 10% of the lowest empirical average $d_{H}(\times 10^{2})$. |  |  |  |  | $\hat{N}-N$ |  |  |  |  |  |
---|---|---|---|---|---|---|---|---|---|---|---
Model | Method |  $\leq$-3 | -2 | -1 | 0 | 1 | 2 |  $\geq$3 | MSE | $d_{H}(\times 10^{2})$ | time
(M5) | TS | 0 | 0 | 1 | 71 | 24 | 4 | 0 | 0.031 | 1.42 | 1.49
NOT | 0 | 0 | 99 | 1 | 0 | 0 | 0 | 0.040 | 1.20 | 0.29
ID | 0 | 0 | 1 | 2 | 14 | 32 | 51 | 0.277 | 8.28 | 0.30
TF | 0 | 0 | 0 | 0 | 0 | 0 | 100 | 0.116 | 6.17 | 43.13
CPOP | 0 | 0 | 0 | 11 | 22 | 39 | 28 | 0.023 | 1.41 | 5.12
BUP | 0 | 0 | 10 | 45 | 37 | 7 | 1 | 0.090 | 4.78 | 3.64
(M6) | TS | 0 | 0 | 0 | 96 | 4 | 0 | 0 | 0.013 | 0.05 | 1.65
NOT | 63 | 22 | 4 | 2 | 3 | 0 | 6 | 0.240 | 15.51 | 0.28
ID | 3 | 16 | 0 | 9 | 44 | 1 | 27 | 0.151 | 16.37 | 0.37
TF | 0 | 0 | 0 | 0 | 0 | 0 | 100 | 0.134 | 10.98 | 48.19
CPOP | 0 | 0 | 0 | 20 | 41 | 24 | 15 | 0.034 | 0.13 | 5.11
BUP | 0 | 0 | 0 | 0 | 0 | 0 | 100 | 0.135 | 10.17 | 4.00
(M7) | TS | 0 | 5 | 21 | 40 | 28 | 6 | 0 | 0.119 | 7.02 | 0.65
NOT | 1 | 1 | 8 | 56 | 31 | 3 | 0 | 0.065 | 2.62 | 0.25
ID | 3 | 0 | 16 | 14 | 26 | 13 | 28 | 0.320 | 10.87 | 0.12
TF | 0 | 0 | 0 | 0 | 0 | 0 | 100 | 0.097 | 6.11 | 23.19
CPOP | 0 | 0 | 1 | 1 | 3 | 17 | 78 | 0.055 | 3.37 | 1.19
BUP | 70 | 25 | 5 | 0 | 0 | 0 | 0 | 0.277 | 11.89 | 1.58
(M8) | TS | 0 | 0 | 0 | 100 | 0 | 0 | 0 | 0.001 | 0.00 | 1.01
NOT | 0 | 0 | 0 | 100 | 0 | 0 | 0 | 0.001 | 0.00 | 0.17
ID | 0 | 0 | 0 | 100 | 0 | 0 | 0 | 0.001 | 0.00 | 0.59
TF | 0 | 0 | 0 | 78 | 5 | 2 | 15 | 0.002 | 9.08 | 35.79
CPOP | 0 | 0 | 0 | 100 | 0 | 0 | 0 | 0.001 | 0.00 | 12.96
BUP | 0 | 0 | 0 | 0 | 0 | 0 | 100 | 0.011 | 46.34 | 2.63

###  4.3 Results

The summary of the results for all models and methods can be found in Tables 2
and 3. We run 100 simulations and as a measure of accuracy of estimators, we
use Monte-Carlo estimates of the Mean Squared Error of the estimated signal
defined as MSE=$\mathbb{E}\\{(1/T)\sum_{t=1}^{T}(f_{t}-\hat{f}_{t})^{2}\\}$.
The empirical distribution of $\hat{N}-N$ is also reported where $\hat{N}$ is
the estimated number of change-points and $N$ is the true one. In addition to
this, for comparing the accuracy of the locations of the estimated change-
points $\hat{\eta}_{i}$, we show estimates of the scaled Hausdorff distance
given by

$d_{H}=\frac{1}{T}\mathbb{E}\max\bigg{\\{}\max_{i}\min_{j}\big{|}\eta_{i}-\hat{\eta}_{j}\big{|},\quad\max_{j}\min_{i}\big{|}\hat{\eta}_{j}-\eta_{i}\big{|}\bigg{\\}},$
(18)

where $i=0,\ldots,N+1$ and $j=0,\ldots,\hat{N}+1$ with the convention
$\eta_{0}=\hat{\eta}_{0}=0,\eta_{N+1}=\hat{\eta}_{N+1}=T$ and $\hat{\eta}$ and
$\eta$ denote estimated and true locations of the change-points. The smaller
the Hausdorff distance, the better the estimation of the change-point
locations. For each method, the average computation time in seconds is shown.

The results for (M1) and (M2) are similar. TrendSegment shows comparable
performance to NOT, ID and CPOP in terms of the estimation of the number of
change-points, while it is slightly less attractive in terms of the estimated
locations of change-points. TF tends to overestimate the number of change-
points throughout all models. When the signal is a mix of constant and linear
trends as in (M3), TrendSegment, NOT and ID still perform well in terms of the
estimation of the number of change-points, while CPOP tends to overestimate.
We see that TrendSegment has a particular advantage over the other methods
especially in (M4) and (M5), when point anomalies exist or in the case of
frequent change-points. TrendSegment shows its relative robustness in
estimating the number and the location of change-points while ID and CPOP
significantly underperform and NOT ignores the point anomalies, as expected.
(M6) is another example where only TrendSegment exhibits good performance. For
the estimation of the piecewise-constant signal (M7), no method performs well
and NOT, ID and TrendSegment tend to underestimate the number of change-points
while CPOP and TF overestimate. In the case of the no-change-point signal
(M8), all methods except TF perform well.

In summary, TrendSegment is always among the best methods, and is particularly
attractive for signals with point anomalies or short segments. With respect to
computation time, NOT and ID are very fast in all cases, TrendSegment is
slower than these two but is faster than TF, CPOP and BUP, especially when the
length of the time series is larger than 2000.

##  5 Data applications

###  5.1 Average January temperatures in Iceland

We analyse a land temperature dataset available from http://berkeleyearth.org,
consisting of average temperatures in January recorded in Reykjavik recorded
from 1763 to 2013. Figure 3(a) shows the data; the point corresponding to 1918
appears to be a point anomaly, and we comment on this aspect later on in this
section.

(a) data

(b) TrendSegment

(c) NOT and ID

(d) TF and CPOP

Figure 3: Change-point analysis for January average temperature in Reykjavik
from 1763 to 2013 in Section 5.1. (a) the data series, (b) the data series
(grey dots) and estimated signal with change-points returned by TrendSegment(
), (c) estimated signal with change-points returned by NOT ( ) and ID ( ), (d)
estimated signal with change-points returned by TF ( ) and CPOP ( ).

The TrendSegment estimate of the piecewise-linear trend is shown in Figure
3(b). It identifies 2 change-points, 1917 and 1918, where the temperature in
1918 is fitted as a single point as it is much lower than in other years.
Figures 3(c) and 3(d) show that NOT and CPOP detect the change of slope in
1974, ID returns an increasing function with no change-points and TF reports 6
points with the most recent one in 1981, but none of them detect the point
anomaly. To assess the goodness of fit of the TrendSegment estimate, we
computed the sample autocorrelation and partial autocorrelation functions of
the empirical residuals from the TrendSegment fit and both were consistent
with white noise.

Regarding the 1918 observation, Moore and Babij,  (2017) report that “[t]he
winter of 1917/1918 is referred to as the Great Frost Winter in Iceland. It
was the coldest winter in the region during the twentieth century. It was
remarkable for the presence of sea ice in Reykjavik Harbour as well as for the
unusually large number of polar bear sightings in northern Iceland.” This
example illustrates the flexibility of the TrendSegment as it detects not only
change-points in linear trend but it can identify a point anomaly at the same
time, which the competing methods do not achieve.

###  5.2 Monthly average sea ice extent of Arctic and Antarctic

We analyse the average sea ice extent of the Arctic and the Antarctic
available from https://nsidc.org to estimate the change-points in its trend.
As mentioned in Serreze and Meier,  (2018), sea ice extent is the most common
measure for assessing the condition of high-latitude oceans and it is defined
as the area covered with an ice concentration of at least 15$\%$. Here we use
the average ice extent in February and September as it is known that the
Arctic has the maximum ice extent typically in February while the minimum
occurs in September and the Antarctic experiences the opposite.

Serreze and Meier,  (2018) indicate that the clear decreasing trend of sea ice
extent of the Arctic in September is one of the most important indicator of
climate change. In contrast to the Arctic, the sea ice extent of the Antarctic
has been known to be stable in the sense that it shows a weak increasing trend
in the decades preceding 2016 (Comiso et al., , 2017; Serreze and Meier, ,
2018). However, Rintoul et al.,  (2018) warn of a possible collapse of the
past stability by citing a significant decline of the sea ice extent in 2016.
We now use the most up-to-date records (to 2018) and re-examine the concerns
expressed in Rintoul et al.,  (2018) with the help of our change-point
detection methodology.

Figures 4(a) and 4(c) show the well-known decreasing trend of the average sea
ice extent in the Arctic both in its winter (February) and summer (September).
In Figure 4(a), the TrendSegment estimate identifies change-points in 2004 and
2007 and detects a sudden drop during 2005-2007. One change-point in 2006 is
identified in Figure 4(c), which differentiates the decreasing speed of winter
ice extent in the Arctic before and after 2006. As observed in the above-
mentioned literature, the sea ice extent of the Antarctic shows a modest
increasing trend up until recently (Figures 4(b) and 4(d)); however, we
observe a strong decreasing trend from the detected change-point in 2016 for
the Antarctic summer (February) and from 2015 for the Antarctic winter
(September), which is in line with the message of Rintoul et al.,  (2018). The
results for the other competing methods can be found in Section C of the
supplementary materials.

(a) Arctic in February

(b) Antarctic in February

(c) Arctic in September

(d) Antarctic in September

Figure 4: The TrendSegment estimate of piecewise-linear trend for the monthly
average sea ice extent from 1979 to 2018 in Section 5.2. (a) the data series
(grey dots); the TrendSegment estimate ( ) for average sea ice extent of the
Arctic in February, (b) Antarctic in February, (c) Arctic in September, (d)
Antarctic in September.

##  Appendix A Technical proofs

The proof of Theorem 1-3 are below and Lemmas 1 and 2 can be found in Section
A of the supplementary materials.

Proof of Theorem 1. Let $\mathcal{S}^{1}_{j}$ and $\mathcal{S}^{0}_{j}$ be as
in Lemma 2. From the conditional orthonormality of the unbalanced wavelet
transform, on the set $A_{T}$ defined in Lemma 1, we have

$\displaystyle\|\tilde{f}-f\|_{T}^{2}\;=\;$
$\displaystyle\frac{1}{T}\sum_{j=1}^{J}\sum_{k=1}^{K(j)}\Big{(}d^{(j,k)}\cdot\mathbb{I}\big{\\{}\,\exists(j^{\prime},k^{\prime})\in\mathcal{C}_{j,k}\quad|d^{(j^{\prime},k^{\prime})}|>\lambda\,\big{\\}}-\mu^{(j,k)}\Big{)}^{2}+\;T^{-1}(s^{1}_{1,T}-\mu^{(0,1)})^{2}+\;T^{-1}(s^{2}_{1,T}-\mu^{(0,2)})^{2}$
$\displaystyle\;\leq\;$
$\displaystyle\frac{1}{T}\sum_{j=1}^{J}\Bigg{(}\sum_{k\in\mathcal{S}^{0}_{j}}+\sum_{k\in\mathcal{S}^{1}_{j}}\Bigg{)}\Big{(}d^{(j,k)}\cdot\mathbb{I}\big{\\{}\,\exists(j^{\prime},k^{\prime})\in\mathcal{C}_{j,k}\quad|d^{(j^{\prime},k^{\prime})}|>\lambda\,\big{\\}}-\mu^{(j,k)}\Big{)}^{2}+\;4C_{1}^{2}T^{-1}\log
T$ $\displaystyle\;=:$
$\displaystyle\;\mathit{I}+\mathit{II}+4C_{1}^{2}T^{-1}\log T.$ (19)

where $\mu^{(0,1)}=\langle f,\psi^{(0,1)}\rangle$ and $\mu^{(0,2)}=\langle
f,\psi^{(0,2)}\rangle$. For $j=1,\ldots,J,k\in\mathcal{S}^{0}_{j}$, we have
$|d^{(j,k)}|\leq\lambda$, where $\lambda$ is as in Theorem 1. By Lemma 2,
$\mathbb{I}\big{\\{}\,\exists(j^{\prime},k^{\prime})\in\mathcal{C}_{j,k}\quad|d^{(j^{\prime},k^{\prime})}|>\lambda\big{\\}}=0$
for $k\in\mathcal{S}^{0}_{j}$ and also by the fact that $\mu^{(j,k)}=0$ for
$j=1,\ldots,J,k\in\mathcal{S}^{0}_{j}$, we have $\mathit{I}=0$. For
$\mathit{II}$, we denote
$\mathcal{B}=\big{\\{}\,\exists(j^{\prime},k^{\prime})\in\mathcal{C}_{j,k}\quad|d^{(j^{\prime},k^{\prime})}|>\lambda\,\big{\\}}$
and have

$\displaystyle\big{(}d^{(j,k)}\cdot\mathbb{I}\big{\\{}\mathcal{B}\big{\\}}-\mu^{(j,k)}\big{)}^{2}\;=\;$
$\displaystyle\big{(}d^{(j,k)}\cdot\mathbb{I}\big{\\{}\mathcal{B}\big{\\}}-d^{(j,k)}+d^{(j,k)}-\mu^{(j,k)}\big{)}^{2}$
(20) $\displaystyle\;\leq\;$
$\displaystyle\big{(}d^{(j,k)}\big{)}^{2}\mathbb{I}\big{(}|d^{(j^{\prime},k^{\prime})}|\leq\lambda\big{)}+2|d^{(j,k)}|\;\mathbb{I}\big{(}|d^{(j^{\prime},k^{\prime})}|\leq\lambda\big{)}\;|d^{(j,k)}-\mu^{(j,k)}|+\;\big{(}d^{(j,k)}-\mu^{(j,k)}\big{)}^{2}$
$\displaystyle\;\leq\;$ $\displaystyle\lambda^{2}+2\lambda C_{1}\\{2\log
T\\}^{1/2}+2C_{1}^{2}\log T.$

Combining with the upper bound of $J$,
$\lceil\log(T)/\log(1-\rho)^{-1}\rceil$, and the fact that
$|\mathcal{S}^{1}_{j}|\leq N$, we have $\mathit{II}\leq
8C_{1}^{2}NT^{-1}\lceil\log(T)/\log(1-\rho)^{-1}\rceil\log T$, and therefore
$\|\tilde{f}-f\|_{T}^{2}\;\leq\;C_{1}^{2}\;\frac{1}{T}\;\log(T)\;\Big{\\{}4+8N\;\lceil\;\log(T)/\log(1-\rho)^{-1}\;\rceil\;\Big{\\}}$.
Also, at each scale, the estimated change-points are obtained up to size $N$,
combining it with the largest scale $J$, the number of change-points in
$\tilde{f}$ returned from the inverse TGUW transformation is up to
$CN\text{log}T$ where $C$ is a constant.

Proof of Theorem 2. Let $\tilde{B}$ and $\tilde{\tilde{B}}$ the unbalanced
wavelet bases corresponding to $\tilde{f}$ and $\tilde{\tilde{f}}$,
respectively. As the change-points in $\tilde{\tilde{f}}$ are a subset of
those in $\tilde{f}$, establishing $\tilde{\tilde{f}}$ can be regarded as
applying the TGUW transform again to $\tilde{f}$, which is just a repetition
of the estimation procedure $\tilde{f}$ but performed in a greedy way. Thus
$\tilde{\tilde{B}}$ is classified into two categories, 1) all basis vectors
$\psi^{(j,k)}\in\tilde{B}$ such that $\psi^{(j,k)}$ is not associated with the
change-points in $\tilde{f}$ and $|\langle
X,\psi^{(j,k)}\rangle|=|d^{(j,k)}|<\lambda$ and 2) all vectors $\psi^{(j,1)}$
produced in Stage 1 of post-processing.

We now investigate how many scales are used for this particular transform.
Firstly, the detail coefficients $d^{(j,k)}$ corresponding to the basis
vectors $\psi^{(j,k)}\in\tilde{B}$ live on no more than $J=O(\log T)$ scales
and we have $|\mathcal{S}^{1}_{j}|\leq N$ by the argument used in the proof of
Theorem 1. In addition, the vectors $\psi^{(j,1)}$ in the second category
above correspond to different change-points in $\tilde{f}$ and there exist at
most $\tilde{N}=O(N\text{log}T)$ change-points in $\tilde{f}$ which we examine
one at once (i.e. $|\mathcal{S}^{1}_{j}|\leq 1$), thus at most $\tilde{N}$
scales are required for $d^{(j,1)}$. Combining the results of the two
categories, the equivalent of quantity $\mathit{II}$ in the proof of Theorem 1
for $\tilde{\tilde{f}}$ is bounded by $\mathit{II}\leq C_{3}NT^{-1}\log^{2}T$
and this completes the proof of the $L_{2}$ result,
$\big{\|}\tilde{\tilde{f}}-f\big{\|}_{T}^{2}\;=\;O\big{(}NT^{-1}\log^{2}(T)\big{)}$
where $C_{3}$ is a large enough positive constant.

Finally, we show that there exist at most two change-points in
$\tilde{\tilde{f}}$ between true change points $(\eta_{i},\eta_{i+1})$ for
$i=0,\ldots,N$ where $\eta_{0}=0$ and $\eta_{N+1}=T$. Consider the case where
three change-point for instance
($\tilde{\tilde{\eta}}_{l},\tilde{\tilde{\eta}}_{l+1},\tilde{\tilde{\eta}}_{l+2}$)
lie between a pair of true change-points, $(\eta_{i},\eta_{i+1})$. In this
case, by Lemma 2, the maximum magnitude of two detail coefficients computed
from the adjacent intervals,
$[\tilde{\tilde{\eta}}_{l}+1,\tilde{\tilde{\eta}}_{l+1}]$ and
$[\tilde{\tilde{\eta}}_{l+1}+1,\tilde{\tilde{\eta}}_{l+2}]$, is less than
$\lambda$ and $\tilde{\tilde{\eta}}_{l+1}$ would get removed from the set of
estimated change-points. This leads to $\tilde{\tilde{N}}\leq 2(N+1)$.

Proof of Theorem 3. From the assumptions of Theorem 3, 1) given any
$\epsilon>0$ and $C>0$, for some $T_{1}$ and all $T>T_{1}$, it holds that
$\mathbb{P}\Big{(}\big{\|}\tilde{\tilde{f}}-f\big{\|}_{T}^{2}>\frac{C^{3}}{4}R_{T}\Big{)}\leq\epsilon$
where $\tilde{\tilde{f}}$ is the estimated signal specified in Theorem 2 and
2) For some $T_{2}$, and all $T>T_{2}$, it holds that
$C^{1/3}T^{1/3}R_{T}^{1/3}(\Delta^{f}_{i})^{-2/3}<\Delta^{\eta}_{i}$ for all
$i=1,\ldots,N$. Similar to the argument of Theorem 19 in Lin et al.,  (2016),
we take $T\geq T^{*}$ where $T^{*}=\max\\{T_{1},T_{2}\\}$ and let
$r_{i,T}=\lfloor C^{1/3}T^{1/3}R_{T}^{1/3}(\Delta^{f}_{i})^{-2/3}\rfloor$ for
$i=1,\ldots,N$. Suppose that there exist at least one $\eta_{i}$ whose closest
estimated change-point is not within the distance of $r_{i,T}$. Then there are
no estimated change-points in $\tilde{\tilde{f}}$ within $r_{i,T}$ of
$\eta_{i}$ which means that $\tilde{\tilde{f}}_{j}$ displays a linear trend
over the entire segment $j\in\\{\eta_{i}-r_{i,T},\ldots,\eta_{i}+r_{i,T}\\}$.
Hence,
$\frac{1}{T}\sum_{j=\eta_{i}-r_{i,T}}^{\eta_{i}+r_{i,T}}(\tilde{\tilde{f}}_{j}-f_{j})^{2}\geq\frac{13r_{i,T}^{3}}{24T}\bigg{(}\Delta^{f}_{i}\bigg{)}^{2}>\frac{C^{3}}{4}R_{T}$.
We see that assuming that at least one $\eta_{i}$ does not have any estimated
change-point within the distance of $r_{i,T}$ implies the estimation error
exceeds $\frac{C^{3}}{4}R_{T}$ which is a contradiction as it is an event that
we know occurs with probability at most $\epsilon$. Therefore, there must
exist at least one estimated change-point within the distance of $r_{i,T}$
from each true change point $\eta_{i}$.

Throughout Stage 2 of post processing, $\tilde{\tilde{\eta}}_{i_{0}}$ is
either the closest estimated change-point of any $\eta_{i}$ or not. If
$\tilde{\tilde{\eta}}_{i_{0}}$ is not the closest estimated change-point to
the nearest true change-point on either its left or its right, by the
construction of detail coefficients in Stage 2 of post processing, Lemma 2
guarantees that the corresponding detail coefficient has the magnitude less
than $\lambda$ and $\tilde{\tilde{\eta}}_{i_{0}}$ gets removed. Suppose
$\tilde{\tilde{\eta}}_{i_{0}}$ is the closest estimated change-point of a true
change-point $\eta_{i}$ and it is within the distance of
$CT^{1/3}R_{T}^{1/3}\Big{(}\Delta^{f}_{i}\Big{)}^{-2/3}$ from $\eta_{i}$. If
the corresponding detail coefficient has the magnitude less than $\lambda$ and
$\tilde{\tilde{\eta}}_{i_{0}}$ is removed, there must exist another
$\tilde{\tilde{\eta}}_{i}$ within the distance of
$CT^{1/3}R_{T}^{1/3}\Big{(}\Delta^{f}_{i}\Big{)}^{-2/3}$ from $\eta_{i}$. If
there are no such $\tilde{\tilde{\eta}}_{i}$, then by the construction of the
detail coefficient, the order of magnitude of
$\big{|}d_{p_{i_{0}},q_{i_{0}},r_{i_{0}}}\big{|}$ would be such that
$\big{|}d_{p_{i_{0}},q_{i_{0}},r_{i_{0}}}\big{|}>\lambda$ thus
$\tilde{\tilde{\eta}}_{i_{0}}$ would not get removed. Therefore, after Stage 2
of post processing is finished, each true change-point $\eta_{i}$ has its
unique estimator within the distance of
$CT^{1/3}R_{T}^{1/3}\Big{(}\Delta^{f}_{i}\Big{)}^{-2/3}$.

## References

  * Anastasiou and Fryzlewicz, (2018) Anastasiou, A. and Fryzlewicz, P. (2018).  Detecting multiple generalized change-points by isolating single ones.  Preprint.
  * Bai and Perron, (1998) Bai, J. and Perron, P. (1998).  Estimating and testing linear models with multiple structural changes.  Econometrica, 66:47–78.
  * Bai and Perron, (2003) Bai, J. and Perron, P. (2003).  Computation and analysis of multiple structural change models.  Journal of applied econometrics, 18:1–22.
  * Baranowski et al., (2016) Baranowski, R., Chen, Y., and Fryzlewicz, P. (2016).  Narrowest-over-threshold detection of multiple change-points and change-point-like features.  arXiv preprint arXiv:1609.00293.
  * Bardwell et al., (2017) Bardwell, L., Fearnhead, P., et al. (2017).  Bayesian detection of abnormal segments in multiple time series.  Bayesian Analysis, 12:193–218.
  * Chandola et al., (2009) Chandola, V., Banerjee, A., and Kumar, V. (2009).  Anomaly detection: A survey.  ACM computing surveys (CSUR), 41:Article 15.
  * Comiso et al., (2017) Comiso, J. C., Gersten, R. A., Stock, L. V., Turner, J., Perez, G. J., and Cho, K. (2017).  Positive trend in the antarctic sea ice cover and associated changes in surface temperature.  Journal of Climate, 30:2251–2267.
  * Fisch et al., (2018) Fisch, A. T. M., Eckley, I. A., and Fearnhead, P. (2018).  A linear time method for the detection of point and collective anomalies.  arXiv preprint arXiv:1806.01947.
  * Fryzlewicz, (2018) Fryzlewicz, P. (2018).  Tail-greedy bottom-up data decompositions and fast mulitple change-point detection.  The Annals of Statistics, 46:3390–3421.
  * Hampel, (1974) Hampel, F. R. (1974).  The influence curve and its role in robust estimation.  Journal of the american statistical association, 69:383–393.
  * Jamali et al., (2015) Jamali, S., Jönsson, P., Eklundh, L., Ardö, J., and Seaquist, J. (2015).  Detecting changes in vegetation trends using time series segmentation.  Remote Sensing of Environment, 156:182–195.
  * James et al., (2016) James, N. A., Kejariwal, A., and Matteson, D. S. (2016).  Leveraging cloud data to mitigate user experience from ‘breaking bad’.  In Big Data (Big Data), 2016 IEEE International Conference on, pages 3499–3508. IEEE.
  * Jeng et al., (2012) Jeng, X. J., Cai, T. T., and Li, H. (2012).  Simultaneous discovery of rare and common segment variants.  Biometrika, 100:157–172.
  * Keogh et al., (2004) Keogh, E., Chu, S., Hart, D., and Pazzani, M. (2004).  Segmenting time series: A survey and novel approach.  In Data mining in time series databases, pages 1–21. World Scientific.
  * Kim et al., (2009) Kim, S.-J., Koh, K., Boyd, S., and Gorinevsky, D. (2009).  $\ell_{1}$ trend filtering.  SIAM review, 51:339–360.
  * Lin et al., (2016) Lin, K., Sharpnack, J., Rinaldo, A., and Tibshirani, R. J. (2016).  Approximate recovery in changepoint problems, from $\ell_{2}$ estimation error rates.  arXiv preprint arXiv:1606.06746.
  * Lin et al., (2017) Lin, K., Sharpnack, J. L., Rinaldo, A., and Tibshirani, R. J. (2017).  A sharp error analysis for the fused lasso, with application to approximate changepoint screening.  In Advances in Neural Information Processing Systems, pages 6884–6893.
  * Maeng and Fryzlewicz, (2019) Maeng, H. and Fryzlewicz, P. (2019).  Detecting linear trend changes and point anomalies in data sequences: Simulation code.  URL https://github.com/hmaeng/trendsegment.
  * Maidstone et al., (2017) Maidstone, R., Fearnhead, P., and Letchford, A. (2017).  Detecting changes in slope with an $l\\_0$ penalty.  arXiv preprint arXiv:1701.01672.
  * Matteson and James, (2014) Matteson, D. S. and James, N. A. (2014).  A nonparametric approach for multiple change point analysis of multivariate data.  Journal of the American Statistical Association, 109:334–345.
  * Matteson et al., (2013) Matteson, D. S., James, N. A., Nicholson, W. B., and Segalini, L. C. (2013).  Locally stationary vector processes and adaptive multivariate modeling.  In Acoustics, Speech and Signal Processing (ICASSP), 2013 IEEE International Conference on, pages 8722–8726. IEEE.
  * Moore and Babij, (2017) Moore, G. and Babij, M. (2017).  Iceland’s great frost winter of 1917/1918 and its representation in reanalyses of the twentieth century.  Quarterly Journal of the Royal Meteorological Society, 143:508–520.
  * Olshen et al., (2004) Olshen, A. B., Venkatraman, E., Lucito, R., and Wigler, M. (2004).  Circular binary segmentation for the analysis of array-based dna copy number data.  Biostatistics, 5:557–572.
  * Rintoul et al., (2018) Rintoul, S., Chown, S., DeConto, R., England, M., Fricker, H., Masson-Delmotte, V., Naish, T., Siegert, M., and Xavier, J. (2018).  Choosing the future of antarctica.  Nature, 558:233–241.
  * Robbins et al., (2011) Robbins, M. W., Lund, R. B., Gallagher, C. M., and Lu, Q. (2011).  Changepoints in the north atlantic tropical cyclone record.  Journal of the American Statistical Association, 106:89–99.
  * Robinson et al., (2010) Robinson, L. F., Wager, T. D., and Lindquist, M. A. (2010).  Change point estimation in multi-subject fmri studies.  Neuroimage, 49:1581–1592.
  * Serreze and Meier, (2018) Serreze, M. C. and Meier, W. N. (2018).  The arctic’s sea ice cover: trends, variability, predictability, and comparisons to the antarctic.  Annals of the New York Academy of Sciences.
  * Spiriti et al., (2013) Spiriti, S., Eubank, R., Smith, P. W., and Young, D. (2013).  Knot selection for least-squares and penalized splines.  Journal of Statistical Computation and Simulation, 83:1020–1036.
  * Tibshirani et al., (2014) Tibshirani, R. J. et al. (2014).  Adaptive piecewise polynomial estimation via trend filtering.  The Annals of Statistics, 42:285–323.

Supplementary materials for “Detecting linear trend changes and point
anomalies in data sequences”

Hyeyoung Maeng and Piotr Fryzlewicz

London School of Economics and Political Science

##  Appendix B. Lemmas

###### Lemma 1

Let $\psi^{(j,k)}=\sum_{i=1}^{I^{(j,k)}}\phi_{i}^{(j,k)}g_{i}^{(j,k)}$ where
$\phi_{i}^{(j,k)}$ are constants and $g_{i}^{(j,k)}$ are vectors of equal
length with $\psi_{i}^{(j,k)}$ where
$I^{(j,k)}\in\\{3,4\\},j=1,\ldots,J,\;k=1,\ldots,K(j)$. If we define the set
$G=\\{g_{l}\\}$ where there is a unique correspondence between
$\bigg{\\{}{g_{i}^{(j,k)}}_{i=1,\ldots,I^{(j,k)},j=1,\ldots,J,\,k=1,\ldots,K(j)}\bigg{\\}}$
and $\\{g_{l}\\}$, we then have $P(A_{T})\geq 1-C_{2}T^{-1}$ where

$A_{T}=\bigg{\\{}\max_{g_{l}\in
G}\;|g_{l}^{T}\boldsymbol{\varepsilon}|\leq\lambda\bigg{\\}},$ (A.1)

$\lambda=C_{1}\sqrt{2\log T}$ and $C_{1}$ is a positive constant large enough
and $C_{2}$ is a positive constant.

Proof. We firstly show that for any fixed ${(j,k)}$, $g_{i}^{(j,k)}$ and
$\phi_{i}^{(j,k)}$ satisfy the conditions,
$\Big{(}g_{i}^{(j,k)}\Big{)}^{T}g_{i}^{(j,k)}=1$,
$\Big{(}g_{i}^{(j,k)}\Big{)}^{T}g_{i^{\prime}}^{(j,k)}=0$ and
$\sum_{i}\big{(}\phi_{i}^{(j,k)}\big{)}^{2}=1$, where
$\psi^{(j,k)}=\sum_{i=1}^{I^{(j,k)}}\phi_{i}^{(j,k)}g_{i}^{(j,k)}$. Depending
on the type of merge, $\psi^{(j,k)}$ fall into one of the followings,

$\displaystyle\begin{split}\text{Type 1:
}\psi^{(j,k)}_{p,q,r}&=\alpha_{1}e_{p}+\alpha_{2}e_{p+1}+\alpha_{3}e_{p+2},\\\
\text{Type 2:
}\psi^{(j,k)}_{p,q,r}&=\beta_{1}e_{p}+\beta_{2}(\underbrace{0,\ldots,0}_{p\times
1},\boldsymbol{\ell}_{1,p+1,r}^{T},\underbrace{0,\ldots,0}_{(T-r)\times
1})+\beta_{3}(\underbrace{0,\ldots,0}_{p\times
1},\boldsymbol{\ell}_{2,p+1,r}^{T},\underbrace{0,\ldots,0}_{(T-r)\times
1}),\\\ \psi^{(j,k)}_{p,q,r}&=\beta_{4}(\underbrace{0,\ldots,0}_{(p-1)\times
1},\boldsymbol{\ell}_{1,p,r-1}^{T},\underbrace{0,\ldots,0}_{(T-r+1)\times
1})+\beta_{5}(\underbrace{0,\ldots,0}_{(p-1)\times
1},\boldsymbol{\ell}_{2,p,r-1}^{T},\underbrace{0,\ldots,0}_{(T-r+1)\times
1})+\beta_{6}e_{r},\\\ \text{Type 3:
}\psi^{(j,k)}_{p,q,r}&=\gamma_{1}(\underbrace{0,\ldots,0}_{(p-1)\times
1},\boldsymbol{\ell}_{1,p,q}^{T},\underbrace{0,\ldots,0}_{(T-q)\times
1})+\gamma_{2}(\underbrace{0,\ldots,0}_{(p-1)\times
1},\boldsymbol{\ell}_{2,p,q}^{T},\underbrace{0,\ldots,0}_{(T-q)\times 1})\\\
&+\gamma_{3}(\underbrace{0,\ldots,0}_{q\times
1},\boldsymbol{\ell}_{1,q+1,r}^{T},\underbrace{0,\ldots,0}_{(T-r)\times
1})+\gamma_{4}(\underbrace{0,\ldots,0}_{q\times
1},\boldsymbol{\ell}_{2,q+1,r}^{T},\underbrace{0,\ldots,0}_{(T-r)\times
1}),\\\ \end{split}$ (A.2)

where $e_{i}$ is a vector of length $T$ with one $1$ at $i^{th}$ position and
zero for the others. As will be shown in Section E.,
$\boldsymbol{\ell}_{1,i,j}$ and $\boldsymbol{\ell}_{2,i,j}$ are an arbitrary
orthonormal basis of the subspace
$\\{(x_{1},x_{2},\ldots,x_{j-i+1})\;|\;x_{1}-x_{2}=x_{2}-x_{3}=\cdots=x_{j-i}-x_{j-i+1}\\}$
of $\mathbb{R}^{j-i+1}$.

In any case, we can obtain the representation
$\psi^{(j,k)}=\sum_{i=1}^{I^{(j,k)}}\phi_{i}^{(j,k)}g_{i}^{(j,k)}$ from (A.2)
if the constants $\phi_{i}^{(j,k)}$ correspond to $\\{\alpha_{i}\\}_{i=1}^{3}$
in Type 1, $\\{\beta_{i}\\}_{i=1}^{3}$ or $\\{\beta_{i}\\}_{i=4}^{6}$ in Type
2 and $\\{\gamma_{i}\\}_{i=1}^{4}$ in Type 3. From the orthonormality of the
basis ($\boldsymbol{\ell}_{1,\cdot,\cdot},\boldsymbol{\ell}_{2,\cdot,\cdot}$),
we see that the conditions, $\Big{(}g_{i}^{(j,k)}\Big{)}^{T}g_{i}^{(j,k)}=1$
and $\Big{(}g_{i}^{(j,k)}\Big{)}^{T}g_{i^{\prime}}^{(j,k)}=0$, are satisfied
for any $(i,i^{\prime},j,k)$ where $i\neq i^{\prime}$. In addition, as
$\psi^{(j,k)}$ keeps orthonormality, we can argue that $\phi_{i}^{(j,k)}$ is
bounded by the condition $\sum_{i}\big{(}\phi_{i}^{(j,k)}\big{)}^{2}=1$ for
any ${(i,j,k)}$ which implies
$\sum_{i=1}^{3}\alpha_{i}^{2}=\sum_{i=1}^{3}\beta_{i}^{2}=\sum_{i=4}^{6}\beta_{i}^{2}=\sum_{i=1}^{4}\gamma_{i}^{2}=1$
in (A.2).

If we predefine the pairs
($\boldsymbol{\ell}_{1,i,j},\boldsymbol{\ell}_{2,i,j}$) for all $i$ and $j$ by
choosing an orthonormal basis of the subspace
$\\{(x_{1},x_{2},\ldots,x_{j-i+1})\;|\;x_{1}-x_{2}=x_{2}-x_{3}=x_{j-i}-x_{j-i+1}\\}$
of $\mathbb{R}^{j-i+1}$, then there exist at most $2T^{2}$ vectors $g_{l}$ in
the set $G$. Now we are in position to show that $P(A_{T})\geq 1-C_{2}T^{-1}$.
Using a simple Bonferroni inequality, we have

$1-P(A_{T})\leq\sum_{G}P(|Z|>\lambda)\leq
2T^{2}\frac{\phi_{Z}(\lambda)}{\lambda}\leq\frac{C_{2}}{T}$ (A.3)

where $\phi_{Z}$ is the p.d.f. of a standard normal $Z$. This completes the
proof.

###### Lemma 2

Let $\mathcal{S}^{1}_{j}=\\{1\leq k\leq K(j):d^{(j,k)}$ is $d_{p,q,r}$ such
that $p<\eta_{i}-0.5<r$ for some $i=1,\ldots,N$ $\\}$, and
$\mathcal{S}^{0}_{j}=\\{1,\ldots,K(j)\\}\setminus\mathcal{S}^{1}_{j}$. On the
set $A_{T}$ which satisfies $P(A_{T})\rightarrow 1$ as $T\rightarrow\infty$,
we have

$\max_{\begin{subarray}{c}j=1,\ldots,J,\\\
k\in\mathcal{S}^{0}_{j}\end{subarray}}\big{|}d^{(j,k)}\big{|}\leq\lambda,$
(A.4)

where $\lambda$ is as in Theorem 1.

Proof. On the set $A_{T}$, the following holds for
$j=1,\ldots,J,k\in\mathcal{S}^{0}_{j}$,

$\displaystyle\big{|}d^{(j,k)}\big{|}=\;$
$\displaystyle\big{|}(\psi^{(j,k)})^{T}\boldsymbol{\varepsilon}\big{|}$
$\displaystyle=\;$
$\displaystyle\Big{|}\phi_{1}^{(j,k)}\Big{(}g_{1}^{(j,k)}\Big{)}^{T}\boldsymbol{\varepsilon}+\phi_{2}^{(j,k)}\Big{(}g_{2}^{(j,k)}\Big{)}^{T}\boldsymbol{\varepsilon}+\phi_{3}^{(j,k)}\Big{(}g_{3}^{(j,k)}\Big{)}^{T}\boldsymbol{\varepsilon}+\phi_{4}^{(j,k)}\Big{(}g_{4}^{(j,k)}\Big{)}^{T}\boldsymbol{\varepsilon}\Big{|}$
$\displaystyle\leq\;$
$\displaystyle\max_{j,\;k}\Big{(}\big{|}\phi_{1}^{(j,k)}\big{|}+\big{|}\phi_{2}^{(j,k)}\big{|}+\big{|}\phi_{3}^{(j,k)}\big{|}+\big{|}\phi_{4}^{(j,k)}\big{|}\Big{)}\cdot\bigg{(}\max_{l:\;g_{l}\in
G}\;\big{|}g_{l}^{T}\boldsymbol{\varepsilon}\big{|}\bigg{)}$

where $\boldsymbol{\varepsilon}=(\varepsilon_{1},\ldots,\varepsilon_{T})^{T}$.
The condition, $\sum_{i}\Big{(}\phi_{i}^{(j,k)}\Big{)}^{2}=1$ for any fixed
${(i,j,k)}$, given in the proof of Lemma 1 implies
$\max_{i}\big{|}\phi_{i}^{(j,k)}\big{|}\leq 1$ for any ${(j,k)}$, thus we have
(A.4) when the constant $C_{1}$ for $\lambda$ in (A.4) is larger than or equal
to $4$ times $C_{1}$ used in (A.1).

##  Appendix C. Additional simulation results

In addition to the simulations in Section 4, here we present the results for
three different distributions of the noise $\varepsilon_{t}$, (a)
$\varepsilon_{t}$ follows a stationary Gaussian AR(1) process of $\phi=0.3$,
with zero-mean and unit-variance, (b) $\varepsilon_{t}\sim$ i.i.d. scaled
$t_{5}$ distribution with unit-variance and (c) $\varepsilon_{t}\sim$ i.i.d.
Laplace(0, $1/\sqrt{2}$). We use $C=1.8$ as a default thresholding constant
for TS. Among other competitors, only ID provides the option for heavy-tailed
noise in their R package IDetect and other methods are set to their default
settings.

Table 4: Distribution of $\hat{N}-N$ for models (M1)-(M4) and all methods with the noise term $\epsilon_{t}$ being $AR(1)$ process of $\phi=0.3$ over 100 simulation runs. Also the average MSE (Mean Squared Error) of the estimated signal $\hat{f}_{t}$, the average Hausdorff distance $d_{H}$ and the average computational time in seconds using an Intel Core i5 2.9 GHz CPU with 8 GB of RAM, all over 100 simulations. Bold: methods within 10% of the highest empirical frequency of $\hat{N}-N=0$ or within 10% of the lowest empirical average $d_{H}(\times 10^{2})$. |  |  |  |  | $\hat{N}-N$ |  |  |  |  |  |
---|---|---|---|---|---|---|---|---|---|---|---
Model | Method |  $\leq$-3 | -2 | -1 | 0 | 1 | 2 |  $\geq$3 | MSE | $d_{H}(\times 10^{2})$ | time
(M1) | TS | 0 | 0 | 4 | 93 | 3 | 0 | 0 | 0.081 | 3.58 | 1.27
NOT | 0 | 0 | 0 | 91 | 9 | 0 | 0 | 0.072 | 2.95 | 0.26
ID | 0 | 0 | 0 | 82 | 14 | 4 | 0 | 0.067 | 2.65 | 0.38
TF | 0 | 0 | 0 | 0 | 0 | 0 | 100 | 0.532 | 5.00 | 36.61
CPOP | 0 | 0 | 0 | 7 | 15 | 6 | 72 | 0.080 | 3.69 | 4.95
BUP | 0 | 0 | 8 | 86 | 6 | 0 | 0 | 0.077 | 3.66 | 2.75
(M2) | TS | 1 | 6 | 23 | 69 | 1 | 0 | 0 | 0.195 | 2.44 | 1.12
NOT | 0 | 0 | 8 | 83 | 6 | 2 | 1 | 0.182 | 2.11 | 0.31
ID | 0 | 0 | 0 | 69 | 24 | 5 | 2 | 0.155 | 1.75 | 0.40
TF | 0 | 0 | 0 | 0 | 0 | 0 | 100 | 0.600 | 2.38 | 32.03
CPOP | 0 | 0 | 0 | 1 | 6 | 8 | 85 | 0.163 | 1.98 | 1.50
BUP | 100 | 0 | 0 | 0 | 0 | 0 | 0 | 0.717 | 4.63 | 2.39
(M3) | TS | 0 | 1 | 5 | 88 | 6 | 0 | 0 | 0.052 | 4.16 | 1.56
NOT | 0 | 0 | 0 | 89 | 7 | 4 | 0 | 0.042 | 3.40 | 0.31
ID | 0 | 0 | 3 | 77 | 16 | 3 | 1 | 0.064 | 3.12 | 0.50
TF | 0 | 0 | 0 | 0 | 0 | 0 | 100 | 0.259 | 6.24 | 44.94
CPOP | 0 | 0 | 0 | 1 | 4 | 10 | 85 | 0.068 | 4.67 | 9.57
BUP | 0 | 0 | 0 | 0 | 3 | 18 | 79 | 0.056 | 5.57 | 3.56
(M4) | TS | 0 | 6 | 23 | 53 | 18 | 0 | 0 | 0.058 | 2.41 | 1.53
NOT | 0 | 93 | 6 | 1 | 0 | 0 | 0 | 0.086 | 2.91 | 0.31
ID | 2 | 6 | 30 | 49 | 10 | 2 | 1 | 0.165 | 2.99 | 0.48
TF | 0 | 0 | 0 | 0 | 0 | 0 | 100 | 0.218 | 6.22 | 45.99
CPOP | 0 | 0 | 0 | 1 | 3 | 9 | 87 | 0.066 | 4.02 | 5.40
BUP | 0 | 0 | 0 | 11 | 35 | 37 | 17 | 0.109 | 5.64 | 3.77
Table 5: Distribution of $\hat{N}-N$ for models (M5)-(M8) and all methods with the noise term $\epsilon_{t}$ being $AR(1)$ process of $\phi=0.3$ over 100 simulation runs. Also the average MSE (Mean Squared Error) of the estimated signal $\hat{f}_{t}$, the average Hausdorff distance $d_{H}$ and the average computational time in seconds using an Intel Core i5 2.9 GHz CPU with 8 GB of RAM, all over 100 simulations. Bold: methods within 10% of the highest empirical frequency of $\hat{N}-N=0$ or within 10% of the lowest empirical average $d_{H}(\times 10^{2})$. |  |  |  |  | $\hat{N}-N$ |  |  |  |  |  |
---|---|---|---|---|---|---|---|---|---|---|---
Model | Method |  $\leq$-3 | -2 | -1 | 0 | 1 | 2 |  $\geq$3 | MSE | $d_{H}(\times 10^{2})$ | time
(M5) | TS | 0 | 0 | 19 | 54 | 21 | 6 | 0 | 0.055 | 1.87 | 1.54
NOT | 0 | 0 | 91 | 6 | 3 | 0 | 0 | 0.060 | 1.94 | 0.28
ID | 0 | 0 | 9 | 23 | 23 | 18 | 27 | 0.402 | 9.47 | 0.46
TF | 0 | 0 | 0 | 0 | 0 | 0 | 100 | 0.182 | 6.21 | 42.65
CPOP | 0 | 0 | 0 | 0 | 2 | 4 | 94 | 0.068 | 3.70 | 4.08
BUP | 0 | 0 | 0 | 15 | 37 | 32 | 16 | 0.112 | 5.25 | 3.53
(M6) | TS | 0 | 0 | 0 | 98 | 2 | 0 | 0 | 0.018 | 0.06 | 1.70
NOT | 68 | 9 | 10 | 4 | 1 | 3 | 5 | 0.257 | 21.63 | 0.25
ID | 20 | 10 | 0 | 0 | 11 | 0 | 59 | 0.164 | 12.83 | 0.63
TF | 0 | 0 | 0 | 0 | 0 | 0 | 100 | 0.332 | 11.04 | 47.43
CPOP | 0 | 0 | 0 | 5 | 11 | 17 | 67 | 0.056 | 4.86 | 5.31
BUP | 0 | 0 | 0 | 0 | 0 | 0 | 100 | 0.170 | 10.18 | 3.95
(M7) | TS | 11 | 38 | 31 | 15 | 3 | 2 | 0 | 0.217 | 11.52 | 0.68
NOT | 5 | 12 | 19 | 24 | 22 | 7 | 11 | 0.158 | 7.69 | 0.24
ID | 32 | 1 | 18 | 26 | 14 | 5 | 4 | 0.511 | 17.54 | 0.03
TF | 3 | 0 | 0 | 0 | 0 | 0 | 97 | 0.623 | 7.01 | 23.25
CPOP | 0 | 0 | 0 | 0 | 0 | 1 | 99 | 0.162 | 5.27 | 0.85
BUP | 54 | 43 | 3 | 0 | 0 | 0 | 0 | 0.283 | 11.92 | 1.55
(M8) | TS | 0 | 0 | 0 | 100 | 0 | 0 | 0 | 0.003 | 0.00 | 1.09
NOT | 0 | 0 | 0 | 93 | 3 | 3 | 1 | 0.005 | 2.02 | 0.19
ID | 0 | 0 | 0 | 100 | 0 | 0 | 0 | 0.003 | 0.00 | 0.51
TF | 0 | 0 | 0 | 0 | 0 | 0 | 100 | 0.551 | 49.94 | 35.81
CPOP | 0 | 0 | 0 | 30 | 10 | 3 | 57 | 0.035 | 19.71 | 7.55
BUP | 0 | 0 | 0 | 0 | 0 | 0 | 100 | 0.025 | 46.73 | 2.72
Table 6: Distribution of $\hat{N}-N$ for models (M1)-(M4) and all methods with the noise term $\epsilon_{t}$ being i.i.d. $t_{5}$ over 100 simulation runs. Also the average MSE (Mean Squared Error) of the estimated signal $\hat{f}_{t}$, the average Hausdorff distance $d_{H}$ and the average computational time in seconds using an Intel Core i5 2.9 GHz CPU with 8 GB of RAM, all over 100 simulations. Bold: methods within 10% of the highest empirical frequency of $\hat{N}-N=0$ or within 10% of the lowest empirical average $d_{H}(\times 10^{2})$. |  |  |  |  | $\hat{N}-N$ |  |  |  |  |  |
---|---|---|---|---|---|---|---|---|---|---|---
Model | Method |  $\leq$-3 | -2 | -1 | 0 | 1 | 2 |  $\geq$3 | MSE | $d_{H}(\times 10^{2})$ | time
(M1) | TS | 0 | 0 | 3 | 55 | 32 | 5 | 5 | 0.083 | 3.41 | 1.09
NOT | 0 | 0 | 0 | 98 | 2 | 0 | 0 | 0.037 | 2.21 | 0.26
ID | 0 | 0 | 0 | 85 | 10 | 4 | 1 | 0.036 | 1.87 | 0.29
TF | 0 | 0 | 0 | 0 | 0 | 0 | 100 | 0.017 | 4.36 | 36.49
CPOP | 0 | 0 | 0 | 21 | 20 | 20 | 39 | 0.064 | 2.28 | 5.69
BUP | 0 | 4 | 13 | 78 | 5 | 0 | 0 | 0.071 | 3.84 | 2.62
(M2) | TS | 0 | 3 | 11 | 70 | 10 | 5 | 1 | 0.164 | 2.18 | 1.03
NOT | 0 | 0 | 3 | 85 | 11 | 0 | 1 | 0.098 | 1.69 | 0.29
ID | 0 | 0 | 0 | 77 | 21 | 2 | 0 | 0.102 | 1.36 | 0.38
TF | 0 | 0 | 0 | 0 | 0 | 0 | 100 | 0.067 | 2.29 | 31.41
CPOP | 0 | 0 | 0 | 14 | 23 | 25 | 38 | 0.119 | 1.54 | 1.66
BUP | 100 | 0 | 0 | 0 | 0 | 0 | 0 | 0.752 | 4.69 | 2.18
(M3) | TS | 0 | 1 | 11 | 41 | 25 | 10 | 12 | 0.073 | 4.90 | 1.44
NOT | 0 | 0 | 0 | 96 | 3 | 1 | 0 | 0.021 | 2.54 | 0.31
ID | 0 | 0 | 1 | 73 | 19 | 3 | 4 | 0.053 | 2.72 | 0.44
TF | 0 | 0 | 0 | 0 | 0 | 0 | 100 | 0.024 | 5.92 | 46.35
CPOP | 0 | 0 | 0 | 9 | 10 | 11 | 70 | 0.065 | 3.57 | 11.71
BUP | 0 | 0 | 0 | 1 | 21 | 40 | 38 | 0.043 | 5.44 | 3.52
(M4) | TS | 0 | 3 | 14 | 34 | 23 | 16 | 10 | 0.075 | 3.10 | 1.46
NOT | 0 | 97 | 3 | 0 | 0 | 0 | 0 | 0.066 | 2.45 | 0.28
ID | 1 | 12 | 22 | 48 | 10 | 3 | 4 | 0.159 | 2.42 | 0.42
TF | 0 | 0 | 0 | 0 | 0 | 0 | 100 | 0.081 | 6.06 | 45.74
CPOP | 0 | 0 | 0 | 4 | 4 | 15 | 77 | 0.062 | 3.37 | 5.15
BUP | 0 | 2 | 7 | 28 | 47 | 12 | 4 | 0.095 | 5.30 | 3.56
Table 7: Distribution of $\hat{N}-N$ for models (M5)-(M8) and all methods with the noise term $\epsilon_{t}$ being i.i.d. $t_{5}$ over 100 simulation runs. Also the average MSE (Mean Squared Error) of the estimated signal $\hat{f}_{t}$, the average Hausdorff distance $d_{H}$ and the average computational time in seconds using an Intel Core i5 2.9 GHz CPU with 8 GB of RAM, all over 100 simulations. Bold: methods within 10% of the highest empirical frequency of $\hat{N}-N=0$ or within 10% of the lowest empirical average $d_{H}(\times 10^{2})$. |  |  |  |  | $\hat{N}-N$ |  |  |  |  |  |
---|---|---|---|---|---|---|---|---|---|---|---
Model | Method |  $\leq$-3 | -2 | -1 | 0 | 1 | 2 |  $\geq$3 | MSE | $d_{H}(\times 10^{2})$ | time
(M5) | TS | 0 | 0 | 9 | 40 | 19 | 21 | 11 | 0.069 | 2.72 | 1.49
NOT | 0 | 0 | 95 | 4 | 1 | 0 | 0 | 0.042 | 1.29 | 0.26
ID | 0 | 0 | 1 | 16 | 24 | 23 | 36 | 0.372 | 9.86 | 0.43
TF | 0 | 0 | 0 | 0 | 0 | 0 | 100 | 0.118 | 6.15 | 43.23
CPOP | 0 | 0 | 0 | 3 | 8 | 12 | 77 | 0.060 | 2.97 | 3.55
BUP | 0 | 0 | 10 | 40 | 43 | 6 | 1 | 0.083 | 4.76 | 3.51
(M6) | TS | 0 | 0 | 0 | 46 | 2 | 39 | 13 | 0.035 | 3.16 | 1.66
NOT | 54 | 21 | 4 | 8 | 5 | 1 | 7 | 0.244 | 17.30 | 0.23
ID | 8 | 8 | 0 | 0 | 6 | 0 | 78 | 0.125 | 6.99 | 0.62
TF | 0 | 0 | 0 | 0 | 0 | 0 | 100 | 0.138 | 10.99 | 48.53
CPOP | 0 | 0 | 0 | 9 | 11 | 17 | 63 | 0.059 | 4.68 | 3.48
BUP | 0 | 0 | 0 | 0 | 0 | 0 | 100 | 0.145 | 10.27 | 3.92
(M7) | TS | 14 | 28 | 32 | 14 | 8 | 4 | 0 | 0.204 | 11.21 | 0.65
NOT | 0 | 6 | 16 | 30 | 36 | 11 | 1 | 0.079 | 5.12 | 0.22
ID | 14 | 8 | 12 | 17 | 24 | 13 | 12 | 0.421 | 16.22 | 0.04
TF | 0 | 0 | 0 | 0 | 0 | 0 | 100 | 0.098 | 6.08 | 23.86
CPOP | 0 | 0 | 0 | 0 | 4 | 5 | 91 | 0.102 | 3.01 | 0.81
BUP | 69 | 28 | 3 | 0 | 0 | 0 | 0 | 0.266 | 12.12 | 1.47
(M8) | TS | 0 | 0 | 0 | 49 | 0 | 43 | 8 | 0.030 | 14.86 | 1.01
NOT | 0 | 0 | 0 | 100 | 0 | 0 | 0 | 0.001 | 0.00 | 0.17
ID | 0 | 0 | 0 | 99 | 1 | 0 | 0 | 0.001 | 0.00 | 0.45
TF | 0 | 0 | 0 | 65 | 12 | 9 | 14 | 0.003 | 14.63 | 36.03
CPOP | 0 | 0 | 0 | 35 | 0 | 34 | 31 | 0.042 | 20.53 | 3.91
BUP | 0 | 0 | 0 | 0 | 0 | 0 | 100 | 0.014 | 46.80 | 2.62
Table 8: Distribution of $\hat{N}-N$ for models (M1)-(M4) and all methods with the noise term $\epsilon_{t}$ being i.i.d. Laplace over 100 simulation runs. Also the average MSE (Mean Squared Error) of the estimated signal $\hat{f}_{t}$, the average Hausdorff distance $d_{H}$ and the average computational time in seconds using an Intel Core i5 2.9 GHz CPU with 8 GB of RAM, all over 100 simulations. Bold: methods within 10% of the highest empirical frequency of $\hat{N}-N=0$ or within 10% of the lowest empirical average $d_{H}(\times 10^{2})$. |  |  |  |  | $\hat{N}-N$ |  |  |  |  |  |
---|---|---|---|---|---|---|---|---|---|---|---
Model | Method |  $\leq$-3 | -2 | -1 | 0 | 1 | 2 |  $\geq$3 | MSE | $d_{H}(\times 10^{2})$ | time
(M1) | TS | 5 | 17 | 36 | 33 | 8 | 1 | 0 | 0.167 | 5.68 | 1.10
NOT | 0 | 1 | 10 | 86 | 3 | 0 | 0 | 0.086 | 3.56 | 0.28
ID | 0 | 0 | 0 | 89 | 9 | 2 | 0 | 0.065 | 2.33 | 0.32
TF | 0 | 0 | 0 | 0 | 0 | 0 | 100 | 0.034 | 4.44 | 36.66
CPOP | 0 | 0 | 0 | 28 | 30 | 16 | 26 | 0.091 | 2.21 | 7.63
BUP | 0 | 0 | 0 | 0 | 0 | 4 | 96 | 0.080 | 4.64 | 2.64
(M2) | TS | 69 | 16 | 9 | 4 | 2 | 0 | 0 | 0.466 | 5.81 | 1.03
NOT | 5 | 8 | 36 | 45 | 4 | 2 | 0 | 0.240 | 2.80 | 0.33
ID | 0 | 1 | 0 | 84 | 10 | 3 | 2 | 0.158 | 1.87 | 0.39
TF | 0 | 0 | 0 | 0 | 0 | 0 | 100 | 0.102 | 2.30 | 31.92
CPOP | 0 | 0 | 0 | 26 | 16 | 20 | 38 | 0.182 | 1.71 | 2.20
BUP | 64 | 21 | 14 | 1 | 0 | 0 | 0 | 0.375 | 3.71 | 2.23
(M3) | TS | 2 | 16 | 30 | 33 | 17 | 2 | 0 | 0.127 | 7.75 | 1.45
NOT | 0 | 0 | 11 | 83 | 5 | 0 | 1 | 0.049 | 3.74 | 0.33
ID | 0 | 0 | 4 | 83 | 11 | 2 | 0 | 0.061 | 3.29 | 0.48
TF | 0 | 0 | 0 | 0 | 0 | 0 | 100 | 0.038 | 5.75 | 46.08
CPOP | 0 | 0 | 0 | 17 | 13 | 16 | 54 | 0.100 | 3.90 | 16.09
BUP | 0 | 0 | 0 | 0 | 0 | 0 | 100 | 0.072 | 5.98 | 3.48
(M4) | TS | 11 | 42 | 32 | 9 | 5 | 1 | 0 | 0.154 | 4.94 | 1.46
NOT | 1 | 89 | 8 | 2 | 0 | 0 | 0 | 0.092 | 3.12 | 0.31
ID | 11 | 8 | 44 | 27 | 7 | 1 | 2 | 0.188 | 3.77 | 0.49
TF | 0 | 0 | 0 | 0 | 0 | 0 | 100 | 0.099 | 5.91 | 45.56
CPOP | 0 | 2 | 5 | 13 | 13 | 20 | 47 | 0.118 | 3.31 | 8.35
BUP | 0 | 0 | 0 | 0 | 0 | 0 | 100 | 0.132 | 6.00 | 3.51
Table 9: Distribution of $\hat{N}-N$ for models (M5)-(M8) and all methods with the noise term $\epsilon_{t}$ being i.i.d. Laplace over 100 simulation runs. Also the average MSE (Mean Squared Error) of the estimated signal $\hat{f}_{t}$, the average Hausdorff distance $d_{H}$ and the average computational time in seconds using an Intel Core i5 2.9 GHz CPU with 8 GB of RAM, all over 100 simulations. Bold: methods within 10% of the highest empirical frequency of $\hat{N}-N=0$ or within 10% of the lowest empirical average $d_{H}(\times 10^{2})$. |  |  |  |  | $\hat{N}-N$ |  |  |  |  |  |
---|---|---|---|---|---|---|---|---|---|---|---
Model | Method |  $\leq$-3 | -2 | -1 | 0 | 1 | 2 |  $\geq$3 | MSE | $d_{H}(\times 10^{2})$ | time
(M5) | TS | 0 | 6 | 43 | 22 | 23 | 5 | 1 | 0.141 | 3.49 | 1.49
NOT | 0 | 0 | 86 | 13 | 1 | 0 | 0 | 0.074 | 2.12 | 0.30
ID | 0 | 4 | 12 | 38 | 20 | 17 | 9 | 0.448 | 9.50 | 0.46
TF | 0 | 0 | 0 | 0 | 0 | 0 | 100 | 0.168 | 6.14 | 44.16
CPOP | 0 | 0 | 1 | 4 | 12 | 16 | 67 | 0.113 | 3.10 | 5.21
BUP | 0 | 0 | 0 | 0 | 0 | 0 | 100 | 0.148 | 5.91 | 3.55
(M6) | TS | 4 | 13 | 0 | 57 | 3 | 22 | 1 | 0.077 | 6.14 | 1.64
NOT | 76 | 11 | 11 | 1 | 1 | 0 | 0 | 0.260 | 27.86 | 0.27
ID | 25 | 13 | 1 | 0 | 23 | 0 | 38 | 0.181 | 16.77 | 0.64
TF | 2 | 0 | 0 | 0 | 0 | 0 | 98 | 0.230 | 11.40 | 48.91
CPOP | 0 | 0 | 0 | 6 | 13 | 16 | 65 | 0.092 | 3.98 | 3.63
BUP | 0 | 0 | 0 | 0 | 0 | 0 | 100 | 0.211 | 10.59 | 3.97
(M7) | TS | 52 | 35 | 11 | 1 | 1 | 0 | 0 | 0.353 | 12.30 | 0.65
NOT | 42 | 41 | 8 | 5 | 2 | 2 | 0 | 0.232 | 11.71 | 0.25
ID | 45 | 5 | 21 | 12 | 11 | 4 | 2 | 0.607 | 24.48 | 0.04
TF | 0 | 0 | 0 | 0 | 0 | 0 | 100 | 0.147 | 6.04 | 22.68
CPOP | 0 | 0 | 3 | 9 | 12 | 15 | 61 | 0.260 | 9.59 | 1.39
BUP | 0 | 0 | 4 | 35 | 42 | 18 | 1 | 0.232 | 5.64 | 1.54
(M8) | TS | 0 | 0 | 0 | 67 | 0 | 31 | 2 | 0.026 | 8.72 | 1.01
NOT | 0 | 0 | 0 | 99 | 1 | 0 | 0 | 0.003 | 0.01 | 0.19
ID | 0 | 0 | 0 | 99 | 1 | 0 | 0 | 0.003 | 0.01 | 0.49
TF | 0 | 0 | 0 | 61 | 10 | 6 | 23 | 0.006 | 16.88 | 36.09
CPOP | 0 | 0 | 0 | 47 | 2 | 25 | 26 | 0.047 | 15.32 | 4.04
BUP | 0 | 0 | 0 | 0 | 0 | 0 | 100 | 0.056 | 48.08 | 2.66

##  Appendix D. Additional data application results

(a) NOT and ID

(b) TF and CPOP

Figure 5: Change-point analysis for the monthly average sea ice extent of the
Arctic in February from 1979 to 2018 in Section 5.2. (a) the data series (grey
dots) and estimated signal with change-points returned by NOT ( ) and ID ( ),
(b) estimated signal with change-points returned by TF ( ) and CPOP ( ).

(a) NOT and ID

(b) TF and CPOP

Figure 6: Change-point analysis for the monthly average sea ice extent of the
Arctic in September from 1979 to 2018 in Section 5.2. (a) the data series
(grey dots) and estimated signal with change-points returned by NOT ( ) and ID
( ), (b) estimated signal with change-points returned by TF ( ) and CPOP ( ).

(a) NOT and ID

(b) TF and CPOP

Figure 7: Change-point analysis for the monthly average sea ice extent of the
Antarctic in February from 1979 to 2018 in Section 5.2. (a) the data series
(grey dots) and estimated signal with change-points returned by NOT ( ) and ID
( ), (b) estimated signal with change-points returned by TF ( ) and CPOP ( ).

(a) NOT and ID

(b) TF and CPOP

Figure 8: Change-point analysis for the monthly average sea ice extent of the
Antarctic in September from 1979 to 2018 in Section 5.2. (a) the data series
(grey dots) and estimated signal with change-points returned by NOT ( ) and ID
( ), (b) estimated signal with change-points returned by TF ( ) and CPOP ( ).

##  Appendix E. Geometric interpretation of the TGUW transformation

In this section, we explore the shape of the unbalanced wavelet basis and
offer the geometric interpretation of the TGUW transformation through an
illustrative example. We recall that the TGUW transformation can be summarised
as an orthonormal transformation of initial data $X$ using orthonormal basis
for $\mathbb{R}^{T}$ as follows,

$\begin{pmatrix}s^{1}_{1,T}\\\ s^{2}_{1,T}\\\ \begin{pmatrix}\\\
{d^{(j,k)}}_{j=1,\ldots,J,k=1,\ldots,K(j)}\\\ \\\
\end{pmatrix}\end{pmatrix}_{T\times 1}=\quad\begin{pmatrix}\psi^{(0,1)}\\\
\psi^{(0,2)}\\\ \begin{pmatrix}\\\
{\psi^{(j,k)}}_{j=1,\ldots,J,k=1,\ldots,K(j)}\\\ \\\
\end{pmatrix}\end{pmatrix}_{T\times T}\begin{pmatrix}X_{1}\\\ X_{2}\\\
\vdots\\\ X_{T}\\\ \end{pmatrix}_{T\times 1}=\quad\Psi_{T\times
T}\begin{pmatrix}X_{1}\\\ X_{2}\\\ \vdots\\\ X_{T}\\\ \end{pmatrix},$ (A.5)

where $\Psi$ is an orthogonal matrix, $\\{\psi^{(j,k)}\\}$ is the set of
vectors such that $d^{(j,k)}=\langle X,\psi^{(j,k)}\rangle$ and
($\psi^{(0,1)},\psi^{(0,2)}$) are vectors satisfying $s^{1}_{1,T}=\langle
X,\psi^{(0,1)}\rangle$ and $s^{2}_{1,T}=\langle X,\psi^{(0,2)}\rangle$,
respectively. We note that there is a unique correspondence between
$\psi_{p,q,r}^{(j,k)}$ and $\psi^{(j,k)}$ where $j=1,\ldots,J$,
$k=1,\ldots,K(j)$ as is also the case for $d_{p,q,r}^{(j,k)}$ and $d^{(j,k)}$.

###  E.1 Shape of the unbalanced wavelet basis

The orthogonal transformation matrix $\Psi$ in (A.5) contains $T$ orthonormal
basis in its rows which can be categorised into two: 1)
$\Big{\\{}{\psi^{(j,k)}}_{j=1,\ldots,J,k=1,\ldots,K(j)}\Big{\\}}$
corresponding to detail coefficients
$\Big{\\{}{d^{(j,k)}}_{j=1,\ldots,J,k=1,\ldots,K(j)}\Big{\\}}$ and 2)
$\psi^{(0,1)}$ and $\psi^{(0,2)}$ corresponding to two smooth coefficients,
$s^{1}_{1,T}=\langle X,\psi^{(0,1)}\rangle$ and $s^{2}_{1,T}=\langle
X,\psi^{(0,2)}\rangle$, respectively.

We firstly investigate the shape of wavelet basis for detail coefficients as
it plays an important role in the proof of theorems in Section B. The basis
vector $\psi_{p,q,r}^{(j,k)}$ has length $T$ but has non-zero values only from
its $p^{\text{th}}$ element to $r^{\text{th}}$ element where $r-p\geq 2$ and
the shape of non-zero region of $\psi_{p,q,r}^{(j,k)}$ depends on the type of
merging which produces the corresponding detail coefficient. The general form
in (A.2) implies that $\psi_{p,q,r}^{(j,k)}$ always has a shape of linear
trend in previously merged region and this linearity is preserved as long as
the later merges are performed under the “two together” rule. This is due to
that the orthonormal transforms continue in a way of extending the dimension
of space in which an orthonormal basis lives We now give a geometric
interpretation of the TGUW transformation through an example to show how the
linear trend in wavelet basis is preserved through the consecutive orthonormal
transforms.

###  E.2 Example

We use a simple example to explain how the basis $\psi_{p,q,r}^{(j,k)}$ in
(A.2) keep its linearity in subregions which are previously merged. Suppose we
have the initial data sequence $\boldsymbol{s}^{0}=(X_{1},\ldots,X_{5})$ and
the initial weight vectors for constancy and linearity are
$\boldsymbol{w}^{c,0}=(1,1,1,1,1)^{\text{T}}$ and
$\boldsymbol{w}^{l,0}=(1,2,3,4,5)^{\text{T}}$, respectively. As we have 5
initial observations, there would be 3 orthonormal transformations and the
most important task for each transform is finding the appropriate orthonormal
matrix as follows,

$\Lambda=\begin{pmatrix}&\ell_{1,1}&\ell_{1,2}&\ell_{1,3}&\\\
&\ell_{2,1}&\ell_{2,2}&\ell_{2,3}&\\\
&a&b&c&\end{pmatrix}=\begin{pmatrix}&\boldsymbol{\ell}_{1}^{\text{T}}&\\\
&\boldsymbol{\ell}_{2}^{\text{T}}&\\\
&\boldsymbol{h}^{\text{T}}\end{pmatrix}.$ (A.6)

##### First merge

Assume that $(X_{3},X_{4},X_{5})$ is chosen as the first triplet to be merged.
To find the values of the matrix $\Lambda$, we firstly seek the detail filer,
$\boldsymbol{h}$, which satisfies the conditions (1)
$\boldsymbol{h}^{\text{T}}\boldsymbol{w}^{c,0}_{3,4,5}=0$, (2)
$\boldsymbol{h}^{\text{T}}\boldsymbol{w}^{l,0}_{3,4,5}=0$ and (3)
$\boldsymbol{h}^{\text{T}}\boldsymbol{h}=1$, where
$\boldsymbol{w}^{\cdot}_{p,q,r}$ is $(p,q,r)^{th}$ element of
$\boldsymbol{w}^{\cdot}$. Thus, $\boldsymbol{h}$ is obtained as a normal
vector of the plane $\\{(x,y,z)\;|\;x-2y+z=0\\}$ as
$\boldsymbol{w}^{c,0}_{3,4,5}$ and $\boldsymbol{w}^{l,0}_{3,4,5}$ are vectors
of that plane. From the detail filter $\boldsymbol{h}$, two low filter vectors
($\boldsymbol{\ell}_{1}$ and $\boldsymbol{\ell}_{2}$) are obtained under the
conditions, (1) $\boldsymbol{\ell}_{1}^{\text{T}}\boldsymbol{h}=0$, (2)
$\boldsymbol{\ell}_{2}^{\text{T}}\boldsymbol{h}=0$, (3)
$\boldsymbol{\ell}_{1}^{\text{T}}\boldsymbol{\ell}_{2}=0$ and (4)
$\boldsymbol{\ell}_{1}^{\text{T}}\boldsymbol{\ell}_{1}=\boldsymbol{\ell}_{2}^{\text{T}}\boldsymbol{\ell}_{2}=1$.
Therefore, $\boldsymbol{\ell}_{1}$ and $\boldsymbol{\ell}_{2}$ are obtained as
an orthonormal basis of the plane $\\{(x,y,z)\;|\;x-2y+z=0\\}$ and this
guarantees the linear trend of $\boldsymbol{\ell}_{1}$ and
$\boldsymbol{\ell}_{2}$. In summary,
$\boldsymbol{w}^{c,0}_{3,4,5},\boldsymbol{w}^{l,0}_{3,4,5},\boldsymbol{\ell}_{1},\boldsymbol{\ell}_{2}$
are vectors of the plane $\\{(x,y,z)\;|\;x-2y+z=0\\}$ and $\boldsymbol{h}$ is
a normal vector which is orthogonal to that plane. Now, the orthonormal
transform updates the data sequence and weight vectors as follows,

$\displaystyle\begin{split}\boldsymbol{s}^{0}=(X_{1},\ldots,X_{5})&\quad\rightarrow\quad\boldsymbol{s}=(X_{1},X_{2},s^{1}_{3,4,5},s^{2}_{3,4,5},d_{3,4,5}),\\\
\boldsymbol{w}^{c,0}=(1,1,1,1,1)^{\text{T}}&\quad\rightarrow\quad{\boldsymbol{w}^{c}}=(1,1,e_{c_{1}},e_{c_{2}},0)^{\text{T}},\\\
\boldsymbol{w}^{l,0}=(1,2,3,4,5)^{\text{T}}&\quad\rightarrow\quad{\boldsymbol{w}^{l}}=(1,2,e_{l_{1}},e_{l_{2}},0)^{\text{T}},\\\
\end{split}$ (A.7)

where constants $(e_{c_{1}},e_{c_{2}})$ and $(e_{l_{1}},e_{l_{2}})$ are
obtained from
$\Lambda\boldsymbol{w}^{c,0}=_{3,4,5}=(e_{c_{1}},e_{c_{2}},0)^{\text{T}}$ and
$\Lambda\boldsymbol{w}^{l,0}_{3,4,5}=(e_{l_{1}},e_{l_{2}},0)^{\text{T}}$,
respectively. As $\boldsymbol{\ell}_{1}$ and $\boldsymbol{\ell}_{2}$ are an
orthonormal basis of the plane $\\{(x,y,z)\;|\;x-2y+z=0\\}$,
$e_{c_{1}},e_{c_{2}}$ and $e_{l_{1}},e_{l_{2}}$ are unique constants which
represent $\boldsymbol{w}^{c,0}_{3,4,5}$ and $\boldsymbol{w}^{l,0}_{3,4,5}$ as
a linear span of basis $\boldsymbol{\ell}_{1}$ and $\boldsymbol{\ell}_{2}$ as
follows,

$\displaystyle\begin{split}\boldsymbol{w}^{c,0}_{3,4,5}=e_{c_{1}}\boldsymbol{\ell}_{1}+e_{c_{2}}\boldsymbol{\ell}_{2},\quad\boldsymbol{w}^{l,0}_{3,4,5}=e_{l_{1}}\boldsymbol{\ell}_{1}+e_{l_{2}}\boldsymbol{\ell}_{2}.\\\
\end{split}$ (A.8)

To show how the linear trend in the subregion of $\psi$ is preserved through
recursive orthonormal transforms, we now describe the process of obtaining the
matrix $\Psi_{T\times T}$ in (A.5) (an orthonormal basis of $\mathbb{R}^{5}$
in this example) from the initial input $\Psi^{0}=\mathbf{I}$ where
$\mathbf{I}$ is $5\times 5$ identity matrix. When $(p,q,r)^{th}$ values in
$\boldsymbol{s}$ are selected to be merged, we extract the corresponding three
columns of $\Psi^{\text{T}}$ and update them by matrix multiplication with the
orthogonal matrix $\Lambda$ defined for that merge. For example, after the
first orthonormal transform for $(X_{3},X_{4},X_{5})$ is performed as in
(A.7), we have the following updated matrix,

$\Psi^{\text{T}}=\begin{pmatrix}&1&0&0&0&0&\\\ &0&1&0&0&0&\\\
&0&0&\ell_{1,1}&\ell_{2,1}&a&\\\ &0&0&\ell_{1,2}&\ell_{2,2}&b&\\\
&0&0&\ell_{1,3}&\ell_{2,3}&c&\end{pmatrix}.$ (A.9)

The $5^{\text{th}}$ column of $\Psi^{\text{T}}$ is now fixed (not going to be
updated again) as it corresponds to the detail coefficient but other columns
would be updated as the mergings continue.

##### Second merge

From now on, the “two together” rule is applied. suppose that
$(X_{2},s^{1}_{3,4,5},s^{2}_{3,4,5})$ are selected to be merged at next scale.
Then we firstly select the corresponding weight vectors,
$\boldsymbol{w}^{c}_{2,3,4}=(1,e_{c_{1}},e_{c_{2}})^{\text{T}}$ and
$\boldsymbol{w}^{l}_{2,3,4}=(2,e_{l_{1}},e_{l_{2}})^{\text{T}}$ and construct
the detail filter ${\boldsymbol{h}^{*}}^{\text{T}}=(a^{*},b^{*},c^{*})$ which
satisfies the conditions (1)
${\boldsymbol{h}^{*}}^{\text{T}}\boldsymbol{w}^{c}_{2,3,4}=0$, (2)
${\boldsymbol{h}^{*}}^{\text{T}}\boldsymbol{w}^{l}_{2,3,4}=0$ and (3)
${\boldsymbol{h}^{*}}^{\text{T}}{\boldsymbol{h}^{*}}=1$. The detail filter is
designed as a weight vector for investigating how strong linear trend exist in
$(X_{2},X_{3},X_{4},X_{5})$ as $(e_{c_{1}},e_{c_{2}})$ and
$(e_{l_{1}},e_{l_{2}})$ contain the information of three values
$(X_{3},X_{4},X_{5})$. Two low filters, ${\boldsymbol{\ell}_{1}^{*}}$ and
${\boldsymbol{\ell}_{2}^{*}}$, are obtained by satisfying the conditions,
${\boldsymbol{\ell}_{1}^{*}}^{\text{T}}{\boldsymbol{h}^{*}}=0$,
${\boldsymbol{\ell}_{2}^{*}}^{\text{T}}{\boldsymbol{h}^{*}}=0$,
${\boldsymbol{\ell}_{1}^{*}}^{\text{T}}{\boldsymbol{\ell}_{2}^{*}}=0$ and
${\Lambda^{*}}^{\text{T}}\Lambda^{*}=\mathbf{I}$, and the orthonormal
transform is defined as follows,

$\Lambda^{*}=\begin{pmatrix}&\ell^{*}_{1,1}&\ell^{*}_{1,2}&\ell^{*}_{1,3}&\\\
&\ell^{*}_{2,1}&\ell^{*}_{2,2}&\ell^{*}_{2,3}&\\\
&a^{*}&b^{*}&c^{*}&\end{pmatrix}=\begin{pmatrix}&{\boldsymbol{\ell}_{1}^{*}}^{\text{T}}&\\\
&{\boldsymbol{\ell}_{2}^{*}}^{\text{T}}&\\\
&{\boldsymbol{h}^{*}}^{\text{T}}&\end{pmatrix},$ (A.10)

Now the data sequence and weight vectors are updated as follows,

$\displaystyle\begin{split}\boldsymbol{s}=(X_{1},X_{2},s^{1}_{3,4,5},s^{2}_{3,4,5},d_{3,4,5})&\quad\rightarrow\quad\boldsymbol{s}=(X_{1},s^{1}_{2,2,5},s^{2}_{2,2,5},d_{2,2,5},d_{3,4,5}),\\\
{\boldsymbol{w}_{c}}=(1,1,e_{c_{1}},e_{c_{2}},0)^{\text{T}}&\quad\rightarrow\quad{\boldsymbol{w}_{c}}=(1,e^{*}_{c_{1}},e^{*}_{c_{2}},0,0)^{\text{T}},\\\
{\boldsymbol{w}_{l}}=(1,2,e_{l_{1}},e_{l_{2}},0)^{\text{T}}&\quad\rightarrow\quad{\boldsymbol{w}_{l}}=(1,e^{*}_{l_{1}},e^{*}_{l_{2}},0,0)^{\text{T}},\\\
\end{split}$ (A.11)

and the wavelet matrix is also updated as follows,

$\Psi^{\text{T}}=\begin{pmatrix}1&0&0&0&0\\\
0&\ell^{*}_{1,1}&\ell^{*}_{2,1}&a^{*}&0\\\ \begin{matrix}&\;0&\\\ &\;0&\\\
&\;0&\end{matrix}&\begin{pmatrix}\\\
\ell^{*}_{1,2}\boldsymbol{\ell}_{1}+\ell^{*}_{1,3}\boldsymbol{\ell}_{2}\\\ \\\
\end{pmatrix}&\begin{pmatrix}\\\
\ell^{*}_{2,2}\boldsymbol{\ell}_{1}+\ell^{*}_{2,3}\boldsymbol{\ell}_{2}\\\ \\\
\end{pmatrix}&\begin{pmatrix}\\\
b^{*}\boldsymbol{\ell}_{1}+c^{*}\boldsymbol{\ell}_{2}\\\ \\\
\end{pmatrix}&\begin{matrix}&\;a&\\\ &\;b&\\\ &\;c&\end{matrix}\end{pmatrix}.$
(A.12)

At this scale, the $4^{\text{th}}$ column of $\Psi^{\text{T}}$ is fixed (not
going to be updated again) and it corresponds to the second row of (A.2) whose
non-zero subregion is composed of a single point ($a^{*}$) and a linear trend
in three values ($b^{*}\boldsymbol{\ell}_{1}+c^{*}\boldsymbol{\ell}_{2}$).

Importantly, the orthonormal transform is performed in a way of returning an
orthonormal basis of the expanded subspace e.g. $2^{\text{nd}}$ and
$3^{\text{rd}}$ columns of (A.12) (which are referred to as
${\boldsymbol{\ell}_{1}^{**}}$ and ${\boldsymbol{\ell}_{2}^{**}}$ in (A.13))
are obtained as an arbitrary orthonormal basis of the subspace
$\\{(w,x,y,z)\;|\;w-x=x-y=y-z\\}$ of $\mathbb{R}^{4}$. This is due to the
semi-orthogonality of the transformation matrix $\mathbf{S}$ in (A.13) which
extends the dimension of the subspace from 3 to 4 but preserves the
orthonormality of $\boldsymbol{\ell}_{1}^{*}$ and $\boldsymbol{\ell}_{2}^{*}$
which guarantees the properties,
${\boldsymbol{\ell}_{1}^{**}}^{\text{T}}{\boldsymbol{\ell}_{2}^{**}}=0$ and
${\boldsymbol{\ell}_{1}^{**}}^{\text{T}}{\boldsymbol{\ell}_{1}^{**}}={\boldsymbol{\ell}_{2}^{**}}^{\text{T}}{\boldsymbol{\ell}_{2}^{**}}=1$,
where

${\boldsymbol{\ell}_{1}^{**}}=\begin{pmatrix}\ell^{*}_{1,1}\\\
\begin{pmatrix}\\\
\ell^{*}_{1,2}\boldsymbol{\ell}_{1}+\ell^{*}_{1,3}\boldsymbol{\ell}_{2}\\\ \\\
\end{pmatrix}\end{pmatrix}=\mathbf{S}\begin{pmatrix}\ell^{*}_{1,1}\\\
\ell^{*}_{1,2}\\\
\ell^{*}_{1,3}\end{pmatrix},\;{\boldsymbol{\ell}_{2}^{**}}=\begin{pmatrix}\ell^{*}_{2,1}\\\
\begin{pmatrix}\\\
\ell^{*}_{2,2}\boldsymbol{\ell}_{1}+\ell^{*}_{2,3}\boldsymbol{\ell}_{2}\\\ \\\
\end{pmatrix}\end{pmatrix}=\mathbf{S}\begin{pmatrix}\ell^{*}_{2,1}\\\
\ell^{*}_{2,2}\\\
\ell^{*}_{2,3}\end{pmatrix},\;\mathbf{S}=\begin{pmatrix}1&0&0\\\
\begin{matrix}\;0\\\ \;0\\\ \;0\end{matrix}&\begin{pmatrix}\\\
\boldsymbol{\ell}_{1}\\\ \\\ \end{pmatrix}&\begin{pmatrix}\\\
\boldsymbol{\ell}_{2}\\\ \\\ \end{pmatrix}\end{pmatrix},$ (A.13)

and $\mathbf{S}$ is obtained from the $2^{\text{nd}}$ to $4^{\text{th}}$
columns of (A.9) corresponding to the indices of smooth coefficients merged by
the orthonormal matrix in (A.10).

As the linearity of $\boldsymbol{w}^{c,0}_{3,4,5}$ and
$\boldsymbol{w}^{l,0}_{3,4,5}$ is expressed by a linear span of basis with
unique constants in (A.8), now the extended subregion of the original weight
vectors, $\boldsymbol{w}^{c,0}_{2:5}$ and $\boldsymbol{w}^{l,0}_{2:5}$, can
also be written as a linear span of $\boldsymbol{\ell}_{1}^{**}$ and
$\boldsymbol{\ell}_{2}^{**}$ (an orthonormal basis of the subspace
$\\{(w,x,y,z)\;|\;w-x=x-y=y-z\\}$ of $\mathbb{R}^{4}$) as follows,

$\displaystyle\begin{split}\boldsymbol{w}^{c,0}_{2:5}=e^{*}_{c_{1}}{\boldsymbol{\ell}_{1}^{**}}+e^{*}_{c_{2}}{\boldsymbol{\ell}_{2}^{**}},\quad\boldsymbol{w}^{l,0}_{2:5}=e^{*}_{l_{1}}{\boldsymbol{\ell}_{1}^{**}}+e^{*}_{l_{2}}{\boldsymbol{\ell}_{2}^{**}}.\\\
\end{split}$ (A.14)

This can be simply shown by 1) expressing the weight vectors as a linear span
of two low filters,

$\displaystyle\begin{split}&\boldsymbol{w}^{c}_{2,3,4}=(1,e_{c_{1}},e_{c_{2}})^{\text{T}}=e^{*}_{c_{1}}{\boldsymbol{\ell}_{1}^{*}}+e^{*}_{c_{2}}{\boldsymbol{\ell}_{2}^{*}},\\\
&\boldsymbol{w}^{l}_{2,3,4}=(2,e_{l_{1}},e_{l_{2}})^{\text{T}}=e^{*}_{l_{1}}{\boldsymbol{\ell}_{1}^{*}}+e^{*}_{l_{2}}{\boldsymbol{\ell}_{2}^{*}},\end{split}$
(A.15)

and 2) performing the matrix multiplication with $\mathbf{S}$ in (A.13) to
both sides of (A.15),

$\displaystyle\begin{split}&\text{LHS}:\mathbf{S}\boldsymbol{w}^{c}_{2,3,4}=(1,e_{c_{1}}\boldsymbol{\ell}_{1}+e_{c_{2}}\boldsymbol{\ell}_{2})^{\text{T}}=(1,1,1,1)^{\text{T}}=\boldsymbol{w}^{c,0}_{2:5},\quad\text{RHS}:e^{*}_{c_{1}}{\boldsymbol{\ell}_{1}^{**}}+e^{*}_{c_{2}}{\boldsymbol{\ell}_{2}^{**}},\\\
&\text{LHS}:\mathbf{S}\boldsymbol{w}^{l}_{2,3,4}=(2,e_{l_{1}}\boldsymbol{\ell}_{1}+e_{l_{2}}\boldsymbol{\ell}_{2})^{\text{T}}=(2,3,4,5)^{\text{T}}=\boldsymbol{w}^{l,0}_{2:5},\quad\text{RHS}:e^{*}_{l_{1}}{\boldsymbol{\ell}_{1}^{**}}+e^{*}_{l_{2}}{\boldsymbol{\ell}_{2}^{**}}.\end{split}$
(A.16)

##### Last merge

In the same manner, if the last orthonormal transform is applied to
$(X_{1},s^{1}_{2,3,5},s^{2}_{2,3,5})$, we end up with the finalised
$\Psi^{\text{T}}$ which contains an orthonormal basis of the subspace
$\\{(v,w,x,y,z)\;|\;v-w=w-x=x-y=y-z\\}$ in $\mathbb{R}^{5}$ in its first and
second columns which correspond to two vectors, $\psi^{(0,1)}$ and
$\psi^{(0,2)}$, in (A.5). In other words, regardless of the length of data
$T$, first two wavelet vectors in $\Psi$ which build two final smooth
coefficients ($s^{1}_{1,T},s^{2}_{1,T}$) always have a linear trend with
length $T$ and other wavelet vectors corresponding to detail coefficients
follow one of the forms in (A.2).

Although two low filters ($\boldsymbol{\ell}_{1}$ and $\boldsymbol{\ell}_{2}$
in $\Lambda$) are not uniquely determined at each orthonormal transform, this
has no effect on preserving the linearity of non-zero regions in a wavelet
basis vector as stated above. In simulation study, we empirically find that
the choice of low filters has no effect on the results, thus use a fixed type
of filter function rather than choosing an arbitrary filter every run which
also saves the computational costs.