# Robust Narrowest Significance Pursuit: inference for multiple change-points
in the median

Piotr Fryzlewicz Department of Statistics, London School of Economics,
Houghton Street, London WC2A 2AE, UK. Email<EMAIL_ADDRESS>
###### Abstract

We propose Robust Narrowest Significance Pursuit (RNSP), a methodology for
detecting localised regions in data sequences which each must contain a
change-point in the median, at a prescribed global significance level. RNSP
works by fitting the postulated constant model over many regions of the data
using a new sign-multiresolution sup-norm-type loss, and greedily identifying
the shortest intervals on which the constancy is significantly violated. By
working with the signs of the data around fitted model candidates, RNSP is
able to work under minimal assumptions, requiring only sign-symmetry and
serial independence of the signs of the true residuals. In particular, it
permits their heterogeneity and arbitrarily heavy tails. The intervals of
significance returned by RNSP have a finite-sample character, are
unconditional in nature and do not rely on any assumptions on the true signal.
Code implementing RNSP is available at https://github.com/pfryz/nsp.

Keywords: confidence intervals, structural breaks, post-selection inference,
post-inference selection, narrowest-over-threshold.

##  1 Introduction

There is a growing interest in the statistical literature in the problem of
uncertainty quantification for possibly multiple parameter changes in time- or
space-ordered data, motivated by the practical question of whether any
suspected changes are “real” (i.e. can be attributed to structural changes in
the underlying stochastic model) or due to random fluctuations. General types
of approaches to this problem include confidence sets associated with
simultaneous multiscale change-point estimation (Frick et al., 2014; Pein et
al., 2017), post-selection inference (Hyun et al., 2018, 2021; Jewell et al.,
2020; Duy et al., 2020), inference without selection and post-inference
selection via Narrowest Significance Pursuit (Fryzlewicz, 2021), asymptotic
confidence intervals conditional on the estimated change-point locations (Bai
and Perron, 1998; Eichinger and Kirch, 2018; Cho and Kirch, 2021), False
Discovery Rate (Hao et al., 2013; Li and Munk, 2016; Cheng et al., 2020), and
Bayesian inference (Lavielle and Lebarbier, 2001; Fearnhead, 2006). These
approaches go beyond mere change-point detection and offer statistical
significance statements regarding the existence and locations of change-points
in the statistical model underlying the data.

However, the vast majority of the existing works, with only a handful of
exceptions which we review in more detail later, do so under strong
assumptions on the innovations entering the model, such as their homogeneity
over time, and/or finite moments up to a certain order, and/or a specific
known distribution or distributional family. This is a limitation, because
change-point modelling and inference for data that do not necessarily fall
into these categories is of interest in various domains, including
bioinformatics (Muggeo and Adelfio, 2011; Arlot and Celisse, 2011), biophysics
(Pein et al., 2017), cancer medicine (Kim et al., 2000), economics (Bai and
Perron, 2003) and finance (Oomen, 2019), amongst others.

In this paper, we are concerned with the following problem: given a sequence
of noisy independent observations, automatically determine localised regions
of the data which each must contain a change-point in the median, at a
prescribed global significance level. The methodology we introduce, referred
to as Robust Narrowest Significance Pursuit (RNSP), achieves this for the
piecewise-constant median model, capturing changes in the level of the median.
It does so with practically no distributional assumptions on the data, other
than serial independence of the signs of the true residuals. In contrast to
the existing literature, RNSP requires no knowledge of the distribution on the
part of the user, and permits arbitrarily heavy tails, heterogeneity over time
and/or lack of symmetry. In addition, RNSP is able to handle distributions
that are continuous, continuous with mass points, or discrete, where these
properties may also vary over the signal.

To situate RNSP in the context of the existing literature, we now review the
few existing approaches that permit uncertainty statements for multiple
change-point locations under the assumption of possible heterogeneity in the
distribution of the noise.

H-SMUCE (Pein et al., 2017) is a change-point detector in the heterogeneous
Gaussian piecewise-constant model $Y_{i}=f_{i}+\sigma_{i}\varepsilon_{i}$,
where $f_{i}$ is a piecewise-constant signal, $\varepsilon_{i}$ are i.i.d.
$N(0,1)$ variables, and $\sigma_{i}$ is also piecewise-constant in such a way
that it can only change when $f_{i}$ does. In their Lemma C.7, the authors
provide a theoretical construction of confidence intervals for the locations
of the change-points in $f_{i}$ which works for a fixed significance level
$\alpha$. However, their result assumes that the number of change-points has
not been underestimated, and therefore necessarily assumes certain (and in
this particular case somewhat strong) minimum distance and minimum jump size
conditions, which are typically difficult to verify in practice. In addition,
it offers confidence intervals of the same length regardless of the individual
prominence of each change-point. Finally, the practical use of this result
requires the knowledge of the number of change-points, a lower bound on the
minimum spacing between each pair of change-points, and the minimum jump size,
all of which will often be unknown to the analyst.

However, Pein et al. (2017) also provide a different, practical, algorithmic
construction of confidence intervals, in their Section A.1. This requires
having an estimate of the number of change-points, and involves screening the
data for short intervals over which a constant signal fit is unsuitable and
they must therefore contain change-points. Crucially, this algorithmic
construction relies on the knowledge of scale-dependent critical values (for
measuring the unsuitability of a locally constant fit), which are not
available analytically but only by simulation, and therefore the method would
not extend automatically to unknown noise distributions (as the analyst needs
to know what distribution to sample from). In Section 4, we show that H-SMUCE
suffers from inflated type I error rates in the sense that the thus-
constructed confidence intervals, in the examples of Gaussian models shown, do
not all contain at least one true change-point each in more than
$100(1-\alpha)\%$ of the cases, contrary to what this algorithmic construction
sets out to do. H-SMUCE is also prone to failing if the model is mis-
specified, e.g. if the distribution of the data has a mass point (which is
unsurprising in view of its assumption of Gaussianity).

In contrast to H-SMUCE, RNSP does not assume Gaussianity, and works with any
(possibly heterogeneous) noise distributions, as long as they are sign-
symmetric around the median, a very weak requirement which is immaterial for
continuous distributions (as all, even non-symmetric, continuous distributions
are sign-symmetric). The noise distribution does not need to remain constant
between each consecutive pair of change-points, can change arbitrarily, and
can be continuous, discrete, or continuous with mass points. The execution of
RNSP does not rely on having an estimate of the number of change-points.
Critical values needed by RNSP do not depend on the noise distribution (which
can be unknown to the analyst), and can be accurately approximated
analytically. RNSP explicitly targets the shortest possible intervals of
significance, and its algorithmic construction ensures exact finite-sample
coverage guarantees.

The non-robust predecessor of RNSP, Narrowest Significance Pursuit (NSP;
Fryzlewicz, 2021) is able to handle heterogeneous data in a variety of models,
including the piecewise-constant signal plus independent noise model. However,
NSP requires that the noise, if heterogeneous, satisfies the central limit
theorem (i.e. is within the domain of attraction of the normal distribution)
and is symmetric, neither of which is assumed in RNSP. Crucially, noise
heterogeneity is dealt with in NSP via self-normalisation using the
functional-analytic framework of Rac̆kauskas and Suquet (2003). The self-
normalised statistic used in NSP includes a term resembling an estimate of the
local variance of the noise, which, however, is only unbiased under the null
hypothesis of no change-point being present locally. The fact that the same
term over-estimates the variance under the alternative hypothesis, reduces the
power of the detection statistic, which leads to typically long intervals of
significance. We illustrate this issue with the self-normalised non-robust NSP
for heterogeneous data in Section 4 and show that RNSP offers a significant
improvement.

The Bayesian approach of Fearnhead (2006) provides an inference framework for
independent observations drawn from a $f(\cdot|\theta_{j})$ density within the
$j$th segment, with $\theta_{j}$ possibly vector-valued, which permits some
types of heterogeneity. Operating in this framework requires the choice of
either of two types of priors: one involving a prior on the number of change-
points plus a conditional prior on their locations, and the other – an
instance of the product-partition model. In addition, the distribution family
$f$ is assumed known, and a prior need to be specified for the $\theta_{j}$
parameters. While self-contained and not overly computationally intensive,
this approach suffers from the fact that the estimation of the number of
change-points is a difficult statistical problem and therefore so is choosing
a good prior for this quantity; for the same reason, it is a choice that can
significantly affect the inferential conclusions.

Bai and Perron (1998) and Bai and Perron (2003), in seminal papers studying
least-squares estimation of multiple change-points in regression models under
possible heterogeneity of the errors, describe a procedure for computing
confidence intervals conditional on detection. For an unknown distribution of
the errors, the limiting distribution of each estimated change-point location
converges to an appropriate functional of the Wiener process only under the
assumption that the corresponding break size goes to zero with the sample
size. The asymptotic validity of the resulting confidence interval relies on
the consistency of the estimators of the unknown quantities involved (such as
the local variance of the innovations or the break size); it is therefore a
large-sample, asymptotic construction. Crucially, it does not take into the
account the uncertainty associated with detection, which can be considerable
especially for the more difficult problems (for an example, see the “US ex-
post real interest rate” case study in Bai and Perron (2003), where there is
genuine uncertainty between models with 2 and 3 change-points; we revisit this
example in Section 5.1). By contrast, RNSP produces intervals of significant
change in the median that are not conditional on detection, have a finite-
sample nature and are valid regardless of the size of the breaks.

The paper is organised as follows. Section 2 motivates RNSP and sets out its
general algorithmic framework. Section 3 describes how RNSP measures the local
deviation from model constancy, its crucial ingredient. It also gives
theoretical performance guarantees for RNSP. Section 4 contains numerical
examples and comparisons. Section 5 includes two examples showing the
practical usefulness of RNSP. Section 6 concludes with a brief discussion.
Software implementing RNSP is available at https://github.com/pfryz/nsp.

##  2 Motivation and review of RNSP algorithmic setting

###  2.1 RNSP: context and modus operandi

Robust NSP (RNSP) differs from NSP as defined in Fryzlewicz (2021) in terms of
the data functional it targets, and in how it measures local deviation of this
data functional from the baseline model. RNSP discovers regions in the data in
which the median departs from the baseline model. This is in contrast to NSP,
which targets the mean. In RNSP, we do not make any moment assumptions about
the data, and therefore the median is a natural measure of data centrality. In
RNSP, the baseline model is constant, and therefore the algorithm will look to
discover segments in the data in which the median must experience abrupt
departures from constancy (change-points in the level), at a certain global
significance level.

In this section, we review the components of the algorithmic framework that
are shared between NSP and RNSP, with the generic measurement of local
deviation from the constant model as one of its building blocks. In Section 3,
we will introduce the particular way in which local deviation from the
constant model is measured in RNSP, which is appropriate for the median and
hence fundamentally different from NSP.

RNSP operates in the signal plus noise model

$Y_{t}=f_{t}+Z_{t},\quad t=1,\ldots,T,$ (1)

where the variables $\\{Z_{t}\\}_{t=1}^{T}$ and the signal
$\\{f_{t}\\}_{t=1}^{T}$ satisfy, respectively, the assumptions below. Before
stating them, we define the $\text{sign}(x)$ function as follows.

$\text{sign}(x)=\left\\{\begin{array}[]{cc}1&x>0\\\ 0&x=0\\\
-1&x<0\end{array}\right..$ (2)

###### Assumption 2.1

  1. (a)

The variables $\\{\mathrm{sign}(Z_{t})\\}_{t=1}^{T}$ are mutually independent.

  2. (b)

$\forall\,t$, $M(Z_{t})=0$, where $M$ is the median operator. (If the median
is non-unique, we require $0\in M(Z_{t})$.)

  3. (c)

The variables $\\{Z_{t}\\}_{t=1}^{T}$ are sign-symmetric in the sense that
$P(Z_{t}>0)=P(Z_{t}<0)$, $\forall\,t$.

$\\{Z_{t}\\}_{t=1}^{T}$ do not have to be identically distributed, and can
have arbitrary mass atoms, or none, as long as their distributions satisfy (b)
and (c) in Assumption 2.1. We do not impose any moment assumptions. Any
continuous distribution (even one with an asymmetric density function) is
sign-symmetric. The distribution(s) of $\\{Z_{t}\\}_{t=1}^{T}$ can be unknown
to the analyst. Note that the requirement of the independence of
$\\{\text{sign}(Z_{t})\\}_{t=1}^{T}$ is significantly weaker than that of the
independence of $\\{Z_{t}\\}_{t=1}^{T}$ itself: for example, if $Z_{t}$ is a
(G)ARCH process driven by symmetric, independent innovations, then
$\text{sign}(Z_{t})$ is serially independent, while $Z_{t}$ is not.

###### Assumption 2.2

In (1), $f_{t}$ is a piecewise-constant vector with an unknown number $N$ and
locations $0=\eta_{0}<\eta_{1}<\ldots<\eta_{N}<\eta_{N+1}=T$ of change-points.
(The location $\eta_{j}$ is a change-point if $f_{\eta_{j}}\neq
f_{\eta_{j}+1}$.)

In this work, RNSP does not venture beyond the piecewise-constant model; this
is in contrast to the non-robust NSP (Fryzlewicz, 2021), applicable also to
piecewise-polynomial models of higher degrees, and piecewise-stationary
regression models with arbitrary covariates (including with autoregression).
The reason for it is that for robustness, RNSP works with the signs of the
data – both at the stage of fitting the best constant model on each interval,
and in examining the empirical residuals from the fits. Due to the non-
linearity of the sign operator, RNSP cannot use linear programming for the
constant model fits (unlike the non-robust NSP), and has to manually try all
possible constant model fits on each interval, to be able to select one that
gives the minimum deviation, as this is required for its coverage guarantees.
This is a problem of dimensionality one in the piecewise-constant model, which
makes it computationally straightforward and fast (details are in Section
3.2). However, in models of parameter dimension two or higher (such as the
piecewise-linear model), an analogous exact solution would be too
computationally demanding, hence the restriction to the piecewise-constant
model in the current work. As we explain in detail later on, our coverage
guarantees for RNSP in the piecewise-constant median model are exact and hold
for any sample size.

RNSP achieves the high level of generality in terms of the permitted
distributions of $Z_{t}$ stated in Assumption 2.1 thanks to its use of the
sign transformation. The use of 0 in $\text{sign}(x)$ is critical for RNSP’s
objective to provide correct coverage guarantees also for discrete
distributions and continuous distributions with mass points. We illustrate the
importance of this point in Section 3.2.

###  2.2 (R)NSP: generic algorithm

The generic algorithmic framework underlying both RNSP and the non-robust NSP
in Fryzlewicz (2021) is the same and is based on recursively searching for the
shortest sub-samples in the data with globally significant deviations from the
baseline model. In this section, we introduce this shared generic framework.
In the following sections, we show how RNSP diverges from NSP through its use
of a robust measure of deviation from the baseline model, suitable for the
very broad class of distributions specified in Assumption 2.1.

We start with a pseudocode definition of the RNSP algorithm, in the form of a
recursively defined function RNSP. In its arguments, $[s,e]$ is the current
interval under consideration and at the start of the procedure, we have
$[s,e]=[1,T]$; $Y$ (of length $T$) is as in the model formula (1); $M$ is the
minimum guaranteed number of sub-intervals of $[s,e]$ drawn (unless the number
of all sub-intervals of $[s,e]$ is less than $M$, in which case drawing $M$
sub-intervals would mean repetition); $\lambda_{\alpha}$ is the threshold
corresponding to the global significance level $\alpha$ (typical values for
$\alpha$ would be 0.05 or 0.1) and $\tau_{L}$ (respectively $\tau_{R}$) is a
functional parameter used to specify the degree of overlap of the left
(respectively right) child interval of $[s,e]$ with respect to the region of
significance identified within $[s,e]$, if any. The no-overlap case would
correspond to $\tau_{L}=\tau_{R}\equiv 0$. In each recursive call on a generic
interval $[s,e]$, RNSP adds to the set $\mathcal{S}$ any globally significant
local regions (intervals) of the data identified within $[s,e]$ on which $Y$
is deemed to depart significantly (at global level $\alpha$) from the baseline
constant model. We provide more details underneath the pseudocode below. In
the remainder of the paper, the subscript [s,e] relates to a constant indexed
by the interval $[s,e]$, whose value will be clear from the context.

1:function RNSP($s$, $e$, $Y$, $M$, $\lambda_{\alpha}$, $\tau_{L}$,
$\tau_{R}$)

2: if $e-s<1$ then

3: STOP

4: end if

5: if $M\geq\frac{1}{2}(e-s+1)(e-s)$ then

6: ${M}:=\frac{1}{2}(e-s+1)(e-s)$

7: draw all intervals $[s_{m},e_{m}]\subseteq[s,s+1,\ldots,e]$,
$m=1,\ldots,{M}$, s.t. $e_{m}-s_{m}\geq 1$

8: else

9: draw a representative (see description below) sample of intervals
$[s_{m},e_{m}]\subseteq[s,s+1,\ldots,e]$, $m=1,\ldots,{M}$, s.t.
$e_{m}-s_{m}\geq 1$

10: end if

11: for $m\leftarrow 1,\ldots,M$ do

12: $D_{[s_{m},e_{m}]}:=$ DeviationFromConstantModel$(s_{m},e_{m},Y)$

13: end for

14:
$\mathcal{M}_{0}:=\arg\min_{m}\\{e_{m}-s_{m}\,\,:\,\,m=1,\ldots,M;\,\,D_{[s_{m},e_{m}]}>\lambda_{\alpha}\\}$

15: if $|\mathcal{M}_{0}|=0$ then

16: STOP

17: end if

18:
$m_{0}:=$AnyOf$(\arg\max_{m}\\{D_{[s_{m},e_{m}]}\,\,:\,\,m\in\mathcal{M}_{0}\\})$

19:
$[\tilde{s},\tilde{e}]:=$ShortestSignificantSubinterval$(s_{m_{0}},e_{m_{0}},Y,M,\lambda_{\alpha})$

20: add $[\tilde{s},\tilde{e}]$ to the set $\mathcal{S}$ of significant
intervals

21:
RNSP$(s,\tilde{s}+\tau_{L}(\tilde{s},\tilde{e},Y),Y,M,\lambda_{\alpha},\tau_{L},\tau_{R})$

22:
RNSP$(\tilde{e}-\tau_{R}(\tilde{s},\tilde{e},Y),e,Y,M,\lambda_{\alpha},\tau_{L},\tau_{R})$

23:end function

The RNSP algorithm is launched by the pair of calls below.

$\mathcal{S}:=\emptyset$
RNSP$(1,T,Y,M,\lambda_{\alpha},\tau_{L},\tau_{R})$

On completion, the output of RNSP is in the variable $\mathcal{S}$. We now
comment on the RNSP function line by line. In lines 2–4, execution is
terminated for intervals that are too short. In lines 5–10, a check is
performed to see if $M$ is at least as large as the number of all sub-
intervals of $[s,e]$. If so, then $M$ is adjusted accordingly, and all sub-
intervals are stored in $\\{[s_{m},e_{m}]\\}_{m=1}^{M}$. Otherwise, a sample
of $M$ sub-intervals $[s_{m},e_{m}]\subseteq[s,e]$ is drawn in which either
(a) $s_{m}$ and $e_{m}$ are obtained uniformly and with replacement from
$[s,e]$ (random sampling), or (b) $s_{m}$ and $e_{m}$ are all possible pairs
from an (approximately) equispaced grid on $[s,e]$ which permits at least $M$
such sub-intervals (deterministic sampling).

In lines 11–13, each sub-interval $[s_{m},e_{m}]$ is checked to see to what
extent the response on this sub-interval (denoted by $Y_{s_{m}:e_{m}}$)
conforms to the baseline constant model. This core step of the RNSP algorithm
will be described in more detail in Section 3.

In line 14, the measures of deviation obtained in line 12 are tested against
threshold $\lambda_{\alpha}$, chosen to guaranteed global significance level
$\alpha$. How to choose $\lambda_{\alpha}$ is independent of the distribution
of $Z_{t}$ if it is continuous, and there is also a simple distribution-
independent choice of $\lambda_{\alpha}$ for discrete distributions and
continuous distributions with probability masses; see Section 3.3. The
shortest sub-interval(s) $[s_{m},e_{m}]$ for which the test rejects the
baseline model at global level $\alpha$ are collected in set
$\mathcal{M}_{0}$. In lines 15–17, if $\mathcal{M}_{0}$ is empty, then the
procedure decides that it has not found regions of significant deviations from
the constant model on $[s,e]$, and stops on this interval as a consequence.
Otherwise, in line 18, the procedure continues by choosing the sub-interval,
from among the shortest significant ones, on which the deviation from the
baseline constant model has been the largest. The chosen interval is denoted
by $[s_{m_{0}},e_{m_{0}}]$.

In line 19, $[s_{m_{0}},e_{m_{0}}]$ is searched for its shortest significant
sub-interval, i.e. the shortest sub-interval on which the hypothesis of the
baseline model is rejected locally at a global level $\alpha$. Such a sub-
interval certainly exists, as $[s_{m_{0}},e_{m_{0}}]$ itself has this
property. The structure of this search again follows the workflow of the RNSP
procedure; more specifically, it proceeds by executing lines 2–18 of RNSP, but
with $s_{m_{0}},e_{m_{0}}$ in place of $s,e$. The chosen interval is denoted
by $[\tilde{s},\tilde{e}]$. The importance of this second-stage search in
RNSP’s pursuit to produce short intervals is discussed in detail (in the
context of NSP) in Fryzlewicz (2021). In line 20, the selected interval
$[\tilde{s},\tilde{e}]$ is added to the output set $\mathcal{S}$.

In lines 21–22, RNSP is executed recursively to the left and to the right of
the detected interval $[\tilde{s},\tilde{e}]$. However, we optionally allow
for some overlap with $[\tilde{s},\tilde{e}]$. The overlap, if present, is a
function of $[\tilde{s},\tilde{e}]$ and, if it involves detection of the
location of a change-point within $[\tilde{s},\tilde{e}]$, then it is also a
function of $Y$. An example of the relevance of this is given in Section 5.1.

##  3 Robust NSP: measuring deviation from the constant model

###  3.1 Principles and building blocks of the new deviation measure

The main structure of the DeviationFromConstantModel$(s_{m},e_{m},Y)$
operation is as follows:

  1. 1.

Fit the best, in the sense described precisely later, constant model to
$Y_{s_{m}:e_{m}}$.

  2. 2.

Examine the signs of the empirical residuals from this fit. If their
distribution is deemed to contain a change-point (which indicates that the
constant model fit is unsatisfactory on $[s_{m},e_{m}]$ and therefore the
model contains a change-point on that interval), the value returned by
DeviationFromConstantModel$(s_{m},e_{m},Y)$ should be large; otherwise small.

For RNSP to be able to control, globally and non-trivially, the probability of
spurious detection, the deviation measure it uses must meet the two desiderata
below.

Desideratum A.

If there is no change-point on $[s_{m},e_{m}]$, then the value of
DeviationFromConstantModel$(s_{m},e_{m},Y)$ should be bounded from above by
the same deviation measure for the true model on $[s_{m},e_{m}]$ (the latter
is an oracular quantity, unknown to the analyst), and this in turn should be
bounded from above, uniformly over all $[s_{m},e_{m}]$, by a random variable
involving the signs of the true residuals $Z$ only, such that either its
distribution is known or its quantiles can easily be bounded, sharply, from
above.

Desideratum B.

The deviation measure on $[s_{m},e_{m}]$ cannot be made smaller by proposing
unrealistic constant model fits on that interval; otherwise it would be easy
to force non-detection on any interval. This is an important desired property
as our deviation measure will need to ‘try’ all possible constant model fits
and choose one for which the deviation measure is the smallest, to ensure that
Desideratum A (the part relating to boundedness from above by the deviation
measure for the true model) is met.

As in Fryzlewicz (2021), a key ingredient of our measure of deviation, which
helps it achieve the above desired properties, is a multiresolution sup-norm
introduced below. The fundamental difference with Fryzlewicz (2021) is that
RNSP uses it on the sign of the input rather than in the original data domain.
The basic building block of the multiresolution sup-norm is a scaled partial
sum statistic, defined for an arbitrary input sequence $\\{x_{t}\\}_{t=1}^{T}$
by

$U_{s,e}(x)=\frac{1}{(e-s+1)^{1/2}}\sum_{t=s}^{e}x_{t}.$ (3)

We define the multiresolution sup-norm (Nemirovski, 1986; Li, 2016) of an
input vector $x$ (of length $T$) with respect to the interval set
$\mathcal{I}$ as

$\|x\|_{\mathcal{I}}=\max_{[s,e]\in\mathcal{I}}|U_{s,e}(x)|.$ (4)

The set $\mathcal{I}$ used in RNSP contains intervals at a range of scales and
locations. A canonical example of a suitable interval set $\mathcal{I}$ is the
set $\mathcal{I}^{a}$ of all subintervals of $[1,T]$. We will use
$\mathcal{I}^{a}$ in defining the largest acceptable global probability of
spurious detection. However, for computational reasons,
DeviationFromConstantModel will use a smaller interval set (we give the
details below). This will not affect the exactness of our coverage guarantees,
because, naturally, if $\mathcal{J}\subseteq\mathcal{I}$, then
$\|x\|_{\mathcal{J}}\leq\|x\|_{\mathcal{I}}$. We also define the restriction
of $\mathcal{I}$ to an arbitrary interval $[s,e]$ as

$\mathcal{I}_{[s,e]}=\\{[u,v]\subseteq[s,e]\,\,:\,\,[u,v]\in\mathcal{I}\\}.$

Note the trivial inequality

$\|x_{s:e}\|_{\mathcal{I}^{a}_{[s,e]}}\leq\|x\|_{\mathcal{I}^{a}}$ (5)

for any $[s,e]\subseteq[1,T]$.

There are a number of reasons for the use of the multiresolution sup-norm in
our deviation measure. One is that thanks to the particular construction of
the RNSP algorithm and our deviation measure, in order for us to be able to
control spurious detections at a global level $\alpha$, it suffices to know,
or be able to bound from above (as sharply as possible), the $1-\alpha$
quantile of the distribution of $\|\text{sign}(Z)\|_{\mathcal{I}^{a}}$.
However, relatively sharp results of this type are either well-known or easy
to obtain, as we show in Section 3.3.

Another reason is that the multiresolution sup-norm of the signs of the
empirical residuals from a postulated baseline model fit is naturally large
for proposed model fits if they are unrealistic, which penalises them, thereby
ensuring that Desideratum B above is met. This important property would not be
available if, instead of the multiresolution sup-norm of the residual signs,
we chose to use, for example, maximum-likelihood or other contrast-based
approaches to test deviations from the baseline (the simplest example of such
a test is the CUSUM test). We discuss this point in more detail in Section
3.2.1, where we also give further reasons supporting the use of the
multiresolution sup-norm in our context.

In line with Desideratum A, in the following sections we define
$D_{[s_{m},e_{m}]}:=$ DeviationFromConstantModel$(s_{m},e_{m},Y)$ which
satisfies the property that if there is no change-point on the interval
$[s_{m},e_{m}]$, then it is guaranteed that

$D_{[s_{m},e_{m}]}\leq\|\text{sign}(Z_{s_{m}:e_{m}})\|_{\mathcal{I}^{a}_{[s_{m},e_{m}]}}.$
(6)

###  3.2 Deviation measure: definition and properties

We now define the deviation measure $D_{[s_{m},e_{m}]}$, returned by the
DeviationFromConstantModel function from line 12 of the RNSP algorithm.

The discussion below assumes that there is no change-point in $[s_{m},e_{m}]$.
For the true constant signal $f_{s_{m}:e_{m}}$, denote
$f_{[s_{m},e_{m}]}:=f_{s_{m}}=\ldots=f_{e_{m}}$. If we knew $f_{s_{m}:e_{m}}$,
we could use
$\|\text{sign}(Y_{s_{m}:e_{m}}-f_{[s_{m},e_{m}]})\|_{\mathcal{I}^{a}_{[s_{m},e_{m}]}}$
as $D_{[s_{m},e_{m}]}$, since we trivially have

$\|\text{sign}(Y_{s_{m}:e_{m}}-f_{[s_{m},e_{m}]})\|_{\mathcal{I}^{a}_{[s_{m},e_{m}]}}=\|\text{sign}(Z_{s_{m}:e_{m}})\|_{\mathcal{I}^{a}_{[s_{m},e_{m}]}},$

which satisfies the desired condition (6). Of course $f_{[s_{m},e_{m}]}$ is
not known, but a key observation is that there are only at most
$2(s_{m}-e_{m})+3$ different possible constants $\tilde{f}_{[s_{m},e_{m}]}$
leading to different sequences
$\\{\text{sign}(Y_{t}-\tilde{f}_{[s_{m},e_{m}]})\\}_{t=s_{m}}^{e_{m}}$. To see
this, sort the values of $Y_{s_{m}:e_{m}}$ in non-decreasing order to create
$Y_{(1)},Y_{(2)},\ldots,Y_{(e_{m}-s_{m}+1)}$. Take candidate constants
$\tilde{f}_{[s_{m},e_{m}]}^{\\{j\\}}$, $j=1,\ldots,2(s_{m}-e_{m})+3$, defined
as follows.

$\displaystyle\tilde{f}_{[s_{m},e_{m}]}^{\\{1\\}}$ $\displaystyle<$
$\displaystyle Y_{(1)}\quad\text{(but otherwise arbitrary)}$
$\displaystyle\tilde{f}_{[s_{m},e_{m}]}^{\\{2\\}}$ $\displaystyle=$
$\displaystyle Y_{(1)}$ $\displaystyle\tilde{f}_{[s_{m},e_{m}]}^{\\{3\\}}$
$\displaystyle=$ $\displaystyle\frac{1}{2}(Y_{(1)}+Y_{(2)})$
$\displaystyle\tilde{f}_{[s_{m},e_{m}]}^{\\{4\\}}$ $\displaystyle=$
$\displaystyle Y_{(2)}$ $\displaystyle\tilde{f}_{[s_{m},e_{m}]}^{\\{5\\}}$
$\displaystyle=$ $\displaystyle\frac{1}{2}(Y_{(2)}+Y_{(3)})$
$\displaystyle\vdots$
$\displaystyle\tilde{f}_{[s_{m},e_{m}]}^{\\{2(e_{m}-s_{m})+2\\}}$
$\displaystyle=$ $\displaystyle Y_{(e_{m}-s_{m}+1)}$
$\displaystyle\tilde{f}_{[s_{m},e_{m}]}^{\\{2(e_{m}-s_{m})+3\\}}$
$\displaystyle>$ $\displaystyle Y_{(e_{m}-s_{m}+1)}\quad\text{(but otherwise
arbitrary)}.$ (7)

We have the following simple result.

###### Proposition 3.1

Assume no change-point in $[s_{m},e_{m}]$ and denote
$f_{[s_{m},e_{m}]}:=f_{s_{m}}=\ldots=f_{e_{m}}$. Let the constants
$\tilde{f}_{[s_{m},e_{m}]}^{\\{j\\}}$, $j=1,\ldots,2(s_{m}-e_{m})+3$ be
defined as in (7). There exists a $j_{0}\in\\{1,\ldots,2(s_{m}-e_{m})+3\\}$
such that

$\\{\mathrm{sign}(Y_{t}-f_{[s_{m},e_{m}]})\\}_{t=s_{m}}^{e_{m}}=\\{\mathrm{sign}(Y_{t}-\tilde{f}_{[s_{m},e_{m}]}^{\\{j_{0}\\}})\\}_{t=s_{m}}^{e_{m}}.$
(8)

Proof. If there is a $j_{1}$ such that $f_{[s_{m},e_{m}]}=Y_{(j_{1})}$, then
choose $j_{0}=2j_{1}$ so that
$\tilde{f}_{[s_{m},e_{m}]}^{\\{j_{0}\\}}=Y_{(j_{1})}=f_{[s_{m},e_{m}]}$ and
(8) is trivially satisfied. Otherwise, if there is a $j_{1}$ such that
$f_{[s_{m},e_{m}]}\in(Y_{(j_{1})},Y_{(j_{1}+1)})$, then choose
$j_{0}=2j_{1}+1$ so that

$\displaystyle\\{\mathrm{sign}(Y_{t}-f_{[s_{m},e_{m}]})\\}_{t=s_{m}}^{e_{m}}$
$\displaystyle=$
$\displaystyle\\{\mathrm{sign}(Y_{t}-(Y_{(j_{1})}+Y_{(j_{1}+1)})/2)\\}_{t=s_{m}}^{e_{m}}$
$\displaystyle=$
$\displaystyle\\{\mathrm{sign}(Y_{t}-\tilde{f}_{[s_{m},e_{m}]}^{\\{j_{0}\\}})\\}_{t=s_{m}}^{e_{m}},$

and (8) is satisfied. Otherwise, it must be that either
$f_{[s_{m},e_{m}]}<Y_{(1)}$ or $f_{[s_{m},e_{m}]}>Y_{(e_{m}-s_{m}+1)}$; in the
former case take $j_{0}=1$ and in the latter, $j_{0}=2(e_{m}-s_{m})+3$ and (8)
is satisfied by a similar argument. $\square$

We now define our measure of deviation $D_{[s_{m},e_{m}]}$, and prove its key
property as a corollary to Proposition 3.1.

###### Definition 3.1

Let the constants $\tilde{f}_{[s_{m},e_{m}]}^{\\{j\\}}$,
$j=1,\ldots,2(s_{m}-e_{m})+3$ be defined as in (7). We define

$D_{[s_{m},e_{m}]}:=\min_{j\in\\{1,\ldots,2(s_{m}-e_{m})+3\\}}\|\mathrm{sign}(Y_{s_{m}:e_{m}}-\tilde{f}^{\\{j\\}}_{[s_{m},e_{m}]})\|_{\mathcal{I}^{a}_{[s_{m},e_{m}]}}.$
(9)

Essentially, $D_{[s_{m},e_{m}]}$ tries all possible baseline constant model
fits on $[s_{m},e_{m}]$ and chooses the one for which the multiresolution norm
of the signs of the residuals from the fit is the smallest. Choosing the fit
that minimises the multiresolution fit is essential for guaranteeing coverage
properties, as we will see below. If there is a change-point on
$[s_{m},e_{m}]$, the hope is that even the minimum multiresoltion norm fit as
returned by $D_{[s_{m},e_{m}]}$ will be large enough for RNSP to indicate a
change-point on that interval.

###### Corollary 3.1

For any interval $[s_{m},e_{m}]$ on which there is no change-point, we have

$D_{[s_{m},e_{m}]}\leq\|\mathrm{sign}(Z_{s_{m}:e_{m}})\|_{\mathcal{I}^{a}_{[s_{m},e_{m}]}}.$
(10)

In other words, the deviation measure defined in (9) satisfies the desired
property (6).

Proof. Let the index $j_{0}$ be as in the statement of Proposition 3.1. We
have

$\displaystyle D_{[s_{m},e_{m}]}$ $\displaystyle=$
$\displaystyle\min_{j\in\\{1,\ldots,2(s_{m}-e_{m})+3\\}}\|\mathrm{sign}(Y_{s_{m}:e_{m}}-\tilde{f}^{\\{j\\}}_{[s_{m},e_{m}]})\|_{\mathcal{I}^{a}_{[s_{m},e_{m}]}}$
$\displaystyle\leq$
$\displaystyle\|\mathrm{sign}(Y_{s_{m}:e_{m}}-\tilde{f}^{\\{j_{0}\\}}_{[s_{m},e_{m}]})\|_{\mathcal{I}^{a}_{[s_{m},e_{m}]}}$
$\displaystyle=$
$\displaystyle\|\mathrm{sign}(Y_{s_{m}:e_{m}}-f_{[s_{m},e_{m}]})\|_{\mathcal{I}^{a}_{[s_{m},e_{m}]}}$
$\displaystyle=$
$\displaystyle\|\mathrm{sign}(Z_{s_{m}:e_{m}})\|_{\mathcal{I}^{a}_{[s_{m},e_{m}]}}.$

$\square$

This leads to the following guarantee for the RNSP algorithm.

###### Theorem 3.1

Let $\mathcal{S}=\\{S_{1},\ldots,S_{R}\\}$ be the set of intervals returned by
the RNSP algorithm. The following guarantee holds.

$P\left(\exists\,\,{i=1,\ldots,R}\,\,\,\forall\,\,{j=1,\ldots,N}\,\,\,[\eta_{j},\eta_{j}+1]\not\subseteq
S_{i}\right)\leq P(\|\mathrm{sign}(Z)\|_{\mathcal{I}^{a}}>\lambda_{\alpha}).$

Proof. On the set
$\|\mathrm{sign}(Z)\|_{\mathcal{I}^{a}}\leq\lambda_{\alpha}$, each interval
$S_{i}$ must contain a change-point as if it did not, then by Corollary 3.1
and inequality (5), we would have to have

$D_{S_{i}}\leq\|\mathrm{sign}(Z)\|_{\mathcal{I}^{a}}\leq\lambda_{\alpha}.$
(11)

However, the fact that $S_{i}$ was returned by RNSP means, by line 14 of the
RNSP algorithm, that $D_{S_{i}}>\lambda_{\alpha}$, which contradicts (11).
This completes the proof. $\square$

Theorem 3.1 should be read as follows. Let
$\alpha=P(\|\mathrm{sign}(Z)\|_{\mathcal{I}^{a}}>\lambda_{\alpha})$. For a set
of intervals returned by RNSP, we are guaranteed, with probability of at least
$1-\alpha$, that there is at least one change-point in each of these
intervals. Therefore, $\mathcal{S}=\\{S_{1},\ldots,S_{R}\\}$ can be
interpreted as an automatically chosen set of regions (intervals) of
significance in the data. In the no-change-point case ($N=0$), the correct
reading of Theorem 3.1 is that the probability of obtaining one of more
intervals of significance ($R\geq 1$) is bounded from above by
$P(\|\mathrm{sign}(Z)\|_{\mathcal{I}^{a}}>\lambda_{\alpha})$.

We emphasise that Theorem 3.1 does not promise to detect all the change-
points, or to do so asymptotically as the sample size gets larger: this would
be impossible without assumptions on the strength of the change-points
(involving spacing between neighbouring change-points and the sizes of the
jumps). Such assumptions are typically impossible to verify in practice and we
do not make them in this work. Instead, Theorem 3.1 promises that every
interval of significance returned by RNSP must contain at least one change-
point each, with a certain global probability adjustable by the user.

The intervals of significance returned by RNSP have an “unconditional
confidence interval” interpretation: they are not conditional on any prior
detection event, but indicate regions in the data each of which must
unconditionally contain at least one change in the median of $Y_{t}$, with a
global probability of at least $1-\alpha$. Therefore, as in NSP (Fryzlewicz,
2021), RNSP can be viewed as performing “inference without selection” (where
“inference” refers to producing the RNSP intervals of significance and
“selection” to the estimation of change-point locations, absent from RNSP).
This viewpoint also enables “post-inference selection” or “in-inference
selection” if the exact change-point locations (if any) are to be estimated
within the RNSP intervals of significance after or during the execution of
RNSP. See the discussion in Fryzlewicz (2021) for a further discussion of
these aspects. We also emphasise that the statement of Theorem 3.1 is of a
finite-sample character and holds exactly and for any sample size.

####  3.2.1 Deviation measure: discussion

We now comment on a number of aspects of the deviation measure $D_{\cdot}$.

Method of computation. In the non-robust NSP (Fryzlewicz, 2021), the analogue
of $D_{[s_{m},e_{m}]}$ was amenable to computation via linear programming; the
non-robust NSP with self-normalisation also used linear programming as the
main step in its computation. By contrast, in RNSP, measuring the deviation
from the baseline model involves examining the signs of the residuals rather
than the residuals themselves as was the case in NSP. Because of the non-
linearity of the sign function, it is not possible to use linear programming
for the computation of $D_{[s_{m},e_{m}]}$ in RNSP, and (9) needs to be
computed by manually trying out all candidate constants
$\tilde{f}^{\\{j\\}}_{[s_{m},e_{m}]}$. While this may appear expensive, it is
the only option, as (a) trying only ‘realistic’ constant model fits would
require an additional global parameter describing what it means for a local
constant fit to be realistic, and (b) carrying out the fitting in an ad hoc
way, e.g. by only fitting the empirical median of $Y_{s_{m}:e_{m}}$ (rather
than each of the constants $\tilde{f}^{\\{j\\}}_{[s_{m},e_{m}]}$), would
violate the desired property (6) and therefore not be able to lead to exact
coverage guarantees.

Achieving computational savings without affecting coverage guarantees.
Although the operation of trying each constant
$\tilde{f}^{\\{j\\}}_{[s_{m},e_{m}]}$ in (9) is relatively fast, in order to
accelerate it further, we introduce the two computational savings below, which
do not increase $D_{[s_{m},e_{m}]}$ and therefore respect the crucial
inequality (10) and hence also our coverage guarantees in Theorem 3.1.

Reducing the set $\mathcal{I}^{a}_{[s_{m},e_{m}]}$.

To accelerate the computation of (9), we replace the set
$\mathcal{I}^{a}_{[s_{m},e_{m}]}$ in $D_{[s_{m},e_{m}]}$ with the set
$\mathcal{I}^{lr}_{[s_{m},e_{m}]}$ (with $l$ and $r$ standing for left and
right, respectively), defined as follows.

$\mathcal{I}^{lr}_{[s_{m},e_{m}]}:=\mathcal{I}^{l}_{[s_{m},e_{m}]}\cup\mathcal{I}^{r}_{[s_{m},e_{m}]},$

where

$\displaystyle\mathcal{I}^{l}_{[s_{m},e_{m}]}$ $\displaystyle=$
$\displaystyle\\{[s_{m},s_{m}+1],[s_{m},s_{m}+2],\ldots,[s_{m},e_{m}]\\}$
$\displaystyle\mathcal{I}^{r}_{[s_{m},e_{m}]}$ $\displaystyle=$
$\displaystyle\\{[s_{m},e_{m}],[s_{m}+1,e_{m}],\ldots,[e_{m}-1,e_{m}]\\}.$

This reduces the cardinality of the set of intervals included in
$D_{[s_{m},e_{m}]}$ from $O((e_{m}-s_{m})^{2})$ to $O(e_{m}-s_{m})$. As
$\mathcal{I}^{lr}_{[s_{m},e_{m}]}\subseteq\mathcal{I}^{a}_{[s_{m},e_{m}]}$ and
hence
$\|\cdot\|_{\mathcal{I}^{lr}_{[s_{m},e_{m}]}}\leq\|\cdot\|_{\mathcal{I}^{a}_{[s_{m},e_{m}]}}$,
the results of Corollary 3.1 and Theorem 3.1 remain unchanged for the thus-
reduced $D_{[s_{m},e_{m}]}$. On the other hand,
$\mathcal{I}^{lr}_{[s_{m},e_{m}]}$ has been defined in this particular way so
as not to compromise detection power in the piecewise-constant median model.
To see this, consider the following illustrative example. Suppose
$Y_{t}=f_{t}$ (noiseless case) and $f_{t}=0$ for $t=1,\ldots,50$ and $f_{t}=1$
for $t=51,\ldots,100$. On $[s_{m},e_{m}]=[1,100]$, the baseline constant
signal level fitted is $\tilde{f}_{[1,100]}=1/2$ and we have

$\mbox{sign}(Y_{t}-\tilde{f}_{[1,100]})=\left\\{\begin{array}[]{cc}-1&t=1,\ldots,50\\\
1&t=51,\ldots,100.\end{array}\right.$

In this setting, the two multiresoltion sup-norms:

$\|\mbox{sign}(Y_{t}-\tilde{f}_{[1,100]})\|_{\mathcal{I}^{lr}_{[1,100]}}\quad\text{and}\quad\|\mbox{sign}(Y_{t}-\tilde{f}_{[1,100]})\|_{\mathcal{I}^{a}_{[1,100]}}$

are identical, equal to $\sqrt{50}$ and achieve this value for the intervals
$[1,50]$ and $[51,100]$, members of both $\mathcal{I}^{a}_{[1,100]}$ and
$\mathcal{I}^{lr}_{[1,100]}$. This simple example illustrates the wider
phenomenon that if there is a single change-point in the median on a generic
interval $[s_{m},e_{m}]$ under consideration, then in the noiseless case the
multiresolution norm over the set $\mathcal{I}^{a}_{[s_{m},e_{m}]}$ is
maximised at one of the “left” or “right” intervals in
$\mathcal{I}^{lr}_{[s_{m},e_{m}]}$, and we are happy to sacrifice potential
negligible differences in the noisy case in exchange for the substantial
computational savings.

Limiting interval lengths.

In practice, the analyst may not be interested in excessively long RNSP
intervals of significance, or in other words in change-points so weak that
only a very long interval around them is required to be able to detect them.
With this in mind, our software at https://github.com/pfryz/nsp provides an
option of limiting the length of intervals considered in the sense that
$D_{[s_{m},e_{m}]}$ is able to automatically return zero if $e_{m}-s_{m}>L$
for a user-specified maximum length $L$.

Unsuitability of CUSUM or similar contrasts in deviation measure. We note that
it would be impossible to replace the multiresolution sup-norm in
$D_{[s_{m},e_{m}]}$ by, for example, the CUSUM statistic or a similar contrast
measure such as that described in Ellenberger et al. (2021). Consider, for
example, the hypothetical definition of a deviation measure as follows:

$D_{[s_{m},e_{m}]}^{inv}:=\min_{j\in\\{1,\ldots,2(s_{m}-e_{m})+3\\}}|\mbox{CUSUM}\\{\mathrm{sign}(Y_{s_{m}:e_{m}}-\tilde{f}^{\\{j\\}}_{[s_{m},e_{m}]})\\}|,$

where “inv” stands for invalid. The definition of the CUSUM statistic is well-
established, see e.g. Fryzlewicz (2014) for details. This is an invalid
definition as
$\mathrm{sign}(Y_{s_{m}:e_{m}}-\tilde{f}^{\\{1\\}}_{[s_{m},e_{m}]})$ is a
vector of ones, and the CUSUM statistic returns zero for constant vectors, so
$D_{[s_{m},e_{m}]}^{inv}$ would not be able to offer detection under any
circumstances as it would always equal zero. This is an example of a
construction that violates Desideratum B.

Importance of zero in sign function. For valid coverage guarantees in the
presence of mass points in $Z_{t}$ (or if the distribution of $Z_{t}$ is
discrete), it is crucial for the sign function defined in (2) to return zero
if its argument is zero. To illustrate this point, define

$\text{sign}^{inv}(x)=\left\\{\begin{array}[]{cc}1&x\geq 0\\\
-1&x<0\end{array}\right.,$ (12)

where “inv” stands for “invalid”, and consider the trivial case $Z_{t}\equiv
0$. For no-change-point input data
$Y_{s_{m}:e_{m}}=(f_{[s_{m},e_{m}]},\ldots,f_{[s_{m},e_{m}]})$, there are only
three different constants $\tilde{f}_{[s_{m},e_{m}]}^{\\{j\\}}$ (see
definition (7)). With $\text{sign}^{inv}(x)$ in place of $\text{sign}(x)$,
this would lead to

$\displaystyle\\{\mathrm{sign}^{inv}(Y_{t}-\tilde{f}_{[s_{m},e_{m}]}^{\\{1\\}})\\}_{t=s_{m}}^{e_{m}}$
$\displaystyle=$ $\displaystyle(1,\ldots,1)$
$\displaystyle\\{\mathrm{sign}^{inv}(Y_{t}-\tilde{f}_{[s_{m},e_{m}]}^{\\{2\\}})\\}_{t=s_{m}}^{e_{m}}$
$\displaystyle=$ $\displaystyle(1,\ldots,1)$
$\displaystyle\\{\mathrm{sign}^{inv}(Y_{t}-\tilde{f}_{[s_{m},e_{m}]}^{\\{3\\}})\\}_{t=s_{m}}^{e_{m}}$
$\displaystyle=$ $\displaystyle(-1,\ldots,-1)$

and therefore we would have $D_{[s_{m},e_{m}]}=\sqrt{e_{m}-s_{m}+1}$, the
largest value $D_{[s_{m},e_{m}]}$ can possibly take, which would therefore
have to lead to the (spurious) designation of $[s_{m},e_{m}]$ as containing a
change-point, for $e_{m}-s_{m}$ suitably large. Naturally, the analogous
argument would also apply if $\text{sign}^{inv}(0)=-1$ rather than 1. By
contrast, note that the use of the (correct) sign function leads to

$\\{\mathrm{sign}(Y_{t}-\tilde{f}_{[s_{m},e_{m}]}^{\\{2\\}})\\}_{t=s_{m}}^{e_{m}}=(0,\ldots,0),$

which yields $D_{[s_{m},e_{m}]}=0$ and therefore the (correct) designation of
$[s_{m},e_{m}]$ as not containing a change-point.

Importance of placing $\tilde{f}_{[s_{m},e_{m}]}^{\\{j\\}}$ at data points.
For the same reason, it is important that the set of test levels
$\tilde{f}_{[s_{m},e_{m}]}^{\\{j\\}}$ (definition (7)) includes those placed
at the data points themselves (i.e. those indexed by even values of $j$ in
definition (7)). Indeed, continuing the example directly above, if we were to
exclude the constant $\tilde{f}_{[s_{m},e_{m}]}^{\\{2\\}}$ from the list of
test levels considered, we would then have
$D_{[s_{m},e_{m}]}=\sqrt{e_{m}-s_{m}+1}$ even with the use of the correct
function sign (and not only with $\text{sign}^{inv}$), which would again lead
to spurious detection.

Importance of placing $\tilde{f}_{[s_{m},e_{m}]}^{\\{j\\}}$ in between sorted
data points. It is equally important that the test levels
$\tilde{f}_{[s_{m},e_{m}]}^{\\{j\\}}$ should include those placed in between
the sorted data points (i.e. those indexed by odd values of $j$ in definition
(7)). This can be see e.g. by considering $Z_{t}$ such that
$P(Z_{t}=1)=P(Z_{t}=-1)=1/2$; for brevity, we omit the full discussion.

###  3.3 Evaluation and bounds for $\|\text{sign}(Z)\|_{\mathcal{I}^{a}}$

To make Theorem 3.1 operational, we need to obtain an understanding of the
distribution of $\|\text{sign}(Z)\|_{\mathcal{I}^{a}}$ so we are able to
choose $\lambda_{\alpha}$ such that
$P(\|\text{sign}(Z)\|_{\mathcal{I}^{a}}>\lambda_{\alpha})=\alpha$ (or
approximately so) for a desired global significance level $\alpha$.

Initially we consider $Z_{t}$ such that
$P(\text{sign}(Z_{t})=1)=P(\text{sign}(Z_{t})=-1)=1/2$, i.e.
$P(\text{sign}(Z_{t})=0)=0$ for all $t$ (the general case
$P(\text{sign}(Z_{t})=0)\geq 0$ is covered in the next paragraph). From
Theorem 1.1 in Kabluchko and Wang (2014), we have that

$\lim_{T\to\infty}P(\|\text{sign}(Z)\|_{\mathcal{I}^{a}}>a_{T}+\tau/a_{T})=1-\exp(-2\Lambda\exp(-\tau)),$
(13)

where $a_{T}=\\{2\log(T\log^{-1/2}T)\\}^{1/2}$ and $\Lambda$ is a constant. As
Kabluchko and Wang (2014) do not unambiguously state the value of $\Lambda$,
we use simulation over a range of values of $T$ and $\tau$ to determine a
suitable value of $\Lambda$ as 0.274. The practical choice of the significance
threshold $\lambda_{\alpha}$ then proceed as follows: (a) fix $\alpha$ to the
desired level, for example 0.05 or 0.1; (b) obtain the value of $\tau$ by
equating $1-\exp(-2\Lambda\exp(-\tau))=\alpha$; (c) obtain
$\lambda_{\alpha}=a_{T}+\tau/a_{T}$.

Suppose now that $P(\text{sign}(Z_{t})=0)=\rho_{t}\geq 0$; note that the sign-
symmetry Assumption 2.1(c) implies
$P(\text{sign}(Z_{t})=1)=P(\text{sign}(Z_{t})=-1)=(1-\rho_{t})/2$. Construct
the variable $\tilde{Z}_{t}=Z_{t}\,\,|\,\,Z_{t}\neq 0$. As
$P(\text{sign}(\tilde{Z}_{t})=1)=P(\text{sign}(\tilde{Z}_{t})=-1)=1/2$, the
limiting statement (13) applies to $\text{sign}(\tilde{Z}_{t})$. However, we
have the double inequality

$\|\text{sign}(Z)\|_{\mathcal{I}^{a}}\leq\|\text{sign}(\tilde{Z})\|_{\mathcal{I}^{a}_{I}}\leq\|\text{sign}(\tilde{Z})\|_{\mathcal{I}^{a}},$
(14)

with $I=[1,2,\ldots,T_{1}]$, where
$T_{1}=|\\{t\in[1,\ldots,T]\,\,:\,\,Z_{t}\neq 0\\}|$. The first inequality in
(14) holds because every constituent partial sum of
$\|\text{sign}(Z)\|_{\mathcal{I}^{a}}$ has a corresponding larger or equal in
magnitude partial sum in $\|\text{sign}(\tilde{Z})\|_{\mathcal{I}^{a}_{I}}$
constructed by removing the zeros from its numerator and decreasing (or not
increasing) its denominator as it contains fewer (or the same number of)
terms. As an illustrative example, suppose the sequence of
$\text{sign}(Z_{t})$ starts $-1,0,1,1$. The absolute partial sum
$|-1+0+1+1|/\sqrt{4}$, a constituent of
$\|\text{sign}(Z)\|_{\mathcal{I}^{a}}$, is majorised by the absolute partial
sum $|-1+1+1|/\sqrt{3}$, a constituent of
$\|\text{sign}(\tilde{Z})\|_{\mathcal{I}^{a}_{I}}$, where the latter sum has
been constructed by removing the 0 from $-1,0,1,1$ and adjusting for the
number of terms (now 3 instead of 4). The second inequality in (14) holds
simply because $T_{1}\leq T$.

The implication of (14) is that $\|\text{sign}(Z)\|_{\mathcal{I}^{a}}$ for
$\rho_{t}\geq 0$ is majorised by $\|\text{sign}(Z)\|_{\mathcal{I}^{a}}$ for
$\rho_{t}=0$, the case handled by (13). Therefore, the threshold
$\lambda_{\alpha}$ obtain as a consequence of (13) can also be applied in the
general case $\rho_{t}\geq 0$.

##  4 Numerical illustrations

In this section, we demonstrate numerically that the guarantee offered by
Theorem 3.1 holds for RNSP in practice over a variety of homogeneous and
heterogeneous models for which the variables $Z_{t}$ satisfy Assumption 2.1.
We also investigate the circumstances under which similar guarantees are not
offered by H-SMUCE (Pein et al., 2017) or by the self-normalised version of
NSP (SN-NSP), suitable for heterogeneous data (Fryzlewicz, 2021). In this
section, we use the acronyms RNSP and SN-NSP to denote the versions of these
respective procedures with no interval overlaps, i.e. $\tau_{L}=\tau_{R}=0$.
Later in this section, we introduce notation for versions with overlaps. Both
RNSP and SN-NSP use $M=1000$ intervals, the default setting. For H-SMUCE, the
function call we use is `stepR::stepFit(x, alpha = 0.1, family="hsmuce",
confband=TRUE)`.

We begin with null models, by which we mean models (1) for which $f_{t}$ is
constant throughout, i.e. $N=0$. For null models, Theorem 3.1 promises that
RNSP at level $\alpha$ returns no intervals of significance with probability
at least $1-\alpha$. In this section, we use $\alpha=0.1$. There are similar
parameters in H-SMUCE and SN-NSP, and they are also set to 0.1. The null
models used are listed in Table 1.

model name | sample path execution in R
---|---
Plain Gauss | `rnorm(100)`
Plain Poisson | `as.numeric(rpois(200, 1))`
Heterogeneous Gauss | `c(rep(1, 100), rep(8, 50), rep(1, 100)) * rnorm(250)`
Symmetric Bernoulli | `as.numeric(rbinom(200, 1, .5))`
Plain Cauchy | `rcauchy(100, 0)`
Mix 1 | `sample(3, size=300, replace=TRUE, prob=c(.35, .3, .35)) -> xx`
| `xx[xx != 2] <- rnorm(sum(xx !=2 ))`
Mix 2 | `rpois(200, 5)+rnorm(200)/30`
Table 1: Null models for the comparative simulation study in Section 4. model | RNSP | H-SMUCE | SN-NSP
---|---|---|---
Plain Gauss | 100 | 98 | 100
Plain Poisson | 99 | 0 | 1
Heterogeneous Gauss | 99 | 97 | 94
Symmetric Bernoulli | 90 | 0 | 30
Plain Cauchy | 100 | 100 | 100
Mix 1 | 100 | 0 | 99
Mix 2 | 100 | 82 | 100
Table 2: Numbers of times, out of 100 simulated sample paths of each null
model, that the respective method indicates no intervals of significance.

Table 2 shows the null model results. RNSP keeps the nominal size well across
all the models considered, returning no intervals of significance at least 90%
of the time in all situations. The lowest (but still correct) coverage of
exactly 90% is achieved for the Symmetric Binomial model: this is because that
distribution admits only two values ($-1$, $1$), and therefore the single
level $\tilde{f}^{\\{3\\}}_{[s_{m},e_{m}]}=0$ (see Section 3.2) always
achieves the minimum in $D_{[s_{m},e_{m}]}$. As there are no other candidate
test levels $\tilde{f}_{[s_{m},e_{m}]}^{\\{j\\}}$ to choose from when
minimising
$\|\text{sign}(Y_{s_{m}:e_{m}}-\tilde{f}_{[s_{m},e_{m}]}^{\\{j\\}})\|_{\mathcal{I}^{a}_{[s_{m},e_{m}]}}$,
this results in $D_{[s_{m},e_{m}]}$ being relatively large, and therefore RNSP
keeps its 90% coverage promise for this model exactly, rather than generously
as in the other models. H-SMUCE behaves correctly for the two Gaussian models,
but fails for the discrete distributions and model Mix 1, which contains mass
points. It is unexpectedly successful in the Plain Cauchy model, but this is
perhaps it has very limited detection power in the Cauchy model with change-
points (more on this model below). It also underperforms for model Mix 2,
which is continuous (and within the domain of attraction of the Gaussian
distribution) but multimodal. SN-NSP fails for the discrete distributions,
which is a consequence of the (asymptotically guaranteed) closeness of the
self-normalised deviation measure to the appropriate functional of the Wiener
process not “kicking in” in these instances (due to the relatively small
sample sizes).

We now discuss performance for signals with change-points ($N>0$). Table 3
defines our models. Theorem 3.1 promises that any intervals of significance
returned by RNSP at levels $\alpha$ are such that, with probability at least
$1-\alpha$, they each contain at least one true change-point. In addition to
RNSP, H-SMUCE and SN-NSP, we also test versions of RNSP and SN-NSP with the
following overlap functions:

$\displaystyle\tau_{L}(\tilde{s},\tilde{e})$ $\displaystyle=$
$\displaystyle\lfloor(\tilde{s}+\tilde{e})/2\rfloor-\tilde{s},$
$\displaystyle\tau_{R}(\tilde{s},\tilde{e})$ $\displaystyle=$
$\displaystyle\lfloor(\tilde{s}+\tilde{e})/2\rfloor+1-\tilde{e}.$ (15)

This setting means that upon detecting a generic interval of significance
$[\tilde{s},\tilde{e}]$ within $[s,e]$, the RNSP and SN-NSP algorithms
continue on the left interval $[s,\lfloor(\tilde{s}+\tilde{e})/2\rfloor]$ and
the right interval $[\lfloor(\tilde{s}+\tilde{e})/2\rfloor+1,e]$ (recall that
the no-overlap case results uses the left interval $[s,\tilde{s}]$ and the
right interval $[\tilde{e},e]$). We denote the versions of the two procedures
with the overlaps as above by RNSP-O and SN-NSP-O, respectively. As before, we
set $\alpha=0.1$ for all methods tested.

For each model and method tested, we wish to evaluate the following aspects:
the empirical coverage (i.e. whether at least $(1-\alpha)100\%$ of the
simulated sample paths are such that any intervals of significance returned
contain at least true change-point each); if any intervals are returned, the
proportion of those that are genuine (i.e. the proportion of those intervals
returned that contain at least one true change-point); the number of genuine
intervals (i.e. the number of those intervals returned that contain at least
one true change-point); and the average length of genuine intervals (i.e. the
average length of those intervals returned that contain at least one true
change-point).

model name | sample path execution in R
---|---
Blocks |  `blocks + 10*rnorm(2048)`; `blocks` defined in Fryzlewicz (2014)
Cauchy | `c(rcauchy(100, 1), rcauchy(100, 2), rcauchy(100, 1))`
Bursts | `(c(rep(1, 200), rep(3, 80), rep(1, 200), rep(3, 80), rep(1, 200),`
| `rep(4, 40)) * rnorm(800))^2`
Poisson | `as.numeric(rpois(350, c(rep(1, 50), rep(4, 50), rep(10, 50), rep(2, 200)))`
Table 3: Non-null models for the comparative simulation study in Section 4. model | attribute | RNSP | RNSP-O | H-SMUCE | SN-NSP | SN-NSP-O
---|---|---|---|---|---|---
| coverage | 100 | 100 | 26 | 100 | 100
Blocks | prop. genuine interv. | 1 | 1 | 0.81 | 1 | 1
| no. of genuine interv. | 6.15 | 8.17 | 4.84 | 4.93 | 6.12
| avg. genuine interv. length | 90.86 | 99.92 | 94.729 | 160.19 | 182.19
| coverage | 100 | 100 | 100 | 100 | 100
Cauchy | prop. genuine interv. | 1 | 1 | N/A | 1 | 1
| no. of genuine interv. | 0.62 | 0.79 | 0 | 0.2 | 0.23
| avg. genuine interv. length | 119.62 | 123.94 | N/A | 153.8 | 154.55
| coverage | 100 | 100 | 39 | 100 | 100
Bursts | prop. genuine interv. | 1 | 1 | 0.78 | 1 | 1
| no. of genuine interv. | 2.99 | 4.36 | 3.1 | 4.27 | 5.26
| avg. genuine interv. length | 95.96 | 108.36 | 95.44 | 111.92 | 109.52
| coverage | 100 | 100 | 0 | 0 | 0
Poisson | prop. genuine interv. | 1 | 1 | 0.09 | 0.06 | 0.05
| no. of genuine interv. | 2.97 | 3 | 1.4 | 0.79 | 0.88
| avg. genuine interv. length | 37.12 | 37 | 20.86 | 31.35 | 28.95
Table 4: Results for each model+method combination: “coverage” is the number
of times, out of 100 simulated sample paths, that the respective model+method
combination did not return a spurious interval of significance; “prop. genuine
interv.” is the average (over 100 simulated sample paths) proportion of
genuine intervals out of all intervals returned, if any (if none are returned,
the corresponding 0/0 ratio is ignored in the average; N/A means there were no
detections in any of the 100 sample paths); “no. of genuine intervals” is the
average (over 100 sample paths) number of genuine intervals returned; “avg.
genuine interv. length” is the average (over 100 sample paths) length of a
genuine interval returned in the respective model+method combination.

Table 4 shows the results. H-SMUCE does not perform well in any scenario, not
even in the Blocks model, an instance of the homogeneous Gaussian model, a
simple sub-class of the heterogeneous Gaussian model class for which it was
specifically designed (where it achieves the coverage of 26, well short of the
expected 90). Its coverage of 100 in the Cauchy model is an artefact of the
fact that it does not achieve any detections over the 100 simulated sample
paths (so there are also no spurious detections).

RNSP and RNSP-O significantly outperform SN-NSP and SN-NSP-O in three out of
the four scenarios tested: Blocks, Cauchy and Poisson. In Blocks and Cauchy,
the RNSP methods achieve more detections and shorter intervals of significance
(so better localisation). In Poisson, in addition, they achieve much better
coverage (the SN-NSP methods are misled by the discrete nature of this
relatively low-intensity Poisson dataset, for which their required asymptotics
do not kick in, which results in a very large number of spurious detections).
However, the SN-NSP methods work better for the Bursts data in the sense that
they lead to more detections. The underlying reason is that the signal level
in this model is linearly proportional to the standard deviation of the noise,
which particularly suits the self-normalised SN-NSP methods. Note that the
“no. of genuine interv.” value of 5.26 for SN-NSP-O in Bursts, which is
greater than the true number of change-points (5), should not surprise here
despite the value of “coverage” being 100, because SN-NSP-O is a method that
can produce overlapping intervals, so two different intervals of significance
can on occasion correspond to the same change-point.

##  5 Data examples

###  5.1 The US ex-post real interest rate

Figure 1: The US ex-post real interest rate time series (black); intervals of
significance returned by RNSP (shaded red); their mid-points (blue). See
Section 5.1 for a detailed description.

We re-analyse the time series of US ex-post real interest rate (the three-
month treasury bill rate deflated by the CPI inflation rate) considered in
Garcia and Perron (1996), Bai and Perron (2003) and Fryzlewicz (2021). The
dataset is available at http://qed.econ.queensu.ca/jae/datasets/bai001/. The
time series, shown in Figure 1, is quarterly and the range is 1961:1–1986:3,
so $t=1,\ldots,T=103$.

RNSP appears to be an appropriate tool here, as the data displays
heterogeneity and possibly some heavy-tailed movements towards the latter
part. We run the RNSP algorithm with the default setting of $M=1000$, with
$\alpha=0.1$ and with overlaps as defined in (15). The procedure returns two
intervals of significance: $[23,75]$ and $[65,91]$. These are shown in Figure
1, together with their mid-points, which can serve as estimates of the change-
point locations. As with any RNSP execution with non-zero overlaps, one
question that may be asked is whether the two intervals may be indicating the
same change-point, but this, visually, is extremely unlikely here (the reason
for using non-zero overlaps is simply to provide larger samples for the RNSP
following the detection of the first interval; using zero overlaps means the
samples are too short and RNSP with zero overlaps does not pick up the second
change-point). Therefore, the solution points to a model with at least two
change-points. This is consistent (or at least not inconsistent) with both
Garcia and Perron (1996) (who also settle on a model with two change-points)
and Bai and Perron (2003) (who prefer a three-change-point model, not excluded
by RNSP here).

The difference between those two earlier analyses and ours is that those two
(a) were based on asymptotic arguments (and therefore valid asymptotically,
for unspecified large samples) and (b) were conditional in the sense that the
confidence regions for change-point locations in those two works were
conditional on the detection event. By contrast, our analysis via RNSP has a
finite-sample nature and the intervals of significance have an unconditional
character. Importantly, we do not make any distributional assumptions besides
independence and sign-symmetry, both of which are likely to be acceptable for
this dataset. The analysis via RNSP is unaffected by the likely heterogeneity
in the data.

###  5.2 Interest in the search term “data science”

We analyse the weekly interest in the search term “data science” from Google
Trends, in the US state of California. The link to obtain the data (accessed
on 24 August 2021) was
https://trends.google.com/trends/explore?date=today%205-y&geo=US-
CA&q=data%20science. Google Trends describe the data as follows. “Numbers
represent search interest relative to the highest point on the chart for the
given region and time. A value of 100 is the peak popularity for the term. A
value of 50 means that the term is half as popular. A score of 0 means there
was not enough data for this term.” Weeks in this data series start on Sundays
and the dataset spans the weeks from w/c 28th August 2016 to w/c 15th August
2021 (so almost five years’ worth of data). The observations are discrete
(integers from 22 to 100), which would likely pose difficulties for the
competing methods as outlined earlier.

Figure 2: Weekly relative interest (top=100) in the search term “data science”
in California in weeks from w/c 28 August 2016 to 15 August 2021 (black);
intervals of significance returned by RNSP (shaded red); their mid-points
(blue). See Section 5.2 for a detailed description.

We execute the RNSP procedure with the default setting of $M=1000$ and with
$\alpha=0.1$, with no overlaps, which returns the three intervals of
significance shown in Figure 2. The intervals are: w/c 23 April 2017 – w/c 10
June 2018 (interval 1), w/c 24 June 2018 – w/c 17 November 2019 (interval 2),
w/c 3 May 2020 – w/c 14 March 2021 (interval 3). While intervals 1 and 2
correspond to likely increases in the median interest, interval 3 clearly
corresponds to a likely decrease, and the mid-point of interval 3 (w/c 4
October 2020) visually aligns well with what seems to be a rather sudden drop
of interest.

It is difficult to speculate as to possible reasons for this perhaps
surprising drop of interest. A blog post on the popular website
https://towardsdatascience.com (https://bit.ly/3kqwDWO) reports that “2020 was
the first year since 2016, Data Scientist was not the number one job in
America, according to Glassdoor’s annual ranking. That title would belong to
Front End Engineer, followed by Java Developer, followed by Data Scientist.”
However, visually, a similar decline in interest is observed e.g. in the
analogous Google Trends series for the term “Java” (not shown). To offer some
comparative perspective, there is also a slight decline in the search
intensity for the term “artificial intelligence” (also picked up by RNSP; not
shown), but it is much less pronounced than those in “data science” or “Java”.
The role of the Covid-19 pandemic in these declines, their timings or relative
magnitudes is not immediately clear.

##  6 Discussion

Note that Corollary 3.1, crucial to the success of RNSP, can be rewritten as

$D_{[s_{m},e_{m}]}\leq\|\mathrm{sign}\\{X_{s_{m}:e_{m}}-Q_{1/2}(X_{s_{m}:e_{m}})\\}\|_{\mathcal{I}^{a}_{[s_{m},e_{m}]}},$
(16)

where $Q_{q}(\cdot)$ is the population $q$-quantile. However, the left-hand
side of (16) was not specifically constructed with $q=1/2$ in mind and
therefore the inequality

$D_{[s_{m},e_{m}]}\leq\|\mathrm{sign}\\{X_{s_{m}:e_{m}}-Q_{q}(X_{s_{m}:e_{m}})\\}\|_{\mathcal{I}^{a}_{[s_{m},e_{m}]}}$
(17)

is true for any $q\in(0,1)$. This shows that RNSP can equally be used for
significant change detection in any quantile, and not just the median.
However, for $q\neq 1/2$, the challenge is to obtain the null distribution of
$\|\mathrm{sign}\\{X-Q_{q}(X)\\}\|_{\mathcal{I}^{a}}$. Even if this challenge
is overcome (e.g. by simulation), RNSP as defined in this work may not be
effective for change detection in quantiles “far” from the median, due to the
particular way in which $D_{[s_{m},e_{m}]}$ is constructed (involving
minimisation over all levels). For RNSP to be a successful device for change
detection in other quantiles, the definition of $D_{[s_{m},e_{m}]}$ would have
to be modified to only minimise over ‘realistic’ candidate levels not far from
the population $q$-quantile under the local null.

## Acknowledgements

I wish to thank Shakeel Gavioli-Akilagun and Zakhar Kabluchko for helpful
discussions.

## References

  * Arlot and Celisse (2011) S. Arlot and A. Celisse.  Segmentation of the mean of heteroscedastic data via crossvalidation.  _Statist. Comput._ , 21:613–632, 2011.
  * Bai and Perron (1998) J. Bai and P. Perron.  Estimating and testing linear models with multiple structural changes.  _Econometrica_ , 66:47–78, 1998.
  * Bai and Perron (2003) J. Bai and P. Perron.  Computation and analysis of multiple structural change models.  _Journal of Applied Econometrics_ , 18:1–22, 2003.
  * Cheng et al. (2020) D. Cheng, Z. He, and A. Schwartzman.  Multiple testing of local extrema for detection of change points.  _Electron. J. Statist._ , 14:3705–3729, 2020.
  * Cho and Kirch (2021) H. Cho and C. Kirch.  Bootstrap confidence intervals for multiple change points based on moving sum procedures.  _Preprint_ , 2021.
  * Duy et al. (2020) V.N.L. Duy, H. Toda, R. Sugiyama, and I. Takeuchi.  Computing valid $p$-value for optimal changepoint by selective inference using dynaming programming.  In _Advances in Neural Information Processing Systems_ , volume 33, pages 11356–11367, 2020.
  * Eichinger and Kirch (2018) B. Eichinger and C. Kirch.  A MOSUM procedure for the estimation of multiple random change points.  _Bernoulli_ , 24:526–564, 2018.
  * Ellenberger et al. (2021) D. Ellenberger, B. Lausen, and T. Friede.  Exact change point detection with improved power in small-sample binomial sequences.  _Biometrical Journal_ , 63:558–574, 2021.
  * Fearnhead (2006) P. Fearnhead.  Exact and efficient Bayesian inference for multiple changepoint problems.  _Statistics and Computing_ , 16:203–213, 2006.
  * Frick et al. (2014) K. Frick, A. Munk, and H. Sieling.  Multiscale change-point inference (with discussion).  _Journal of the Royal Statistical Society Series B_ , 76:495–580, 2014.
  * Fryzlewicz (2014) P. Fryzlewicz.  Wild Binary Segmentation for multiple change-point detection.  _Ann. Stat._ , 42:2243–2281, 2014.
  * Fryzlewicz (2021) P. Fryzlewicz.  Narrowest Significance Pursuit: inference for multiple change-points in linear models.  _Preprint_ , 2021.
  * Garcia and Perron (1996) R. Garcia and P. Perron.  An analysis of the real interest rate under regime shifts.  _Review of Economics and Statistics_ , 78:111–125, 1996\.
  * Hao et al. (2013) N. Hao, Y. Niu, and H. Zhang.  Multiple change-point detection via a screening and ranking algorithm.  _Statistica Sinica_ , 23:1553–1572, 2013.
  * Hyun et al. (2018) S. Hyun, M. G’Sell, and R. Tibshirani.  Exact post-selection inference for the generalized lasso path.  _Electronic Journal of Statistics_ , 12:1053–1097, 2018\.
  * Hyun et al. (2021) S. Hyun, K. Lin, M. G’Sell, and R. Tibshirani.  Post-selection inference for changepoint detection algorithms with application to copy number variation data.  _Biometrics, to appear_ , 2021.
  * Jewell et al. (2020) S. Jewell, P. Fearnhead, and D. Witten.  Testing for a change in mean after changepoint detection.  _Preprint_ , 2020.
  * Kabluchko and Wang (2014) Z. Kabluchko and Y. Wang.  Limiting distribution for the maximal standardized increment of a random walk.  _Stoch. Proc. Appl._ , 124:2824–2867, 2014.
  * Kim et al. (2000) H.-J. Kim, M. Fay, E. Feuer, and D. Midthune.  Permutation tests for joinpoint regression with application to cancer rates.  _Statistics in Medicine_ , 19:335–351, 2000.
  * Lavielle and Lebarbier (2001) M. Lavielle and E. Lebarbier.  An application of MCMC methods for the multiple change-points problem.  _Signal Processing_ , 81:39–53, 2001.
  * Li (2016) H. Li.  _Variational Estimators in Statistical Multiscale Analysis_.  PhD thesis, Georg August University of Göttingen, 2016.
  * Li and Munk (2016) H. Li and A. Munk.  FDR-control in multiscale change-point segmentation.  _Electronic Journal of Statistics_ , 10:918–959, 2016.
  * Muggeo and Adelfio (2011) V. Muggeo and G. Adelfio.  Efficient change point detection for genomic sequences of continuous measurements.  _Bioinformatics_ , 27:161–166, 2011.
  * Nemirovski (1986) A. Nemirovski.  Nonparametric estimation of smooth regression functions.  _J. Comput. System Sci._ , 23:1–11, 1986.
  * Oomen (2019) R. Oomen.  Price signatures.  _Quantitative Finance_ , 19:733–761, 2019.
  * Pein et al. (2017) F. Pein, H. Sieling, and A. Munk.  Heterogeneous change point inference.  _J. Royal Stat. Soc. B_ , 79:1207–1227, 2017.
  * Rac̆kauskas and Suquet (2003) A. Rac̆kauskas and C. Suquet.  Invariance principle under self-normalization for nonidentically distributed random variables.  _Acta Applicandae Mathematicae_ , 79:83–103, 2003.