Multiscale and multilevel technique for consistent segmentation of nonstationary time series
Haeran Cho
University of Bristol, Bristol, UK.

Piotr Fryzlewicz
London School of Economics, London, UK.
Summary. The importance of nonstationary modelling in time series analysis was borne out,for example, by the recent financial crisis, which, as some authors argue, has been partly
caused by over-reliance on stationary time series models in financial risk management. Piecewise stationarity is the simplest form of departure from stationarity, and one basic task under
this assumption is time series segmentation.
In this paper, we propose a fast, well-performing and consistent method for segmenting apiecewise-stationary, linear time series with an unknown number of breakpoints. The timeseries model we use is the nonparametric Locally Stationary Wavelet model, in which a complete description of the piecewise-stationary second-order structure is provided by wavelet
periodograms computed at multiple scales and locations. The initial stage of our method is a
new binary segmentation procedure, with a theoretically justified and rapidly computable test
criterion, which detects breakpoints in wavelet periodograms at each scale separately. It isthen followed by within-scale and across-scales post-processing steps, leading to consistentestimation of the number and locations of breakpoints in the second-order structure of theoriginal process.
An extensive simulation study demonstrates good performance of our method in comparison
to the state of the art, and its application to the Dow Jones index indicates two breakpoints,
each corresponding to a significant event in the recent financial crisis. A complete R script,
implementing our methodology, is provided.
Keywords: binary segmentation, breakpoint detection, locally stationary wavelet model, piecewise stationarity, post-processing, wavelet periodogram.

1. Introduction
Stationarity assumption is appealing when analysing short time series. However, it is often
unrealistic in many circumstances, for example when observing time series evolving in naturally nonstationary environments. One such example can be found in econometrics, where
price processes are considered to have time-varying variance in response to events taking
place in the market. For example, given the explosion of market volatility during the recent
financial crisis, it is unlikely that the same stationary time series model can accurately describe the evolution of market prices before and during the crisis. Indeed, Janeway (2009)
goes further and argues that the (over)use of stationary time series models in financial risk
management might have been a contributing factor in the crisis due to those models’ lackof flexibility in reacting to rapid changes in the statistical properties of the markets. As ataster, we note that some interesting findings resulting from the application of our methodology to financial data observed before and during the crisis are presented in Section 5 ofthis paper.

2

Haeran Cho and Piotr Fryzlewicz

Exact or approximate piecewise stationarity is a well studied and arguably the simplest
form of departure from stationarity, and one task of interest when faced with data of thisform is to detect their structural breakpoints so that each segment between two breakpoints
can be modelled as approximately stationary.
Breakpoint detection has received considerable attention and many methods have beendeveloped which can be broadly categorized into two: retrospective (a posteriori) methods
and on-line methods. The on-line approach is employed in areas such as quality control,
where one wishes to detect changes in the process whilst monitoring is still in progress. A
survey of literature on the on-line segmentation of serially independent observations can befound in Bhattacharya (1994). Methods for on-line segmentation of sequentially observed
time series include Ombao et al. (2004), who use approximation via autoregressive (AR)processes, and Choi et al. (2008), who analyse the changing spectral characteristics of theprocess.
On the other hand, the “a posteriori” approach takes into account the entire set ofobservations at once and detects breakpoints which occurred in the past. The resulting
segmentation can be of interest e.g. from the point of view of forecasting, where informationfrom the last (approximately) stationary segment can be useful in forecasting the future.
Our interest in this article lies in the “a posteriori” segmentation category, and we propose
a retrospective segmentation procedure which achieves consistency in identifying multiple
breakpoints for a class of nonstationary processes.
Early segmentation literature was mostly devoted to testing the existence of a single
breakpoint in the mean or variance of independent observations (Chernoff and Zacks, 1964;
Sen and Srivastava, 1975; Hawkins, 1977; Hsu, 1977; Worsley, 1986). When the presenceof more than one breakpoint is suspected, an algorithm for detecting multiple breakpoints
is needed. In Vostrikova (1981), a “binary segmentation” procedure was introduced, acomputationally efficient multilevel version of the CUSUM test, which recursively locates
and tests for multiple breakpoints, producing consistent breakpoint estimators for a class ofrandom processes with piecewise constant means. However, the critical values of the tests
at each stage were difficult to compute in practice due to the stochasticity in previously
selected breakpoints. Venkatraman (1993) employed the same procedure to find multiple
breakpoints in the mean of independent and normally distributed variables and showed theconsistency of the detected breakpoints with the tests depending on the sample size only,
and thus being easier to compute. The binary segmentation procedure was also adopted todetect multiple shifts in the variance of independent observations (Inclán and Tiao, 1994;
Chen and Gupta, 1997).
Various multiple breakpoint detection methods have been proposed for time series ofdependent observations. In Lavielle and Moulines (2000), least squares estimators of breakpoint locations were developed for linear processes, extending the work of Bai and Perron(1998). In the Bayesian framework, a number of procedures were studied for the segmentation of signals using piecewise constant linear regression models (McCulloch and Tsay, 1993;
Punskaya et al., 2002; Fearnhead, 2005). Adak (1998) and Ombao et al. (2001) proposed
methods which divided the time series into dyadic blocks and chose the best segmentation according to suitably tailored criteria. Whitcher et al. (2000, 2002) and Gabbanini
et al. (2004) suggested to segment long memory processes by applying the iterative cumulative sum of squares (ICSS) algorithm (proposed in Inclán and Tiao (1994)) to discrete
wavelet coefficients of the process, which were shown to be approximately Gaussian anddecorrelated. Davis et al. (2006) developed the Auto-PARM procedure which found the
optimal segmentation of piecewise stationary AR processes via the minimum description

Time series segmentation

3

length principle, later extended to the segmentation of non-linear processes in Davis et al.
(2008). Kokoszka and Leipus (2000) and Andreou and Ghysels (2002) studied the problem
of detecting structural shifts in the parameters of ARCH or GARCH models.
The aim of our work is to propose a well-performing, theoretically tractable and fast
procedure for segmenting piecewise stationary time series which are linear but otherwise do
not follow any particular parametric model. (We note that we use the term “segmentation”
interchangeably with “multiple breakpoint detection”.) The nonparametric model we use
for this purpose is the Locally Stationary Wavelet (LSW) model first proposed by Nasonet al. (2000) and later studied by Fryzlewicz and Nason (2006) and Van Bellegem and vonSachs (2008). Our choice of model is motivated by several factors, which we explain in
Section 2 below. In the LSW model, the piecewise constant second-order structure of theprocess is completely described by the local wavelet periodograms at multiple scales, and
it is those basic statistics that we use as a basis of our segmentation procedure.
To achieve the multiple breakpoint detection, we propose a binary segmentation method,which we apply to wavelet periodograms at each scale separately, and then propose a withinscale and across-scales post-processing procedure to obtain consistent estimators of breakpoints in the second-order structure of the process. We note that wavelet periodograms
follow a multiplicative statistical model, but our binary segmentation procedure is different
from previously proposed binary segmentation methods for multiplicative models (Inclán
and Tiao, 1994; Chen and Gupta, 1997) in that it allows correlated data, which is essential
when working with wavelet periodograms. We also mention other unique ingredients ofour breakpoint detection procedure which lead to its good performance and consistency inprobability; these are: the theoretical derivation of our test criterion (which only depends on
the length of the time series and is thus fast to compute); and the novel across-scales postprocessing step, essential in combining the results of the binary segmentation procedures
performed for each wavelet periodogram scale separately.
We note that our method can simultaneously be termed “multiscale” and “multilevel”, as
the basic time series model used for our purpose is a wavelet-based, and thus a “multiscale”
model; and the core methodology to segment each scale of the wavelet periodogram in themodel is based on binary segmentation and is thus a “multilevel” procedure.
The paper is organised as follows. Section 2 explains the LSW model and justifies its
choice. Our breakpoint detection methodology (together with the post-processing steps) is
introduced in Section 3, where we also demonstrate its theoretical consistency. In Section
4, we describe the outcome of an extensive simulation study which demonstrates good performance of our method in comparison with the state of the art. In Section 5, we apply ourtechnique to the segmentation of the Dow Jones index, which results in an exciting discoveryof two breakpoints: one coinciding with the initial period of the recent financial crisis, andthe other coinciding with the recent collapse of Lehman Brothers, a major financial services
firm. The proofs of our theoretical results are in the appendix.
Software (an R script) implementing our methodology is available from: http://www.
maths.bris.ac.uk/~ mahrc/msml_technique.html.
2. Locally stationary wavelet time series
In this section, we first define the Locally Stationary Wavelet (LSW) time series model(noting that our definition is a slight modification of Fryzlewicz and Nason (2006)), and thenjustify its choice as an attractive framework for the purpose of developing our methodology

4

Haeran Cho and Piotr Fryzlewicz

for time series segmentation.
−1
Definition 1. A triangular stochastic array {Xt,T }Tt=0
for T = 1, 2, . . . , is in a classof Locally Stationary Wavelet (LSW) processes if there exists a mean-square representation

Xt,T =−1
X∞
X

Wi (k/T )ψi,t−k ξi,k(1)

i=−∞ k=−∞where i ∈ {−1, −2, . . .} and k ∈ Z are scale and location parameters respectively, ψi =
(ψi,0 , . . . , ψi,Li ) are discrete, real-valued, compactly supported, non-decimated wavelet vectors, and ξi,k are zero-mean, orthonormal, identically distributed random variables. Alsofor each i ≤ −1, Wi (z) : [0, 1] → R is a real-valued, piecewise constant function with a finite(but unknown) number of jumps. Let Li denote the total magnitude of jumps in Wi2 (z).
The functions Wi (z) satisfyP−1
2•i=−∞ Wi (z) < ∞ uniformly in z,P−1−i•i=−I 2 Li = O(log T ) where I = log2 T .

The reader unfamiliar with basic concepts in wavelet analysis is referred, at this point,
to the excellent monograph by Vidakovic (1999). Throughout the paper, ξi,k are assumed
to follow the normal distribution; extensions to non-Gaussianity are possible but would
cause substantial technical difficulties in deriving our theoretical results. Comparing theabove definition with the Cramér’s representation of stationary processes, Wi (k/T ) is a(scale- and location-dependent) transfer function, the wavelet vectors ψi are analogous to
the Fourier exponentials, and the innovations ξi,k correspond to the orthonormal increment
process. Small negative values of the scale parameter i denote “fine” scales where the wavelet
vectors are the most localised and oscillatory; large negative values denote “coarser” scaleswith longer, less oscillatory wavelet vectors. By assuming that Wi (z) is piecewise constant,
we are able to model processes with a piecewise constant second-order structure where,
between any two breakpoints in Wi (z), the autocovariance function remains constant. The
Evolutionary Wavelet Spectrum (EWS) is defined as Si (z) = Wi (z)2 , and is in one-toone correspondence with the time-dependent autocovariance function of the process. Our
primary objective is to develop a consistent method for detecting breakpoints in the EWS,
−1and consequently provide a segmentation of the original time series {Xt,T }Tt=0. We place
the following technical assumption on the breakpoints present in the EWS.
Assumption 1. The set of those locations z where (possibly infinitely many) functions
Si (z) contain a jump, is finite. That is, let B = {z; ∃ i limu→z− Si (u) 6= limu→z+ Si (u)},then B = |B| < ∞.
We further define the wavelet periodogram of the LSW time series.
Definition 2. Let Xt,T be an LSW process as in (1). The triangular stochastic array
2(i)

It,T =

X

Xs,T ψi,s−ts

is called the wavelet periodogram of Xt,T at scale i.

(2)

Time series segmentation

5(i)Fryzlewicz and Nason (2006) show that the expectation EIt,T of the wavelet periodogramis “close” (in the sense of the integrated squared bias converging to zero) to the functionP−1
βi (z) = j=−∞ Sj (z)Ai,j , where A is an invertible matrix related to the wavelet vectors in(1) and (2). Further, each function βi (z) is piecewise constant with at most B jumps, all(i)of which occur in the set B. Besides, EIt,T themselves are piecewise constant by definition,
except intervals of length O(2−i ) around the discontinuities occurring in B.
The finiteness of B implies that there exists a fixed index I ∗ < ⌊log2 T ⌋ such that eachbreakpoint in B can be found in at least one of the functions Si (z) for i = −1, . . . , −I ∗ , i.e.,
no further new breakpoints appear at scales coarser than −I ∗ . Thus, from the invertibility(i)of A and the closeness of βi (z) and EIt,T , we conclude that every breakpoint is encoded in
the wavelet periodogram sequences at scales i = −1, . . . , −I ∗ . Therefore our procedure fordetecting breakpoints in the autocorrelation structure of Xt,T only considers the wavelet periodograms at these scales. Since I ∗ is fixed but unknown, in our theoretical considerations
we permit it to increase slowly to infinity with T .
A further reason for disregarding the coarse scales i < −I ∗ is that the autocorrelation
within each wavelet periodogram sequence becomes stronger at coarser scales; similarly, the(i)intervals on which EIt,T is not piecewise constant become longer. Thus, for coarse scales,
wavelet periodograms provide little useful information about breakpoints and can safely beomitted.
We end this section by briefly summarising our reasons behind the choice of the LSWmodel as a suitable framework for developing our methodology:(i) The entire second-order structure of the process, which varies over time in a piecewiseconstant manner, is encoded in the (asymptotically) piecewise constant expectationsof the wavelet periodogram sequences. Thus, any breakpoints in the second order
structure must by definition be detectable by analysing the wavelet periodograms,
which are relatively easy to handle as they follow a multiplicative model and are“localised” due to the compact support of the underlying wavelets.
(ii) Furthermore, due to the “whitening” properties of wavelets (see e.g. Vidakovic (1999),
Chapter 9), the wavelet periodogram sequences are often much less autocorrelated
than the original process.
(iii) The entire array of the wavelet periodograms at all scales is easily and rapidly computable via the non-decimated wavelet transform, which helps keep the computationalload of our procedure to the minimum.
(iv) Last but not least, the use of the “rescaled time” z = k/T in (1) and the associatedregularity assumptions on the transfer functions Wi (z) permit us to establish rigorous
asymptotic properties of our procedure.
3. Binary segmentation algorithm
In this section, noting that each wavelet periodogram sequence follows a multiplicativemodel as described in Section 3.1 below, we first introduce a binary segmentation algorithm
for such class of sequences. Binary segmentation is a computationally efficient tool which
searches for multiple breakpoints in a recursive manner (and can be classed as a “greedy”
and “multilevel” algorithm). Venkatraman (1993) applied the procedure to a sequence

6

Haeran Cho and Piotr Fryzlewicz

of independent normal variables with multiple breakpoints in its mean and showed thatthe detected breakpoints were consistent in terms of their number and locations. In thefollowing, we aim at extending these consistency results to the multiplicative model where,
in addition, dependence between observations is permitted.
Our binary segmentation procedure is then applied to each wavelet periodogram sequence separately, and we establish a procedure which efficiently combines the breakpoint
detection outcomes across the scales. We show that the consistency attained for each periodogram sequence carries over to the consistency in detecting breakpoints in the overall
second-order structure of the process.
3.1. Generic multiplicative model
Recall that each wavelet periodogram ordinate is simply a squared wavelet coefficient of(i)a zero-mean Gaussian time series, distributed as a scaled χ21 variable and satisfies It,T =
(i)

−1
2
, where {Zt,T }Tt=0
EIt,T · Zt,T
are autocorrelated standard normal variables. Hence we firstdevelop a generic breakpoint detection tool for multiplicative sequences defined by
2
2
2
Yt,T
= σt,T· Zt,T, t = 0, . . . , T − 1,(3)

where the following additional conditions hold.
2(i) σt,T
is deterministic and “close” to a piecewise constant function σ 2 (t/T ) in the sense∗2that σt,Tis piecewise constant apart from intervals of length at most O(2I ) around
P
−1
2
|σt,T
− σ 2 (t/T )|2 = o(log−1 T ), where thethe discontinuities in σ 2 (z), and T −1 Tt=0
latter rate comes from the rate of convergence of the integrated squared bias between(i)βi (t/T ) and EIt,T (see Fryzlewicz and Nason (2006) for details), and from the factthat our attention is limited to the I ∗ finest scales only. Further, σ 2 (z) is boundedfrom above and away from zero, with a finite but unknown number of jumps.
−1(ii) The vector {Zt,T }Tt=0
is multivariate normal with mean zero and variance one, and its
autocorrelation sequence is absolutely summable asymptotically;Pthat is the function
ρ(τ ) = supt,T |corr(Zt,T , Zt+τ,T )| satisfies ρ1∞ < ∞ where ρp∞ = τ ρp (τ ).

Once the breakpoint detection algorithm for the generic model (3) has been established, we2
apply it to the wavelet periodograms; thus, the reader is invited to relate Yt,T
, σ 2 (z) and
(i)

(i)

2
σt,T
to It,T , βi (z) and EIt,T , respectively.

3.2. Algorithm
The first step of the binary segmentation procedure is to find the likely location of a breakpoint. We locate such a point in the interval (0, T − 1) as the one which maximizes theabsolute value of
r
r
ν−1
T
−1
X
T −ν X 2νν
Y0,T −1 =Yt,T −Y2 .
(4)T · ν t=0
T · (T − ν) t=ν t,T
Yν0,T −1 can be interpreted as a scaled difference between the partial means of two segments
ν−1
−1
2
2
{Yt,T
}t=0
and {Yt,T}Tt=ν
, where the scaling is chosen so as to keep the variance Yν0,T −1

Time series segmentation

7

2
constant over ν in the idealised case of Yt,T
being i.i.d. Thus, it is a version of the wellknown CUSUM statistic. Once such a ν has been found, we use Yν0,T −1 (but not only thisquantity; see below for details) to test the null hypothesis of σ 2 (t/T ) being constant over
2 T −1{Yt,T}t=0 . The test statistic and its critical value, or threshold, are established in sucha way that when a breakpoint is present, the null hypothesis is rejected with probability
converging to 1. If the null hypothesis is rejected, we continue the simultaneous locating andtesting of breakpoints on the two segments to the left and to the right of ν in a recursive
manner until no further breakpoints are detected. The algorithm is summarised below,
where j is the level index and l is the location index of the node at each level.

Algorithm
Step 1 Begin with (j, l) = (1, 1). Let sj,l = 0 and ej,l = T − 1.
Step 2 Iteratively compute Ybsj,l ,ej,l as in (4) for b ∈ (sj,l , ej,l ), and find bj,l which maximizes its absolute value. Let nj,l = ej,l − sj,l + 1, dj,l = Ybsj,l ,ej,l , and mj,l =
Pej,l
2 √
t=sj,l Yt,T / nj,l .
p
Step 3 Perform hard thresholding on |dj,l |/mj,l with the threshold defined as tj,l = τ T θ log T /nj,lso that dˆj,l = dj,l if |dj,l | > mj,l · tj,l , and dˆj,l = 0 otherwise. The choice of θ and τ isdiscussed in Section 3.4.
Step 4 If either dˆj,l = 0 or max{bj,l −sj,l +1, ej,l −bj,l } < ∆T for l, stop the algorithm on theinterval [sj,l , ej,l ]. If not, let (sj+1,2l−1 , ej+1,2l−1 ) = (sj,l , bj,l ) and (sj+1,2l , ej+1,2l ) =(bj,l + 1, ej,l ) and update the level j as j → j + 1. Again, the choice of ∆T is discussed
in Section 3.4.
Step 5 Repeat Steps 2–4.
The set of detected breakpoints is {bj,l ; dˆj,l 6= 0}. The test statistic |dj,l |/mj,l is a scaled
version of the test statistics in the ICSS algorithm (Inclán and Tiao, 1994; Whitcher et al.,2000). However, the test criteria in those papers are derived empirically under the assumption of independent observations and there is no guarantee that either algorithm producesconsistent breakpoint estimates. On the other hand, our algorithm permits the sequences
to be autocorrelated as in (3), and its test criterion enables the consistent identification of
the total number and locations of breakpoints as shown in Section 3.3, provided that thetrue breakpoints are sufficiently scattered over time.
We also note the similarity between the statistic |dj,l |/mj,l (with the convention 0/0 = 0)
and the Fisz transform of Fryzlewicz and Nason (2006) and Fryzlewicz et al. (2006), thedifference being that the Fisz transform was only defined for the case b = 21 (ej,l + sj,l + 1)
(meaning the segments were split in half) and that it was not used for the purposes ofbreakpoint detection.

3.2.1. Post-processing within a sequence
We further equip the procedure with an extra step aimed at reducing the risk of overestimating the number of breakpoints. The ICSS algorithm in Inclán and Tiao (1994) has a“fine-tune” step whereby if more than one breakpoint is found, each breakpoint is checkedagainst the adjacent ones to reduce the risk of overestimation. We propose a post-processing

8

Haeran Cho and Piotr Fryzlewicz

procedure performing a similar task within the single-sequence multiplicative model (3). At
each breakpoint, the test statistic is re-calculated over the interval between two neighbouring breakpoints and compared with the threshold again. Denote the breakpoint estimatesas η̂p , p = 1, . . . , N̂ and η̂0 = 0, η̂N̂ +1 = T . For each η̂p , we examine whetherη̂Yη̂pp−1 +1,η̂p+1 > τ T θplog T ·

1(η̂p+1 − η̂p−1 )

η̂p+1

X

2
Yt,T
.

t=η̂p−1 +1

If the above inequality does not hold, η̂p is removed and the same procedure is repeated withthe reduced set of breakpoints until the set does not change. Note that the fine-tune step of
the ICSS algorithm re-calculates both the location and test statistic at each iteration andtherefore the locations of breakpoints are subject to change after tuning. However in our
post-processing procedure, only the test statistic is re-calculated at existing breakpoints,
and thus their locations are preserved. We note that the post-processing does not affect thetheoretical consistency of our procedure as (a) the extra checks above are of the same form
as those done in the original algorithm of Section 3.2, (b) the locations of the breakpoints
that survive the post-processing are unchanged. The next section provides details of ourconsistency result.
3.3. Consistency of detected breakpoints
In a breakpoint detection problem, it is desirable that the proposed procedure should (asymptotically) correctly identify the total number and locations of breakpoints. In thissection, we first show the consistency of our algorithm for a multiplicative sequence as in(3), which corresponds to the wavelet periodogram sequence at a single scale. Later, Theorem 2 will show how this consistency result carries over to the consistency of our procedure
in detecting breakpoints in the entire second-order structure of the input LSW process Xt,T .
Denote the number of breakpoints in σ 2 (t/T ) by N and the breakpoints themselves by
0 < η1 < . . . < ηN < T − 1, with η0 = 0, ηN +1 = T − 1. The following assumption statesthat the breakpoints ηp should be sufficiently scattered over time without being too close
to each other.
Assumption 2. For θ ∈ (1/4, 1/2) and Θ ∈ (θ + 1/2, 1), the length of each segment in
σ 2 (t/T ) is bounded from below by δT = O(T Θ ). Further, the breakpoints cannot be too close
to each other, i.e., there exists some constant c > 0 such that,)(ssηp − ηp−1
ηp+1 − ηp
≤ c.
,
max
1≤p≤N
ηp+1 − ηpηp − ηp−1−1
Theorem 1. Suppose that {Yt,T }Tt=0
follows model (3). Assume that• inf p σ 2 (ηp /T ) − σ 2 (ηp+1 /T ) ≥ ξ for some ξ > 0,
• supt |σ 2 (t/T )| ≤ M for some M > 0.

Under Assumption 2, the number and locationsof detected breakpoints are consistent.
Thatnois, Pr(A) → 1 as T → ∞ where A = N̂ = N ; |η̂p − ηp | ≤ O(ǫT ), 1 ≤ p ≤ N ; η̂p , p =1, . . . , N̂ are detected breakpoints and ǫT = T 1/2 log T . (Interpreting this in the rescaled
time interval [0, 1], ǫT /T = T −1/2 log T .)

Time series segmentation

9

3.3.1. Post-processing across the scales(i)
As mentioned in Section 2, we only consider wavelet periodograms It,T at scales i =∗2−1, . . . , −I ∗ , choosing I ∗ to satisfy 2I ≪ ǫT = T 1/2 log T so that the bias between σt,T2and σ (t/T ) does not cause us harm in deriving the results of Theorem 1. Recall from ourprevious discussion that any breakpoint in the second-order structure of the original process
Xt,T must be reflected in a breakpoint in at least one of the asymptotic wavelet periodogram
expectations βi (z) for i = −1, . . . , −I ∗ , and vice versa: a breakpoint in one of the βi (z)’s
implies a breakpoint in the second-order structure of Xt,T . Thus, it is sensible to combine
the estimated breakpoints across the periodogram scales by, roughly speaking, selecting a
breakpoint as significant if it appears in any of the wavelet periodogram sequences. Thissection describes a precise algorithm for doing this, and states a consistency result for the
final set of breakpoints arising from combining them across scales.
We first provide a general description of the algorithm, and then make it more precise
below. In the first stage, we produce a union of the sets of all breakpoints from scalesi = −1, . . . , −I ∗ . In the second stage, we cluster them into groups depending on theirdistance from one another. If there is a single scale containing representatives in eachgroup, we choose those as our final breakpoint estimates. If not, then the representative
breakpoint chosen from each group is the one that corresponds to the finest periodogram
scale at which it appears.
on(i)
The complete algorithm follows. Let B̂i = η̂p , p = 1, . . . , N̂i be the set of detected(i)∗breakpoints from the sequence It,T . Then the post-processing finds a subset of ∪−I
i=−1 B̂i ,say B̂, as formulated below;

Step 1 Arrange all breakpoints into groups so that those from different sequences andwithin the distance of ΛT from each other are classified as belonging to the samegroup; and denote the groups by G1 , . . . , GB̂ .
onStep 2 Find i0 = max arg max−I ∗ ≤k≤−1 N̂k , i.e., the finest scale with the most breakpoints.
(i )

(i)

Step 3 Check whether there exists η̂p00 for every η̂p , i 6= i0 , 1 ≤ p ≤ N̂i , which satisfies(i)(i )η̂p − η̂p00 < ΛT . If so, let B̂ = B̂i0 and quit the post-processing.
Step 4 Otherwise let B̂ =(finest scale).

noν̂p , p = 1, . . . , B̂ where each ν̂p ∈ Gp with the maximum i

We set ΛT = ⌊ǫT /2⌋ in order to take into account the bias arising in deriving the resultsof Theorem 1. As argued previously, breakpoints detected at coarser scales are likely tobe less accurate than those detected at finer scales; therefore, the above algorithm prefers
the latter. The across-scales post-processing procedure preserves the number of “distinct”
breakpoints and also their locations determined in the algorithm. Hence the breakpoints in
set B̂ are still consistent estimates of true breakpoints in the second-order structure of theoriginal nonstationary process Xt,T .
We acknowledge that the above is not the only possible way of combining the breakpoints
across scales which is still permitted by the theory and preserves consistency; however, we
advocate the above algorithm due to its good practical performance.

10

Haeran Cho and Piotr Fryzlewicz

Denote the set of the true breakpoints in the second-order
structure of oXt,T by B =n{νp , p = 1, . . . , B}, and the estimated breakpoints by B̂ = ν̂p , p = 1, . . . , B̂ .

Theorem 2. Suppose that Xt,T satisfies Assumption 1 and that νp , 1 ≤ p ≤ B, arescattered over time as in Assumption 2. Further assume that the conditions in Theorem 1
hold for each βi (z). Then our estimated number and locations of breakpoints are consistent,
that isno
Pr B̂ = B; |ν̂p − νp | ≤ O(ǫT ), 1 ≤ p ≤ B → 1

as T → ∞.

3.4. Choice of ∆T , θ, τ and I ∗To ensure that each estimated segment is of sufficiently large length so as not to distort ourtheoretical results, ∆T is chosen so that ∆T ≥ O(ǫT ). However, in practice our method
works well√for smaller values of ∆T too; e.g. in the forthcoming simulation experiments,∆T = O( T ) is used. As θ ∈ (1/4, 1/2), we use θ = 0.256 (as we have found that themethod works best when θ is close to the lower end of its permitted range) and elaborateon the choice of τ below (noting that our asymptotic theoretical results hold for any fixed
positive τ ).
The selection of τ is not a straightforward task and to get some insight into the issue, a
set of numerical experiments was conducted. A vector of random variables X ∼ NT (0, Σ)
was generated where X = (X1 , . . . , XT )T , then it was transformed into sequences of wavelet(i)periodograms It,T . The covariance matrix satisfied Σ = (σi,j )Ti,j=1 where σi,j = ρ|i−j| suchthat with varying ρ, the variables {Xt }Tt=1 were either independent or correlated. Then we
found b ∈ (1, T ) which maximised
s
r
b
T
X
T − b X (i)
b
(i)
b
Ii =
It,T −
It,T ,
T · b t=1
T (T − b)
t=b+1

P
√
(i)
and computed Ui,ρ,T = Ibi · {T −1 Tt=1 It,T · T θ log T }−1 . This was repeated with a varying
covariance matrix (ρ = 0, 0.3, 0.6, 0.9) and sample size (T = 512, 1024, 2048), 100 times foreach combination.
The quantity Ui,ρ,T√ is the ratio between our test statistic |d1,1 |/m1,1 and the timedependent factor T θ log T appearing in our threshold defined in the Algorithm of Section
3.2. Ui,ρ,T is computed under the “null hypothesis” of no breakpoints being present in thecovariance structure of Xt , and its magnitude serves as a guideline as to how to select thevalue of τ for each scale i to prevent spurious breakpoint detection in the null hypothesiscase. The results showed that the values of Ui,ρ,T and their range tended to increase forcoarser scales, due to the increasing dependence in the wavelet periodogram sequences. In
comparison to the scale factor i, the parameters ρ or T had relatively little impact on Ui,ρ,T .
Based on the above numerical experiments, we propose to use different values of τ in
Step 3 of the Algorithm of Section 3.2 and in the within-scale post-processing procedure ofSection 3.2.1. Denoting the former by τi,1 and the latter by τi,2 , we choose τi,1 differently
for each i as the 95% quantile of Ui,ρ,T , and τi,2 as its 97.5% quantile. The numerical values
are summarised in Table 3.4.

Time series segmentation

11

Table 1. Values of τ for each scale i = −1, . . . , −6.
scale i
−1
−2
−3
−4
−5
−6
τi,1
0.40 0.50 0.65 0.80 0.95 1.25
τi,20.45 0.60 0.75 0.90 1.10 1.35

Finally, we discuss the choice of I ∗ , i.e. the coarsest wavelet periodogram scale atwhich we still apply our breakpoint detection procedure. Firstly, we detect breakpoints inwavelet periodograms at scales i = −1, . . . , −⌊log2 T /3⌋ and perform the across-scalepostnoprocessing as described in Section 3.3.1, obtaining the set of breakpoints B̂ = ν̂p , p = 1, . . . , B̂ .
Subsequently, for the wavelet periodogram at the next finest scale, we compute the quantity
Vp , p = 1, . . . , B̂ + 1 as

Vp =

max

ν∈(ν̂p−1 ,ν̂p )

q

Pν
(i)
ν̂p −ν
t=ν̂p−1 +1 It,T −(ν̂p −ν̂p−1 )·(ν−ν̂p−1 )
Pν̂p

q

Pν̂p
(i)
ν−ν̂p−1
t=ν+1 It,T(ν̂p −ν̂p−1 )·(ν̂p −ν)(i)ν̂p−1 +1 It,T /(ν̂p − ν̂p−1 )where ν̂0 = −1 and ν̂B̂+1 = T − 1. Note that Vp is again of the same form as our basic√
test statistic. Then Vp is compared to τi,1 · T θ log T to see whether there are any further(i)breakpoints yet to be detected in It,T which have not been included in B̂. (This step issimilar to our within-scale post-processing.) If there is an interval [ν̂p−1 + 1, ν̂p ] where Vp
exceeds the threshold, I ∗ is updated as I ∗ := I ∗ + 1 and the above procedure is repeated
to update B̂ until either no further changes are made or I ∗ ≥ ⌊log2 T /2⌋. This approach istheoretically justified by Lemma 6 in the Appendix and shown to work well in our simulation
study.

4. Simulation study
In Davis et al. (2006), the performance of the Auto-PARM was assessed and comparedwith the Auto-SLEX (Ombao et al., 2001) through simulation in various settings. The
Auto-PARM was shown to be superior to Auto-SLEX in identifying both dyadic and nondyadic breakpoints in piecewise stationary time series. Some examples from the paper wereadopted for the comparative study between our method and the Auto-PARM, alongsidesome other new examples. In the simulations below, wavelet periodograms were computed
using Haar wavelets and both post-processing procedures (Section 3.2.1 and Section 3.3.1)
followed the application of the segmentation algorithm. In our examples, T = 1024 andtherefore I ∗ was set as 3 at the start of each application of the algorithm, then updated
automatically if necessary, as described in Section 3.4. The algorithm was coded in R andtook 0.32 seconds on average to run on a 2.40-GHz Intel Core 2 processor when analysinga realisation of length T = 1024.

4.1. Piecewise stationary AR process with clearly observable changes
This example is taken from Davis et al. (2006). The target nonstationary process was
generated from (5),12

Haeran Cho and Piotr Fryzlewicz

 0.9Xt−1 + ǫt
1.68Xt−1 − 0.81Xt−2 + ǫt
Xt =

1.32Xt−1 − 0.81Xt−2 + ǫt

for 1 ≤ t ≤ 512,
for 513 ≤ t ≤ 768,
for 769 ≤ t ≤ 1024

(5)

where ǫt ∼ i.i.d. N (0, 1) in all examples. As seen in Figure 1 (a), there is a clear visual difference between the three segments in the model. Figure 1 (b) shows the wavelet periodogramat scale −4 and the estimation results, where the dotted lines indicate the true breakpoints(η1 = 512, η2 = 768) while the dashed lines indicate the detected ones (η̂1 = 512, η̂2 = 758).
Note that although initially the procedure returned three breakpoints, the within-sequence
post-processing (Section 3.2.1) successfully removed the false one. The experiment was
repeated 100 times and the final counts of breakpoints are given in Table 4.4.
(b)

0

−10

−5

500

0

5

1000

10

1500

(a) X_t

0

200

400

600

800

1000

0

200

time

400

600

800

1000

time

Fig. 1. (a) A realisation of model (5) and true (dotted) and detected (dashed) breakpoints; (b) Wavelet
periodogram at scale −4 and the breakpoint detection outcome, showing also the spurious breakpoint which was later removed in post-processing.

4.2. Piecewise stationary AR process with less clearly observable changes
In this example, the piecewise stationary AR model is revisited, but its breakpoints are lessclear-cut, as seen in Figure 2 (a).

 0.4Xt−1 + ǫt
−0.6Xt−1 + ǫtXt =

0.5Xt−1 + ǫt

for 1 ≤ t ≤ 400,
for 401 ≤ t ≤ 612,
for 613 ≤ t ≤ 1024

(6)

Figure 2 (b) shows the wavelet periodogram at scale −1 for the realisation in the left
panel and also its breakpoint estimates (η̂1 = 372, η̂2 = 622). Both procedures achieved asimilarly good performance: our method accurately detected the two breakpoints in 97% of
the cases while the Auto-PARM performed well for all cases.

Time series segmentation

13

(b)

0

−3

−2

5

−1

0

10

1

15

2

3

20

(a) X_t

0

200

400

600

800

1000

0

200

400

time

600

800

1000

time

Fig. 2. (a) A realisation of model (6) and true (dotted) and detected (dashed) breakpoints; (b) Wavelet
periodogram at scale −1 and the breakpoint detection outcome.

4.3. Piecewise stationary AR process with a short segment
This example is again from Davis et al. (2006). A single breakpoint occurs and one segment
is much shorter than the other.

Xt =



0.75Xt−1 + ǫt
−0.5Xt−1 + ǫt

for 1 ≤ t ≤ 50,
for 51 ≤ t ≤ 1024.

(7)

A typical realisation of (7), its wavelet periodogram at scale −3, and the estimation
outcome are shown in Figure 3, where the jump at η1 = 50 was correctly identified asη̂1 = 49. Even though one segment is substantially shorter than the other, our procedure
was able to detect exactly one breakpoint in 94% of the cases and underestimation did notoccur even when it failed to detect exactly one.
4.4. Piecewise stationary unit-root-like process with changing variance
Financial time series, such as stock indices, individual share or commodity prices, or currency exchange rates, are for certain purposes (such as e.g. pricing of derivative instruments)
often modelled as random walk with a time-varying variance. Motivated by this, we generated a random-walk-like (but piecewise stationary) example from model (8) below, wherethe variance has two breakpoints over time and the AR parameter remains constant andvery close to 1.
 0.999Xt−1 + ǫt
0.999Xt−1 + 1.5ǫt
Xt =
0.999Xt−1 + ǫt

for 1 ≤ t ≤ 400,
for 401 ≤ t ≤ 750,
for 751 ≤ t ≤ 1024.

(8)

Recall that the Auto-PARM is designed to find the “best” combination of the total
number and locations of breakpoints and adopts a genetic algorithm to traverse the vast

14

Haeran Cho and Piotr Fryzlewicz
(b)

−4

0

−2

5

0

10

2

15

4

(a) X_t

0

200

400

600

time

800

1000

0

200

400

600

800

1000

time

Fig. 3. (a) A realisation of model (7) and true (dotted) and detected (dashed) breakpoints; (b) Wavelet
periodogram at scale −3 and the breakpoint detection outcome.

Table 2. Summary of breakpoint detection from simulations; Our
method (CF) and the Auto-PARM (AP).
number of breakpoints
model (5)model (6)model (7) model (8)CF AP CFAP
CF AP CF AP
0
0
0
1
0
0
1
1
41
1
0
94
99
1
29
1
0
0
2
90
99
97
100
6
0
94
18
3
10
1
1
0
0
0
4
10
4
0
0
0
0
0
0
0
0
0
0
0
0
0
2
5
0
0
total 100 100 100 100 100 100 100 100

parameter space. However, due to the stochastic nature of the algorithm, the Auto-PARM
occasionally fails to return consistent estimates, which is not a desirable feature. Thisinstability was emphasised in this example, each run often returning different breakpoints.
For one typical realisation, it detected t = 21, 797 as breakpoints, and then only t = 741
in the next run on the same sample path. Besides, from Figure 4, which summarises andcompares the empirical distribution of the detected breakpoints by both methods over 100
repetitions, it is clear that the performance of Auto-PARM leaves much to be desired forthis particular example, whereas our method performs very well. We emphasise that thisis not to be taken as a criticism of Auto-PARM in general, which performed very wellin our other examples. Note that it was at scale −1 of the wavelet periodogram thatboth breakpoints were consistently identified the most frequently. The computation of the
wavelet periodogram at scale −1 with Haar wavelets is a differencing operation and naturally“whitens” the almost-unit-root process (8), clearly revealing any changes of variance in thesequence. Finally, though inferior to Haar wavelets, performance achieved by our method
using other wavelets was also superior to Auto-PARM in this example.

Time series segmentation

20
15
0

5

10

frequency

15
10
0

5

frequency

15(b) AP

20(a) CF

0

200

400

600

800

time

1000

0

200

400

600

800

1000

time

Fig. 4. Selection frequency for each t; (a) Our method (CF); (b) Auto-PARM (AP).

5. U.S stock market data analysis
5.1. Dow Jones weekly closing values 1970–1975
The time series of weekly closing values of the Dow Jones Industrial Average index between
July 1971 and August 1974 was studied in Hsu (1979) and revisited in Chen and Gupta(1997). Historical data are available on www.google.com/finance/historical?q=INDEXDJX:.DJI, where daily and weekly prices can be extracted for any time period. Both papers
concluded that there was a change in the variance of the index around the third week ofMarch 1973. For the ease of computation of the wavelet periodogram, we chose the same
weekly index between 1 July 1970 and 19 May 1975 so that the data size was T = 256 andthe above-mentioned time period was contained in this interval. The third week of March
1973 corresponds to η = 141 and our procedure detected η̂ = 142 as a breakpoint, which isillustrated in Figure 5 (b).
5.2. Dow Jones daily closing values 2007–2009
We further investigated more recent daily data from the same source, between 8 January
2007 and 16 January 2009. Over this period, the global financial market experienced oneof the worst crises in history, covered extensively in the media. Our breakpoint detection
algorithm estimated two breakpoints (see the illustration in Figure 6), one occurring in
the last week of July 2007, and the other in mid-September 2008. The first breakpoint(η̂1 = 135) coincided with the outbreak of the worldwide “credit crunch” as subprime
mortgage backed securities were discovered in portfolios of banks and hedge funds aroundthe world. The second breakpoint (η̂2 = 424) coincided with the bankruptcy of Lehman
Brothers, a major financial services firm, (September 14, 2008), an event which broughteven more volatility to the market.
We note that Wikipedia (http://en.wikipedia.org/wiki/Financial_crisis_of_2007-2009)
also mentions these two dates as milestones of the crisis:“The financial crisis of 2007–2009 began in July 2007 when a loss of confidence
by investors in the value of securitized mortgages in the United States resulted

16

Haeran Cho and Piotr Fryzlewicz
(b)

0

600

500

700

1000

800

1500

900

2000

1000

2500

(a) X_t

1971

1973

1975

0

50

time

100

150

200

250

time

Fig. 5. (a) Weekly average values of the Dow Jones IA index (July 1970–May 1975); (b) Wavelet
periodogram at scale −1 and the breakpoint detection outcome.

in a liquidity crisis that prompted a substantial injection of capital into financial markets by the United States Federal Reserve, Bank of England and theEuropean Central Bank. The TED spread, an indicator of perceived credit riskin the general economy, spiked up in July 2007, remained volatile for a year,
then spiked even higher in September 2008, reaching a record 4.65% on October
10, 2008. In September 2008, the crisis deepened, as stock markets worldwide
crashed and entered a period of high volatility, and a considerable number ofbanks, mortgage lenders and insurance companies failed in the following weeks.”

(b)

0e+00

8000

9000

1e+05

2e+05

3e+05

4e+05

10000 11000 12000 13000 14000(a) X_t

2007

2008
time

2009

0

100

200

300

400

500

time

Fig. 6. (a) Daily average values of the Dow Jones IA index (Jan 2007–Jan 2009); (b) Wavelet
periodogram at scale −1 and the breakpoint detection outcome.

Time series segmentation

17

A. The proof of Theorem 1
The consistency of our algorithm is first proved for the sequence below,
2
2
Ỹt,T
= σ 2 (t/T ) · Zt,T
, t = 0, . . . , T − 1.

(9)Note that unlike in (3), the above model features the true piecewise constant σ 2 (t/T ).
Denote n = e − s + 1 and define
√
√
b
e
X
e−b
b−s+1 X 2
2
Ỹbs,e = √ √
Ỹt,T
− √ √
Ỹt,T ,n b − s + 1 t=sn e − b t=b+1
√
√
b
e
X
e−b
b−s+1 X 2
S̃bs,e = √ √σ 2 (t/T ) − √ √σ (t/T ),n b − s + 1 t=sn e − b t=b+1
√
√
b
e
X
e−b
b−s+1 X 2
2
σt,T
− √ √
σt,T .
Sbs,e = √ √n b − s + 1 t=sn e − b t=b+1
Note that the above are simply inner products of the respective sequences and a vector
whose support starts at s, is constant and positive until b, then constant negative until e,and normalised such that it sums to zero and sums to one when squared. Let s, e satisfy
ηp0 ≤ s < ηp0 +1 < . . . < ηp0 +q < e ≤ ηp0 +q+1 for 0 ≤ p0 ≤ B − q, which will always be
the case at all stages of the algorithm. In Lemmas 1–5 below, we impose at least one offollowing conditions:s < ηp0 +r − CδT < ηp0 +r + CδT < e for some 1 ≤ r ≤ q,{(ηp0 +1 − s) ∧ (s − ηp0 )} ∨ {(ηp0 +q+1 − e) ∧ (e − ηp0 +q )} ≤ ǫT ,

(10)
(11)

where (as elsewhere in the paper) C is an arbitrary positive constant and ∧ and ∨ are theminimum and maximum operators, respectively. We remark that both conditions (10) and(11) hold throughout the algorithm for all those segments starting at s and ending at ewhich contain previously undetected breakpoints. As Lemma 6 concerns the case when all
breakpoint have already been detected, it does not use either of these conditions.
The proof of the Theorem is constructed as follows. Lemma 1 is used in the proof of
Lemma 2, which in turn is used alongside Lemma 3 in the proof of Lemma 4. From theresult of Lemma 4, we derive Lemma 5 and finally, Lemmas 5 and 6 are used to prove
Theorem 1.
Lemma 1. Let s and e satisfy (10), then there exists 1 ≤ r ≤ q such that√ 
ηp0 +rS̃s,e= max |S̃ts,e | ≥ O δT / T .
s<t<e

(12)

Proof. The equality in (12) is proved by Lemmas 2.2 and 2.3 in Venkatraman (1993). For
the inequality, note first that in the case of a single breakpoint in σ 2 (z) we can use theconstancy of σ 2 (z) to the left and to the right of the breakpoint to show thatp√ηp0 +r − s + 1 e − ηp0 +r 2
ηp +r
√
=
S̃s,e0σ (ηp0 +r /T ) − σ 2 ((ηp0 +r + 1)/T ) ,n√ which is bounded from below by O δT / T . We remark that the order remains the samein the case of multiple breakpoints.
18

Haeran Cho and Piotr Fryzlewicz
ηp0 +r> 0 for some 1 ≤ r ≤ q.
Lemma 2. Suppose (10) holds, and further assume that S̃s,eηp0 +r≥ S̃bs,e + 2 log T.
Then for any b satisfying |ηp0 +r − b| = O(ǫT ) we have, for a large T , S̃s,e

Proof. Without loss of generality, assume ηp0 +r < b. As in Lemma 1, we first derive theresult in the case of a single breakpoint in σ 2 (z). The following holds;√pηp0 +r − s + 1 e − b ηp0 +rb√(13)S̃s,e , andS̃s,e = √e − ηp0 +r b − s + 1!
√
pηp0 +r − s + 1 e − b
ηp0 +r
ηp +r
b
√
1− √
S̃s,e − S̃s,e =S̃s,e0e − ηp0 +r b − s + 1qqb−η 0 +rb−η− 1 − e−ηpp0 +r
1 + ηp +rp−s+1
ηp +r
0
0 +r
q
=
S̃s,e0b−ηp0 +r
1 + ηp +r −s+10
b−ηp0 +rb−ηp0 +r
1
2 ηp0 +r −s+1 + e−ηp0 +r
ηp0 +r√≥ 2 log TS̃s,e≥2
for large T , applying the Taylor expansion in the last but one step, and Lemma 1 in thelast step. Similar arguments are applicable when b < ηp0 +r . Since the order of (13) remains
the same in the case of multiple breakpoints, the lemma is proved.
Lemma 3. Let n ≥ O(δT ). Assume ∃ b ∈ (s, e) satisfying
(r
)
r
b−s+1
e−b
max
≤ c.
,
e−b
s−b+1
Then Ỹbs,e − S̃bs,e ≤ log T with probability converging to 1 with T , uniformly over s, b, e.
Proof. We need to show that
e

Pr

1 X 2
2
√
σ (t/T )(Zt,T− 1) · ct > log T
n t=s

!

−→ 0,

(14)

√
√
√
√
where ct = e − b/ b − s + 1 for t ∈ [s, b] and ct = b − s + 1/ e − b otherwise. Let{Ut }et=s be i.i.d. standard normal variables, V = (vi,j )ni,j=1 with vi,j = corr (Zi,T , Zj,T ),and W = (wi,j )ni,j=1 be a diagonal matrix with wi,i = σ 2 (t/T ) · ct where i = t − s + 1. Bystandard results
e.g. Johnson and Kotz (1970), √page 151), showing (14) is equivalent to
P(see
e
2
showing that
λ
(U
−
1)
is
boundedbyn log T with probability converging to
t−s+1
t
t=s
1, where λi are eigenvalues of the matrix VW. Due to the Gaussianity of Ut , λt−s+1 (Ut2 −1)
satisfy the Cramér’s condition, i.e., there exists a constant C > 0 such thatE λt−s+1 (Ut2 − 1)p

2

≤ C p−2 p!E λt−s+1 (Ut2 − 1) , p = 3, 4, . . . .

Therefore we can apply Bernstein’s inequality (Bosq, 1998) and bound the probability in(14) by
n log2 T√
2 exp − Pn.
(15)
4 i=1 λ2i + 2 maxi |λi |C n log T

Time series segmentation

Pn

19

2

2
2
4
2
Note thati=1 λi = tr (VW) ≤ c maxz σ (z)nρ∞ . Also it follows that maxi |λi | ≤c maxz σ 2 (z)kVk where k · k denotes the spectral norm of a matrix, and kVk ≤ ρ1∞ since V
is non-negative definite. Then the lemma follows as (15) is bounded from above by
n log2 T√
exp − 2−→ 0,
4c maxz σ 4 (z)nρ2∞ + 2c maxz σ 2 (z) n log T ρ1∞∗

as ρp∞ ≤ O(2I ), which can be made to be of order O(log T ), since the only requirement on
I ∗ is that it converges to infinity but no particular speed is required.

Lemma 4. Assume (10) and (11). For b = arg maxs<t<e |Ỹts,e |, there exists 1 ≤ r ≤ qsuch that, for large T ,|b − ηp0 +r | ≤ ǫT .

(16)

Proof. Let S̃s,e = maxs<t<e |S̃ts,e |. From Lemma 3, Ỹbs,e ≥ S̃s,e − log T and S̃bs,e ≥ Ỹbs,e −log T , hence S̃bs,e ≥ S̃s,e − 2 log T . Assume (16) does not hold that b ∈ (ηp0 +r + ǫT , ηp0 +r+1 −ǫT ) for some r. From Lemma 2.2 in Venkatraman (1993), S̃ts,e is either monotonic orηp0 +r∨decreasing and then increasing between two breakpoints (on [ηp0 +r , ηp0 +r+1 ]), and S̃s,eηηp0 +r+1p+r> S̃bs,e . Suppose S̃s,e0 > S̃bs,e . Then there exists b′ ∈ (ηp0 +r , ηp0 +r + ǫT ] satisfying
S̃s,e′ηp +r+1
ηp +r> S̃bs,e (as S̃ts,e isS̃s,e0 − 2 log T ≥ S̃bs,e from Lemma 2. Since b > b′ , we also get S̃s,e0′′
locally increasing at t = b), and there will again be a b ∈ [ηp0 +r+1 − ǫT , ηp0 +r+1 ) satisfying′′
ηp +rS̃s,e0 − 2 log T ≥ S̃bs,e . Since b′′ > b, it contradicts that S̃bs,e ≥ S̃s,e − 2 log T . Similararguments are applicable when b < ηp0 +r and therefore the lemma follows.

Lemma 5. Under (10) and (11), for b = arg maxs<t<e |Ỹts,e |,
!
e
p
1X 2
b
θ
−→ 0.
Ỹ
Pr Ỹs,e < τ T
log T ·
n t=s t,T
Proof. From Lemma 4, there exists some r such that |b − ηp0 +r | < ǫT . Denote m̃ =
Pe
2
b
˜
˜
˜
t=s Ỹt,T /n and d = Ỹs,e = d1 − d2 where
√
√
b
e
X
e−b
b−s+1 X 22˜˜
d1 = √ √Ỹt,T and d2 = √ √Ỹt,T .
n b − s + 1 t=sn e − b t=b+1

˜
Further, let µi = Ed˜i and wi = var(pd˜i ) for i = 1, 2, and define µ = Ed˜ and w = var(d).
θ˜log T /n. We need to show Pr(|d| ≤ m̃ · tn ) → 0.
Finally, tn denotes the threshold τ T
We first note that wi ≤ c2 supz σ 4 (z)ρ2∞ . Using Markov’s and the Cauchy-Schwarz
inequalities, µ2 Pr(d˜ ≤ m̃ · tn ) is bounded bynoµ2 Pr (d˜1 − µ1 )(ctn − 1) + (d˜2 − µ2 )(ctn + 1) + 2ctn µ1 ≥ (1 + ctn )µ≤ 3(1 + ctn )−2 (ctn − 1)2 w1 + (ctn + 1)2 w2 + 4c2 tn2 µ21



4
2
2 2θ
,
≤ O sup σ (z) ρ∞ + τ T log T
z


√ 
√
and since µ = S̃bs,e = O δT / T > T θ log T , the conclusion follows.

20

Haeran Cho and Piotr Fryzlewicz

Lemma 6. For some positive constants C, C ′ , let s, e satisfy either(i) ∃ 1 ≤ p ≤ B such that s ≤ ηp ≤ e and [ηp − s + 1] ∧ [e − ηp ] ≤ CǫT or(ii) ∃ 1 ≤ p ≤ B such that s ≤ ηp < ηp+1 ≤ e and [ηp − s + 1] ∨ [e − ηp+1 ] ≤ C ′ ǫT .
Then for b = arg maxs<t<e |Ỹts,e | and large T ,
e
p
1X 2
Pr
log T ·> τTỸn t=s t,Tno√P2Proof. Let A = Ỹbs,e > τ T θ log T · et=s Ỹt,T/n andỸbs,e

B=

(

θ

!

−→ 0.

e


(ηp − s + 1)σ12 + (e − ηp )σ22
1 X 22<h=Ỹt,T − EỸt,T
n t=s
2n

)

,

where σ12 = σ 2 (ηp /T ) and σ22 = σ 2 ((ηp + 1)/T ). First we assume (i) holds. We have
Pr (A) = Pr (A |B ) Pr (B) + Pr (A |B c ) Pr (B c) ≤ Pr (A |B ) + Pr (B c ) .
The first part is bounded as
Pr
Pr (A |B ) ≤





P
√
2
−h
Ỹbs,e > τ T θ log T · n1 et=s EỸt,T
Pr(B)

.

(17)

We have2
2n · CǫT
max σ 4 (z)E Ỹbs,e ≤ n S̃bs,e ≤n − CǫT z
and n ≥ O(δT ). Therefore by applying Markov’s inequality, the numerator is bounded by2E Ỹbs,e /(τ h2 T 2θ log T ) ≤ O T 1/2−2θ log T −→ 0.
Turning our attention to the denominator of (17), we need to show that
!
e
X
1
2
σ 2 (t/T )(Zt,T− 1) > h −→ 0.
Pr(B c ) = Prn t=sThis can be shown by applying Bernstein’s inequality as in the proof of Lemma 2, and thelemma follows. Similar arguments are applied to prove the lemma when (ii) holds.

We now prove Theorem 1. At the start of the algorithm, as s = 0 and e = T − 1, allconditions for Lemma 5 are met and it finds a breakpoint within the distance of O(ǫT ) fromthe true breakpoint, by Lemma 4. Under Assumption 2, both (10) and (11) are satisfied
within each segment until every breakpoint in σ 2 (t/T ) is identified. Then, either of two
conditions (i) or (ii) in Lemma 6 is met and therefore no further breakpoint is detected withprobability converging to 1.
(i)2) affects the consistency. First we
Finally we study how the bias present in EIt,T (= σt,TP∞ψψ
define the autocorrelation wavelet
Ψ
(τ
)
=
i
k=−∞ i,k i,k+τ , the autocorrelation wavelet
P
Ψ
(τ
)Ψ
(τ
),
and
the across-scales autocorrelation wavelets
inner productmatrix
A
=
i
j
i,j
τ
P
Ψi,j (τ ) = k ψi,k ψj,k+τ . Then it is shown in Fryzlewicz and Nason (2006) that the inte(i)grated bias between EIt,T and βi (t/T ) converges to zero.

Time series segmentation

Proposition 1 (Propositions 2.1-2.2 (Fryzlewicz and Nason, 2006)).
be the wavelet periodogram at a fixed scale i. Under Assumption 1,
T −1

T
−1
X
t=0

(i)

EIt,T − βi (t/T )2= O(T −1 2−i ) + bi,T ,21(i)
Let It,T(18)

where bi,T depends on the sequence {Li }i . Further, each βi (z) is a piecewise constantfunction with at most B jumps, all of which occur in the set B.
Suppose the interval [s, e] includes a true breakpoint ηp as in (10), and denote b =(i)arg maxt∈(s,e) |S̃ts,e | and b̂ = arg maxt∈(s,e) Sts,e . Recall that EIt,T remains constant within−ieach stationary segment, apart from short (of length O(2 )) intervals around the discontinuities in βi (t/T ). Suppose a jump occurs at ηp in βi (t/T ) yet there is no change in(i)
EIt,T for t ∈ [ηp − O(2−i ), ηp + O(2−i )]. Then the integrated bias is bounded from below
by O(δT /T ) from Assumption 2, and Proposition 1 is violated. Therefore there will be a
(i)
(i)
(i)
change in EIt,T as well on such intervals around ηp and EIt1 ,T 6= EIt2 ,T for t1 ≤ ηp − O(2−i )(i)and t2 ≥ ηp + O(2−i ). Although the bias of EIt,T in relation to βi (t/T ) may cause some∗bias between b̂ and b, we have that |b̂ − b| ≤ O(2I ) < ǫT holds for I ∗ = O(log log T ), which
is an admissible rate for I ∗ . Besides, once one breakpoint is detected in such intervals, thealgorithm does not allow any more breakpoints to be detected within the distance of ∆T(i)from the detected breakpoint, by construction. Hence the bias in EIt,T does not affect theresults of Lemmas 1–6 for wavelet periodograms at finer scales and the consistency still
2
2
2
holds for Yt,T
= σt,T· Zt,T
as in (3).
B. The proof of Theorem 2

From Assumption 1 and the invertibility of the autocorrelation wavelet inner product matrix(i)A, there exists at least one sequence of wavelet periodograms among It,T , i = −1, . . . , −I ∗in which any breakpoint in B is detected. Suppose there is only one such scale, i0 , for(i )νq ∈ B and denote the detected breakpoint as η̂p00 . After the across-scales post-processing,(i)(i )η̂p00 is selected as ν̂q since no other η̂p , i 6= i0 , is within the distance of ΛT = O(ǫT ) from
(i )

(i )

either ν̂q or η̂p00 , and νq − η̂p00∗≤ ǫT with probability converging to 1 from Theorem 1.
(i )

(i )

If there are D(≤ I ) breakpoints detected for νq , denote them as η̂p11 , . . . , η̂pDD . Then for
(i )

(i )

(i )

(i )

any 1 ≤ a < b ≤ D, η̂paa − η̂pbb ≤ η̂paa − νq + η̂pbb − νq ≤ O(ǫT ), and only the one
from the finest scale is selected as ν̂q among them by the post-processing procedure. Hencethe across-scales post-processing preserves the consistency for the breakpoints selected as
its outcome.
References
Adak, S. (1998). Time-dependent spectral analysis of nonstationary time series. J. Am.
Stat. Assoc. 93, 1488–1501.
Andreou, A. and E. Ghysels (2002). Detecting multiple breaks in financial market volatility
dynamics. J. Appl. Econ. 17, 579–600.

22

Haeran Cho and Piotr Fryzlewicz

Bai, J. and P. Perron (1998). Estimating and testing linear models with multiple structuralchanges. Econometrica 66, 47–78.
Bhattacharya, P. (1994). Some aspects of change-point analysis. Change-point problems(eds E. Carlstein, H.-G. Müller and D. Siegmund) 23, 28–56.
Bosq, D. (1998). Nonparametric statistics for stochastic process: estimation and prediction.
Springer.
Chen, J. and A. K. Gupta (1997). Testing and locating variance change-points with application to stock prices. J. Am. Stat. Assoc. 92, 739–747.
Chernoff, H. and S. Zacks (1964). Estimating the current mean of a normal distribution
which is subject to changes in time. Ann. Math. Statist. 35, 999–1028.
Choi, H., H. Ombao, and B. Ray (2008). Sequential change-point detection method in timeseries. Technometrics 50, 40–52.
Davis, R. A., T. C. M. Lee, and G. A. Rodriguez-Yam (2006). Structural break estimation
for non-stationary time series. J. Am. Stat. Assoc. 101, 223–239.
Davis, R. A., T. C. M. Lee, and G. A. Rodriguez-Yam (2008). Break detection for a classof nonlinear time series models. J. Time Ser. Anal. 29, 834–867.
Fearnhead, P. (2005). Exact Bayesian curve fitting and signal segmentation. IEEE Trans.
Signal Process. 53, 2160–2166.
Fryzlewicz, P. and G. Nason (2006). Haar-Fisz estimation of evolutionary wavelet spectra.
J. Roy. Stat. Soc. B 68, 611–634.
Fryzlewicz, P., T. Sapatinas, and S. Subba Rao (2006). A Haar-Fisz technique for locallystationary volatility estimation. Biometrika 93, 687–704.
Gabbanini, F., M. Vannucci, G. Bartoli, and A. Moro (2004). Wavelet packet methods forthe analysis of variance of time series with application to crack widths on the Brunelleschi
Dome. J. Comput. Graph. Stat. 13, 639–658.
Hawkins, D. M. (1977). Testing a sequence of observations for a shift in location. J. Am.
Stat. Assoc. 72, 180–186.
Hsu, D. A. (1977). Tests for variance shifts at an unknown time point. Appl. Statist. 26,
179–184.
Hsu, D. A. (1979). Detecting shifts of parameters in gamma sequences with applications tostock price and air traffic flow analysis. J. Am. Stat. Assoc. 74, 31–40.
Inclán, C. and G. C. Tiao (1994). Use of cumulative sums of squares for retrospective
detection of changes of variance. J. Am. Stat. Assoc. 89, 913–923.
Janeway, W. (2009). Six impossible things before breakfast: Lessons from the crisis. Significance 6, 28–31.
Johnson, N. and S. Kotz (1970). Distributions in Statistics: Continuous Univariate Distributions, Vol. 1. Houghton Mifflin Company.

Time series segmentation

23

Kokoszka, P. and R. Leipus (2000). Change-point estimation in ARCH models. Bernoulli 6,
513–539.
Lavielle, M. and E. Moulines (2000). Least-squares estimation of an unknown number ofshifts in a time series. J. Time Ser. Anal. 21, 33–59.
McCulloch, R. E. and R. S. Tsay (1993). Bayesian inference and prediction for mean andvariance shifts in autoregressive time series. J. Am. Stat. Assoc. 88, 968–978.
Nason, G. P., R. von Sachs, and G. Kroisandt (2000). Wavelet processes and adaptive
estimation of the evolutionary wavelet spectrum. J. Roy. Stat. Soc. B 62, 271–292.
Ombao, H., J. Heo, and D. Stoffer (2004). Online analysis of seismic signals: An almostreal-time approach. New York: Springer.
Ombao, H. C., J. A. Raz, R. von Sachs, and B. A. Malow (2001). Automatic statisticalanalysis of bivariate nonstationary time series. J. Am. Stat. Assoc. 96, 543–560.
Punskaya, E., C. Andrieu, A. Doucet, and W. J. Fitzgerald (2002). Bayesian curve fitting
using MCMC with applications to signal segmentation. IEEE Transactions on Signal
Processing 50, 747–758.
Sen, A. and M. S. Srivastava (1975). On tests for detecting change in mean. Ann. Statist. 3,
98–108.
Van Bellegem, S. and R. von Sachs (2008). Locally adaptive estimation of evolutionary
wavelet spectra. Ann. Statist. 36, 1879–1924.
Venkatraman, E. S. (1993). Consistency results in multiple change-point problems. PhD
thesis, Stanford University.
Vidakovic, B. (1999). Statistical modeling by wavelets. Wiley.
Vostrikova, L. J. (1981). Detecting ‘disorder’ in multidimensional random processes. Sov.
Math. Doklady. 24, 55–59.
Whitcher, B., S. D. Byers, P. Guttorp, and D. B. Percival (2002). Testing for homogeneityof variance in time series: Long memory, wavelets, and the Nile River. Water Resour.
Res. 38, 10–1029.
Whitcher, B., P. Guttorp, and D. B. Percival (2000). Multiscale detection and location ofmultiple variance changes in the presence of long memory. J. Statist. Comput. Simul. 68,
65–87.
Worsley, K. J. (1986). Confidence regions and tests for a change-point in a sequence of
exponential family random variables. Biometrika 73, 91–104.