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

We propose Robust Narrowest Significance Pursuit (RNSP), a methodology for detecting localized 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 fulfils its coverage promises 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.


Introduction
The problem of uncertainty quantification for possibly multiple parameter changes in timeor space-ordered data is motivated by the practical question of whether any suspected changes reflect real structural changes in the underlying stochastic model, or are due to random fluctuations.Approaches to this problem include confidence sets associated with simultaneous multiscale change-point estimation (Frick et al., 2014;Pein et al., 2017;Vanegas et al., 2022), post-selection inference (Hyun et al., 2018(Hyun et al., , 2021;;Jewell et al., 2022;Duy et al., 2020), inference without selection and post-inference selection via Narrowest Significance Pursuit (Fryzlewicz, 2024), asymptotic confidence intervals conditional on the estimated change-point locations (Bai and Perron, 1998;Eichinger and Kirch, 2018;Cho and Kirch, 2022), 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.
In this paper, we are concerned with the following problem: given a sequence of noisy observations, automatically determine localized 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.
By its algorithmic construction, RNSP offers exact finite-sample coverage guarantees with practically no distributional assumptions on the data, other than serial independence of the signs of the true residuals and their sign-symmetry, a weak requirement which is immaterial for continuous distributions (as all, even non-symmetric, continuous distributions are signsymmetric).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.The execution of RNSP does not rely on having an estimate of the number or locations of change-points.Critical values needed by RNSP do not depend on the noise distribution, and can be accurately approximated analytically.RNSP explicitly targets the shortest possible intervals of significance.It is worth noting, however, that our large-sample consistency result for RNSP, shown in Section 4, relies on stronger assumptions than its finite-sample coverage properties, discussed in Sections 2 and 3. We now situate RNSP in the context of the existing literature.
Heterogeneous Simultaneous Multiscale Change Point Estimator, abbreviated as H-SMUCE (Pein et al., 2017), an extension of SMUCE (Frick et al., 2014), is a change-point detector in the heterogeneous Gaussian piecewise-constant model Y i = f i +σ i ε i , where f i is a piecewiseconstant signal, ε i are i.i.d.N (0, 1) variables, and σ i can only change when f i does.In Section A.1 of their work, the authors provide an algorithmic construction of confidence intervals for the locations of the change-points in f i , which 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 does not extend automatically to unknown noise distributions (as the analyst needs to know what distribution to sample from).In Section 5, 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 − α)% 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).Multiscale Quantile Segmentation (Vanegas et al., 2022, MQS), another extension of SMUCE, is a procedure for detecting possibly multiple changes in a given piecewise-constant quantile of the input data sequence, which includes the median as a special case.MQS estimates the number of change-points in the quantile function as the minimum among those candidate fits for which the empirical residuals pass a certain multiscale test at level α, where the empirical residuals are defined as binary exceedance sequences of the data over the level defined by each candidate fit.Working with such binary exceedance sequences means that MQS makes no distributional assumptions on the data other than their serial independence.
Like SMUCE and H-SMUCE, MQS then defines a confidence set around the estimated signals as a set of all feasible (in the same sense) signal fits at level α.This then enables the conceptual and algorithmic construction of asymptotic simultaneous confidence intervals for the change-point locations, which are guaranteed to (each) contain a change-point with probability 1−α+o(1).Chen et al. (2014) provide a critique of SMUCE from the inferential point of view, which also applies to H-SMUCE and MQS.In Section 5, we illustrate the advantages and disadvantages of MQS, as implemented in the R package mqs.In contrast to MQS, the coverage guarantees offered by RNSP are wholly based on exact inequalities and therefore hold for any sample size, which is reflected in its performance shown in Section 5.
The (non-robust) Narrowest Significance Pursuit (NSP; Fryzlewicz, 2024) 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, is within the domain of attraction of the normal distribution and is symmetric, neither of which is assumed in RNSP.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 overestimates 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 in Section 5 and show that RNSP offers a significant improvement.Bai andPerron (1998, 2003), working with least-squares estimation of multiple changepoints in regression models under possible heterogeneity of the errors, describe a procedure for computing confidence intervals conditional on detection, with asymptotic validity guarantees.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 6.1).By contrast, RNSP produces intervals of significant change in the median that are not conditional on detection and have a finite-sample nature.
The paper is organised as follows.Section 2 motivates RNSP and sets out its general al- 2 Motivation and review of RNSP algorithmic setting 2.1 RNSP: context and modus operandi RNSP discovers regions in the data in which the median departs from constancy, at a certain global significance level.This is in contrast to NSP (Fryzlewicz, 2024), which targets the mean.RNSP does not make moment assumptions about the data, and therefore the median is a natural measure of data centrality.We now 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 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 in which the signal {f t } T t=1 and the variables {Z t } T t=1 satisfy the assumptions below.Define sign(x) = I(x > 0) − I(x < 0), where I(•) is the indicator function.
Assumption 2.1 In (1), f t is a piecewise-constant vector with an unknown number N RNSP achieves the high level of generality in terms of the permitted distributions of Z t thanks to its use of the sign transformation.The use of 0 in sign(x) is critical for RNSP's objective to provide exact finite-sample coverage guarantees also for discrete distributions and continuous distributions with mass points, an aspect we discuss in Section 3.1.1.The sign function is a critical building block of various procedures for nonparametric changepoint testing and estimation; key references include Bhattacharyya and Johnson (1968), Carlstein (1988) and Dümbgen (1991).

(R)NSP: generic algorithm
The generic algorithmic framework underlying both RNSP and the non-robust NSP in Fryzlewicz (2024) is the same and is based on recursively searching for the shortest subsamples 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 broad class of distributions specified in Assumption 2.2.
We start with a pseudocode definition of the RNSP algorithm, in the form of a recursively defined function RNSP.
The RNSP algorithm is launched by the pair of calls: On completion, the output of RNSP is in the set S; when the context requires it, we write the output as S{RNSP(1, T, Y, M, λ α , τ L , τ R )}.We now comment on the RNSP function line by line.In lines 2-4, execution is terminated for intervals that are too short.A key ingredient of our measure of deviation is a multiresolution sup-norm introduced below, used on the signs of the input rather than in the original data domain.Its basic building block is a scaled partial sum statistic, defined for an arbitrary input sequence {x t } T t=1 by U s,e (x) = (e − s + 1) −1/2 e t=s x t .We define the multiresolution sup-norm (Nemirovski, 1986;Li, 2016) of an input vector x (of length T ) with respect to the interval set I as The set I used in RNSP contains intervals at a range of scales and locations.A canonical example of a suitable interval set I is the set I a of all subintervals of [1, T ].We will use 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 later).This will not affect the exactness of our coverage guarantees, because, naturally, if J ⊆ I, then ∥x∥ J ≤ ∥x∥ I .We also define the restriction of I to an arbitrary interval [s, e] as Note the trivial inequality When the above multiresolution sup-norm is applied to the signs of the input, as is done in this work, rather than the original input, we refer to it as the sign-multiresolution sup-norm.When applied to the empirical residuals from a candidate constant fit on an interval, it can be viewed as a simple robust multiscale measure of data fidelity of the given candidate fit on all time scales up to the length of the interval.
We now define the deviation measure D [sm,em] := 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 The Take candidate constants We have the following simple result; the proof is trivial and we omit it.
We now define our measure of deviation D [sm,em] , and prove its key property as a corollary to Proposition 3.1.
In other words, the deviation measure defined in (6) satisfies the desired property (3).
Proof.Let the index j 0 be as in the statement of Proposition 3.1.We have

□
This leads to the following guarantee for the RNSP algorithm.
Theorem 3.1 Let Assumption 2.1 hold, and let S = {S 1 , . . ., S R } be the set of intervals re-turned by the RNSP algorithm.We have Proof.On the set ∥sign(Z)∥ I a ≤ λ α , each interval S i must contain a change-point as if it did not, then by Corollary 3.1 and inequality (2), we would have to have However, the fact that S i was returned by RNSP means, by line 14 of the RNSP algorithm, that D S i > λ α , which contradicts (8).This completes the proof.□ Theorem 3.1 should be read as follows.Let α = P (∥sign(Z)∥ I a > λ α ).For a set of intervals returned by RNSP, we are guaranteed, with probability of at least 1 − α, that there is at least one change-point in each of these intervals.Therefore, S = {S 1 , . . ., 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 probability of obtaining one of more intervals of significance (R ≥ 1) is bounded from above by P (∥sign(Z)∥ I a > λ α ).Theorem 3.1 is of a finite-sample character and holds exactly and for any sample size.Moreover, it is independent of the form of the innovations Z. Assumptions on Z will only be needed in controlling the term P (∥sign(Z)∥ I a > λ α ); we defer this to Section 3.2.
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).This aspect of RNSP is investigated in Section 4. 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.
Therefore, one particular implication of Theorem 3.1 is that we must have 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 underlying signal f t , with a global probability of at least 1 − α.Therefore, as in NSP (Fryzlewicz, 2024), 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.

Deviation measure: discussion
Achieving computational savings without affecting coverage guarantees.The operation of trying each constant f {j} [sm,em] in ( 6) is fast, but in order to accelerate it further, we introduce the two computational savings below, which do not increase D [sm,em] and therefore respect the inequality (7) and hence also our coverage guarantees in Theorem 3.1.

Reducing the set I a
[sm,em] .To accelerate the computation of ( 6), we replace the set I a  [sm,em] has been defined in this particular way so as not to compromise the detection power in the piecewise-constant signal model.
To see this, consider the following illustrative example.Suppose Y t = f t (noiseless case) and change-point in f t on a generic interval [s m , e m ] under consideration, then in the noiseless case the multiresolution norm over the set I a [sm,em] is maximised at one of the "left" or "right" intervals in I lr [sm,em] , 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 and may therefore wish to ignore intervals for which e m − s m > L for a user-specified maximum length L.
Validity for non-continuously distributed data.The following three aspects of RNSP ensure the validity of its coverage guarantees in the presence of mass points in Z t or if the distribution of Z t is discrete: (a) the fact that the sign function defined in Section 2.1 returns zero if its argument is zero, (b) the fact that the test levels f {j} [sm,em] (definition ( 4)) are placed at data points, and (c) the fact that these test levels are placed in between the sorted data points.Indeed, in the absence of ( a

Evaluation and bounds for ∥sign(Z)∥ I a
To make Theorem 3.1 operational, we need to obtain an understanding of the distribution of ∥sign(Z)∥ I a so we are able to choose λ α such that P (∥sign(Z)∥ I a > λ α ) = α (or approximately so) for a desired global significance level α.
Initially we consider Z t such that P (sign(Z t ) = 1) = P (sign(Z t ) = −1) = 1/2 (the general case P (sign(Z t ) = 0) ≥ 0 is covered in the next paragraph).One simple way of determining the distribution of ∥sign(Z)∥ I a for any finite T is by simulation; this would only need to be done once for every T and the quantiles stored for fast access.Another approach is asymptotic and proceeds as follows.From Theorem 1.1 in Kabluchko and Wang (2014) (which applies to sequences of serially independent symmetric Bernoulli variables as explained in Section 1.5.1 of that work; our Assumption 2.2 means that result is applicable in our context), we have lim where a T = {2 log(T log −1/2 T )} 1/2 and Λ is a constant.As the theoretical calculation of Λ in Kabluchko and Wang (2014)  While this approach is asymptotic in nature (note the limit as T → ∞ in ( 10)), we observe that the finite-sample agreement of P (∥sign(Z)∥ I a > a T + τ /a T ) with its limit in ( 10) is excellent even for small sample sizes.If the user chooses to pursue this route of obtaining λ α , this will be the only asymptotic component of RNSP.
Suppose now that P (sign(Z t ) = 0) = ρ t ≥ 0; note that the sign-symmetry Assumption 2.2(b) implies P (sign holds because every constituent partial sum of ∥sign(Z)∥ I a has a corresponding larger or equal in magnitude partial sum in ∥sign( Z)∥ 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 sign(Z t ) starts is majorised by the absolute partial sum | − 1 + 1 + 1|/ √ 3, a constituent of ∥sign( Z)∥ 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 (11) holds simply because T 1 ≤ T .The implication of ( 11) is that ∥sign(Z)∥ I a for ρ t ≥ 0 is majorised by ∥sign(Z)∥ I a for ρ t = 0, the case handled by (10).Therefore, the threshold λ α obtained as a consequence of (10) can also meaningfully be applied in the general case ρ t ≥ 0.

Detection consistency and lengths of RNSP intervals
This section shows the large-sample consistency of RNSP in detecting change-points, and the rates at which the lengths of the RNSP intervals contract (relative to T ), as T increases.
To simplify our technical arguments, we consider a version of the RNSP algorithm that considers all subintervals of [1, T ].Our focus on i.i.d.Z t 's in this section is mainly due to our ability to rely on the Dvoretzky-Kiefer-Wolfowitz inequality in the proof of Corollary 4.1; however, note that the finite-sample result of Theorem 4.1 holds regardless of the dependence structure of Z t .We focus on continuously-distributed Z t 's as this results in notationally much less involved arguments regarding the minimum signal strength required.
We first introduce some essential notation, and then state our assumption and the result.
For each change-point η j , define In addition, for any process V t , define ϵ V t (w) = I(V t − w > 0) − P (V t − w > 0).executed with no overlaps and with threshold λ α .We have Note that here, λ α does not depend on α and is of a higher order of magnitude than the (α-dependent) λ α of Section 3.2.Paraphrasing, this is to say that we need the global significance level α to tend to zero with T in order to obtain large-sample consistency; a result in line with analogous results in Pein et al. (2017) and Vanegas et al. (2022).
We briefly comment on what the result of Corollary 4.1 means for the minimum signal strength Assumption 4.1(iii), and for the localisation rates of the RNSP algorithm in detecting the change-points.If the distribution of Z t does not vary with T (this section already assumes that it does not vary with t), and if the jump sizes |f η j − f η j−1 | are bounded from below by a positive constant independent of j and T , then ∆ j (formula ( 12)) is also bounded from below by a positive constant independent of j and T .By formula (13), the assumptions on λ, λ α in Corollary 4.1 then imply dj = Θ(log T ) (where Θ should be read "of the exact order").Assumption 4.1(iii) then requires that the spacings between the change-points be at least of order log T .Corollary 4.1 states that the length of each RNSP interval of significance, e j − s j + 1, is, with global probability approaching 1, at most of order log T .
These minimum-spacing assumptions and the implied lengths of the localisation intervals are near-optimal and the same as those in the non-robust literature, see e.g.Theorem 1 in Baranowski et al. (2019) and Corollary 4 in Fryzlewicz (2024), as well as the associated discussions.However, the results of this section also permit ∆ j → 0 with T , which will, naturally, affect the above minimum-spacing requirements and localisation rates as stipulated by formulae ( 12) and ( 13) and Assumption 4.1(iii).
The consistency of RNSP in the sense of Corollary 4.1 implies the consistency of any pointwise estimators ηj contained within the RNSP intervals of significance [s j , e j ], with the localisation rate of ηj bounded from above by the length of the interval [s j , e j ].In particular, the near-optimality of the lengths of the RNSP intervals [s j , e j ] (as discussed in the preceding paragraph) automatically implies the near-optimality of the localisation rate of ηj .Indeed, since by Corollary 4.1, [s j , e j −1] is guaranteed to contain η j (on an event of high probability), for any estimator ηj ∈ [s j , e j − 1], we must have |η j − η j | ≤ |e j − s j − 1| (on the same event).In other words, the rate with which [s j , e j ] contract (relative to T ) is inherited by any estimator ηj ∈ [s j , e j − 1]; this applies even to naive estimators constructed e.g. as the middle points of their corresponding RNSP intervals [s j , e j ], i.e. ηj = ⌊(s j + e j )/2⌋.
More refined estimators, e.g. one based on the CUSUM maximisation of the signs of the data around their median (Sen and Srivastava, 1975) within each RNSP interval, can also be used and will also automatically inherit the consistency and rate.Sections executed with no overlaps and with threshold λ α .We then have lim inf The computation of the deviation measure D  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 α returns no intervals of significance with probability at least 1−α.In this section, we use α = 0.1.There are analogous parameters in H-SMUCE, MQS and SN-NSP, and they are also set to 0.1.
However, while in both RNSP and SN-NSP, the parameter α is responsible for finite-sample coverage guarantees (for both null and non-null models), in MQS and H-SMUCE it is responsible for similar but only asymptotic coverage guarantees, as T → ∞.For H-SMUCE, this latter point is clarified (jointly) in Theorem 5 of Pein et al. (2017)  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).
This setting means that upon detecting a generic interval of significance [s, ẽ] within [s, e], the RNSP and SN-NSP algorithms continue on the left interval [s, ⌊(s + ẽ)/2⌋] and the right interval [⌊(s + ẽ)/2⌋ + 1, e] (recall that the no-overlap case results uses the left interval [s, s] and the right interval [ẽ, 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 α = 0.1 for all methods tested.
For each model and method tested, we evaluate the following metrics, which collectively promote the detection of genuine, and the non-detection of spurious, intervals of significance: [coverage] the empirical coverage (i.e.whether at least (1 − α)100% of the simulated sample paths are such that any intervals of significance returned contain at least one true change-point each); [prop.gen.int.] 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); [no.gen.int.] the number of genuine intervals (i.e. the number of those intervals returned that contain at least one true change-point); and [av.gen.int.len.] the average length of genuine intervals (i.e. the average length of those intervals returned that contain at least one true change-point).No single measure describes the performance of (any) method accurately, but the four measures in conjunction provide a clear picture.As a cartoon illustration of how naive solutions are able to skew some of these measures but not the others, consider a putative interval estimator that always returns the longest possible interval of [1, T ].For a model with ≥ 1 change-points, "coverage" and "prop.gen.int" will return 100 and 1, respectively.However, this naive solution will be penalised by the next two measures, "no.gen.int." and "av.gen.int.len.".For models with a single changepoint, the [1, T ] solution will return an interval that in most cases will be unreasonably long, and this will be picked up by the "av.gen.int.len." measure.In addition, for models with more than one change-point, "no.gen.int." (equal to one) will be inaccurate.
points as this offers the same amount of evidence on either side of the change-pointshence the frequent empirical closeness of RNSP interval midpoints to the truth.) 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 30, 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 almost any detections over the 100 simulated sample paths (so there are also no spurious detections).
With the exception of MQS.easy, the least challenging model, MQS frequently produces spurious intervals in signals with change-points: the coverage figures for MQS are well short of 90% in all but one models tested.On the upside, in MQS.easy (the only set-up in which MQS achieves the correct coverage), the average length of genuine intervals is around 20% below that of RNSP."coverage" is the percentage of times that the respective model+method combination did not return a spurious interval of significance; "prop.gen.int." is the average 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); "no.gen.int." is the average number of genuine intervals returned; "av.gen.int.len." is the average length of a genuine interval returned in the respective model+method combination; "MSE (•)" is the MSE of ft constructed as described in the text, with the respective pointwise change-point estimation in brackets. (* ) note HSMUCE uses its own pointwise change-point estimation.Significance level α = 0.1 (nominal coverage = 90%).All numbers rounded to two decimal digits except when such rounding takes a positive number to zero. 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 changepoint 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.It is difficult to speculate as to possible reasons for this 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).

Discussion
Other quantiles.Note that Corollary 3.1, crucial to the success of RNSP, can be rewritten as where Q q (•) is the population q-quantile.However, the left-hand side of (15) was not specifically constructed with q = 1/2 in mind and therefore the inequality ,em] is true for any q ∈ (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 ̸ = 1/2, the challenge is to obtain the null distribution of ∥sign{Y −Q q (Y )}∥ 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 [sm,em] is constructed (involving minimisation over all levels).For RNSP to be a successful device for change detection in other quantiles, the definition of D [sm,em] would have to be modified to only minimise over 'realistic' candidate levels not far from the population q-quantile under the local null.Continuing from ( 16), this implies The infimum over w will be achieved if both elements of the maximum are the same (or otherwise it would be possible to alter w slightly to decrease the larger of the two moduli).
But this is only possible if w ∈ (f η 1 +1 , f η 1 ) and hence, bearing in mind that Z t is median-zero, we have Let w 0 be such that P {Z t ∈ (w 0 − f η 1 , 0)} = P {Z t ∈ (0, w 0 − f η 1 +1 )}.Since (w 0 − f η 1 +1 ) − Therefore, Continuing from (17), we therefore have As RNSP looks for the shortest intervals of significance, the length of the RNSP interval will not exceed that of [η 1 − d + 1, η 1 + d], which is 2d.From (19), its length will therefore be bounded from above by 2 2λ+λα 2∆ 1 2 + 1 = 2 d1 .We now turn our attention to the multiple change-point case.Note that even though the RNSP interval of significance around η j is guaranteed to be of length at most 2 dj , it will not necessarily be a sub-interval of [η j − dj + 1, η j + dj ].Therefore in order that an interval detection around η j does not interfere with detections around η j−1 or η j+1 , the distances η j − η j−1 and η j+1 − η j must be suitably long, but this is guaranteed by Assumption 4.1(iii).This completes the proof.□Theorem 4.1 implies the result.□

Funding acknowledgement and disclosure statement
The author acknowledges EPSRC grant EP/V053639/1.There are no competing interests to declare.
gorithmic framework.Section 3 describes how RNSP measures the local deviation from model constancy and gives finite-sample theoretical performance guarantees for RNSP.Section 4 quantifies the large-sample behaviour of RNSP.Section 5 contains numerical examples and comparisons.Section 6 includes examples showing the practical usefulness of RNSP.Section 7 concludes with a brief discussion.Software implementing RNSP is available at https://github.com/pfryz/nsp.

s m 0
, e m 0 in place of s, e.The chosen interval is denoted by[s, ẽ].This second-stage search is important to RNSP's pursuit to produce short intervals: indeed, if the sample of intervals [s m , e m ] contained insufficiently short intervals (perhaps because an insufficiently large M was chosen), then, without the second-stage search in line 19, the intervals of significance returned by RNSP might be overly long.The second-stage search in line 19 can be seen as a guard against a small M , or in other words against an insufficiently fine original grid of interval endpoints.In line 20, the selected interval [s, ẽ] is added to the output set S.Because of the second-stage search in line 19, pre-drawing the intervals [s m , e m ] prior to launching the RNSP procedure (rather than drawing them recursively on each current interval as it done in the algorithm) is not an option for RNSP.Indeed, in the presence of the second-stage search, the selected interval of significance [s, ẽ] may be misaligned with the initial grid of intervals drawn, in which case the grid of intervals to the left and to the right of [s, ẽ] must be re-drawn to avoid leaving un-examined gaps in the data.In lines 21-22, RNSP is executed recursively to the left and to the right of the detected interval[s, ẽ].However, we optionally allow for some overlap with[s, ẽ].The overlap, if present, is a function of[s, ẽ]  and, if it involves detection of the location of a change-point within [s, ẽ], then it is also a function of Y .An example of the relevance of this is given in Section 6.1.3Robust NSP: measuring deviation from the constant model 3.1 Deviation measure: motivation, definition and properties The main structure of the DeviationFromConstantModel(s m , e m , Y ) operation is as follows: (1) Fit the best, in the sense described precisely later, constant model to Y sm:em .(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.
discussion below assumes that there is no change-point in [s m , e m ].For the true constant signal f sm:em , denote f [sm,em] := f sm = . . .= f em .There are only at most 2(e m −s m )+3 different possible constants f[sm,em] leading to different sequences {sign(Y t − f[sm,em] )} em t=sm .To see this, sort the values of Y sm:em in non-decreasing order to create Y (1) , Y (2) , . . ., Y (em−sm+1) .
contains an error, we use simulation over a range of values of T and τ to determine a suitable value of Λ as 0.274.The practical choice of the significance threshold λ α then proceed as follows: (a) fix α to the desired level, for example 0.05 or 0.1; (b) obtain the value of τ by equating 1−exp(−2Λ exp(−τ )) = α; (c) obtain λ α = a T +τ /a T .
5 and 6 illustrate both of these estimators.Trivially, in light of Corollary 4.1, the set of estimated change-points {η j } R j=1 is consistent in the Hausdorff measure for {η j } N j=1 , the set of true change-points in f t .Yet another implication of Theorem 4.1 appears below.Corollary 4.2 Let the assumptions of Corollary 4.1 hold.Let λ α be such that P (∥sign(Z)∥ I a ≤ λ α ) ≥ 1 − α and let λ = (1 + δ) log 1/2 T , for any δ > 0. Let S denote the set of intervals of significance [s 1 , e 1 ] < . . .< [s R , e R ] returned by RNSP algorithm that considers all intervals, [s,e]  for all subintervals [s, e] of [1, T ], as the results of this section require, is a O(T 3 ) operation.The cubic complexity should not surprise in the context of a robust method that considers all intervals, as there are O(T 2 ) intervals to consider and O(T ) binary exceedance levels within each interval.In practice, we use three devices to reduce the computational complexity of RNSP: (a) using a fixed (i.e.unchanging with T ) value of M , (b) restricting the length of intervals under consideration to ≤ L, and (c) reducing the set I a [s,e] .Items (b) and (c) are described in more detail in Section 3.1.1.

Figure 1 :
Figure 1: Left: The US 3-month ex-post real interest rate time series (black); intervals of significance returned by RNSP (transparent pink); their midpoints (red); argument-maxima of absolute CUSUMs of signs of the data around the median in each interval of significance (blue).See Section 6.1 for a detailed description.Right: The US 1-month real interest rate time series (black); intervals of significance returned by RNSP (transparent pink); their midpoints (red); argument-maxima of absolute CUSUMs of signs of the data around the median in each interval of significance (blue).Superimposed on the bottom of the graph is the probability of recession time series (purple), on the scale of 0 (bottom purple dashed line) to 1 (top purple dashed line).RNSP significance level α = 0.1.See Section 6.1 for a detailed description.

6. 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 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.We execute the RNSP procedure with the default setting of M = 1000 and with α = 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 2018w/c 17 November 2019 (interval 2), w/c 3 May 2020 -w/c 14 March 2021 (interval 3).

Figure 2 :
Figure2: 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 (transparent pink); their midpoints (red); argument-maxima of absolute CUSUMs of signs of the data around the median in each interval of significance (blue).RNSP significance level α = 0.1.See Section 6.2 for a detailed description.
Set of feasible signals.It is convenient to define the set of feasible signals f 0 t at level α with respect to the algorithmic execution RNSP(1, T, •, M, λ α , τ L , τ R ) by F α = {f 0 :S{RNSP(1, T, Y − f 0 , M, λ α , τ L , τ R )} = ∅}, where λ α is such that P (∥sign(Z)∥ I a ≤ λ α ) ≥ 1 − α.That is, F α is the set of postulated signals f 0 such that, on fitting f 0 to the data and obtaining the empirical residuals, RNSP cannot distinguish (at level α) the empirical residuals from white noise.For any candidate fit f 0 , it is straightforward to check whether or not it is a member of F α by executing RNSP(1, T, Y − f 0 , M, λ α , τ L , τ R ).A ProofsProof of Theorem 4.1.Consider initially the case of a single change-point η 1 .RNSP will, among others, consider intervals symmetric about the true change-point, i.e. [η 1 −d+1, η 1 + d], for all appropriate d.Take a constant candidate fit w on the interval [η 1 − d + 1, η 1 + d] and define U t (w) := sign(Y t − w) = 2I(Y t − w > 0) − 1 (the latter equality holds due to the continuity of the distribution of Z t ).Assume wlog f η 1 > f η 1 +1 .We have D [η 1 −d+1,η 1 +d]
{Z t } T t=1 do not have to be identically distributed, and can have arbitrary mass atoms, or none, as long as their distributions satisfy Assumption 2.2.The distribution(s) of {Z t } T where M is the median operator.(If the median is non-unique, we require 0 ∈ M (Z t ).)(b) The variables {Z t } T t=1 are sign-symmetric, i.e.P (Z t > 0) = P (Z t < 0), ∀ t. t=1 can be unknown to the analyst, and we do not impose moment assumptions.Any zeromedian continuous distribution (even one with an asymmetric density function) is signsymmetric.The requirement of the independence of {sign(Z t )} T t=1 is weaker than that of the independence of {Z t } T t=1 itself: e.g. if Z t is a (G)ARCH process driven by symmetric, independent innovations, then sign(Z t ) is serially independent, while Z t is not.While the results of Section 3 require the independence of sign(Z t ), those in Section 4 require the independence of Z t ; see Section 4 for details.
M 0 := arg min m {e m − s m : m = 1, . . ., M ; D [sm,em] > λ α } 1: function RNSP(s, e, Y , M , λ α , τ L , τ R ) 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 the shortest sub-interval on which the hypothesis of the baseline model is rejected locally at a global level α.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; it proceeds by executing lines 2-18 of RNSP, but with Definition 3.1 Let the constants f {j} [sm,em] , j = 1, . . ., 2(e m − s m ) + 3 be defined as in (4).
D[sm,em]tries all possible baseline constant model fits on [s m , e m ] and chooses the one that minimises the sign-multiresolution norm of the residuals, to ensure that the finite-sample coverage guarantees hold, as we will see below.Corollary 3.1 Under Assumption 2.1, for any interval [s m , e m ] on which there is no change-point, we have D [sm,em] ≤ ∥sign(Z sm:em )∥ I a [sm,em] .

Table 1 :
Null models for the comparative simulation study in Section 5.

Table 2 :
Percentage of times, out of 200 simulated sample paths of each null model, that the respective method indicates no intervals of significance at level α = 0.1 (nominal coverage = 90%).
Vanegas et al. (2022) the online supplement to that work, and for MQS -in Theorem 2.3 ofVanegas et al. (2022)and in Section S.4 of the online supplement to that work.The null models used are listed in Table1, and Table2shows the associated results.RNSP, MQS-R and MQS-K keep the nominal size well across all the models considered, returning no intervals of significance at least 95% of the time in all situations.H-SMUCE behaves correctly for the three 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 because it has very limited detection power in the Cauchy model with change-points (more on this model below).It also underperforms slightly for model Mix 2, which is continuous (and within the domain of attraction of the Gaussian distribution) Table 3 defines our models; the model labelled MQS.easy is taken from https://github.com/ljvanegas/ mqs/blob/master/mqs.ipynb and MQS.hard and MQS.vhard are its lower signal-to-noise versions.Theorem 3.1 promises that any intervals of significance returned by RNSP at levels α are such that, with probability at least 1 − α, they each contain at least one true change-point.In addition to RNSP, H-SMUCE, MQS-R, MQS-K and SN-NSP, we also test versions of RNSP and SN-NSP with the following overlap functions:

Table 4 :
Results for each model+method combination, out of 200 simulated sample paths: To illustrate RNSP on a more recent dataset of a similar nature, we examine the 1-month US real interest rate, available from https://fred.stlouisfed.org/series/REAINTRATREARAT1MO.The time series, shown in the right plot of Figure1, is monthly and runs from January 1982 to February 2023.We run RNSP with M = 1000, α = 0.1 and no overlaps.The intervals of significance returned by RNSP appear visually plausible and it is interesting (albeit not unexpected) to observe that the periods of peaking probabilities of recession (data available from: https://fred.stlouisfed.org/series/RECPROUSM156N)from year 2000 onwards are wholly contained within RNSP intervals of significance, indicating declines in the 1-month real interest rate.The period of high probability of recession in 1982 does not appear supported in the 1-month real rate data in a way detectable to RNSP.Finally, the period of high probability of recession in 1990 appears to coincide with a falling 1-month rate but RNSP has, in this case, a visually justifiable preference for intervals just before and just after this likely recession period.