1

Detecting Changes in Time Series Data using
Volatility Filters

arXiv:1709.03105v2 [cs.SY] 28 Dec 2017

Alireza Ahrabian, Nazli Farajidavar, Clive Cheong-Took and Payam Barnaghi

Abstract
This work develops techniques for the sequential detection and location estimation of transient changes in the volatility
(standard deviation) of time series data. In particular, we introduce a class of change detection algorithms based on the windowed
volatility filter. The first method detects changes by employing a convex combination of two such filters with differing window
sizes, such that the adaptively updated convex weight parameter is then used as an indicator for the detection of instantaneous
power changes. Moreover, the proposed adaptive filtering based method is readily extended to the multivariate case by using
recent advances in distributed adaptive filters, thereby using cooperation between the data channels for more effective detection
of change points. Furthermore, this work also develops a novel change point location estimator based on the differenced output
of the volatility filter. Finally, the performance of the proposed methods were evaluated on both synthetic and real world data.

Index Terms— Volatility Change Detection, Transient Detection, Multi-Sensory/Multivariate Data.

I. I NTRODUCTION
The piecewise segmentation of time series data into intervals with similar stochastic characteristics is often important in the
exploratory analysis of data and is generally referred to as change detection. With rapid growth in sensing technologies, there is
a need for effective and adaptive methods to identify and analyse activities in time series data [1]. In particular the partitioning
of data into segments of constant volatility (variance/standard deviation) has interesting applications for fields ranging from
the segmentation of biomedical data [2] to audio segmentation [3] [4].
Algorithms proposed for estimating the transition times for data points with transient changes in the variance are generally
divided into two categories: i) Offline algorithms that use the entire data set so as to estimate the location of the change points
and ii) online (sequential) methods, where data is only available up to the current time sample being analysed. Early research
primarily focused on developing statistical significance based change detection algorithms, where in particular the work in [5]
proposed a method based on the cumulative sum of squares in combination with a statistical significance test to provide an
estimate of change points. More recently, algorithms based on Markov Chain Monte Carlo (MCMC) techniques for obtaining
estimates of the change point locations have been developed [6] [7] [8]. However, such methods are computationally expensive
and process the data offline. More computationally efficient block based dynamic programming [9] and Bayesian methods [10]
The authors Alireza Ahrabian and Payam Barnaghi are with the Department of Electrical and Electronic Engineering, Institute for Communication Systems
at the University of Surrey (email: {a.ahrabian, p.barnaghi}@surrey.ac.uk).
Clive Cheong-Took is with the Department of Computer Science at the University of Surrey (email:c.cheongtook@surrey.ac.uk).
Nazli Farajidavar is with the Computational Health Informatics laboratory at the University of Oxford (email:nazli.farajidavar@eng.ox.ac.uk).

2

Fig. 1: Time series data with instantaneous power changes, where the vertical dashed lines show the location of the transition
points.

have been proposed, however, they have limited use for online detection of change points. It should be noted, that block based
multivariate change point detection algorithms have also been developed. In particular the work in [11] developed an efficient
Bayesian approach, while the work in [12] used the group fused LASSO.
Many real world practical change detection problems require sequential algorithms that are not only capable of detecting
with the smallest possible delay the change volatility, but also providing a good estimate of the change point location. The
work in [13] introduced a method based on the likelihood ratio test that includes a simple threshold and test statistic for
both detecting and estimating the location of the change in the volatility of the time series data. Furthermore, the work in [7]
proposed a threshold independent sequential change detection method using the minimum description length (MDL), while
an accurate online kernel-Support Vector Machine (k-SVM) based method was proposed in [3]. However, such methods often
have drawbacks, including: high computational complexity, low accuracy of change point location estimate and high detection
latency.
We introduce a class of change detection and location estimation algorithms based on the windowed volatility filter [14].
The first proposed algorithm seeks to detect such transient power/volatility variations as rapidly as possible. This is achieved
by employing a convex combination of two such filters, where a state change is detected by selecting the appropriate filter
given the statistics of the data points. The proposed method is designed based on the trade off that exists between the steady
state accuracy and convergence speed of volatility filters. By using an adaptive filter [15] [16], the proposed method seeks
to overcome this trade off by updating the convex weight parameter that combines the respective windowed volatility filters,
thereby providing an indicator (in the form of the convex weight parameter) for detecting changes in the volatility of the
data. An advantage of the use of adaptive filters for the change point detection problem is a multivariate extension of the
proposed algorithm that is readily developed by using recent advances in distributed signal processing [17]. Furthermore, this
work analyses the behaviour of the volatility filter weights on the performance of the change detection algorithm, thereby
providing a means of selecting such weights for improved performance. The second proposed method provides an estimate
for the change point location. This is achieved by proposing a new estimator based on the differenced output of the volatility
filter. The performance of both the proposed change detection and location estimation techniques are verified on both synthetic
and real world data sets.

3

II. P ROBLEM O UTLINE
Given a sequence of data points {x1 , x2 , x3 , . . . , xt } drawn from a zero mean Gaussian distribution with piecewise constant
variance; there exists an unknown sequence of transition time instances {τ1 , τ2 , . . . , τt } at which a change in the variance
occurs (an example is shown in Fig. 1). Accordingly, this work seeks to address the following problems: 1) detecting a change
in the volatility of time series data with the shortest possible latency; 2) providing an accurate online estimate of the location
of the transition times {τ̂1 , τ̂2 , . . . τ̂t }.
III. P ROPOSED A LGORITHMS
This work proposes volatility filtering based change detection and location estimation algorithms given a zero mean set of
data points, x = [x1 , x2 , x3 , . . . , xN ]; where an estimate of volatility is given as the output of the volatility filter, σ̂p (t). In
particular the volatility filter is defined as the finite impulse response (FIR) filtered estimate of the past Tp squared data points,
that is
σ̂p2 (t) = wp,1 x2t + ... + wp,Tp x2t−Tp +1

(1)

= wTp xp
where wp = [wp,1 , ...wp,Tp ] correspond to the set of volatility filter weight coefficients with

PTp

i=1

wp,i = 1 (for the unbiased

estimate of variance this equality is approximately satisfied for large samples) and xp = [x2t , ..., x2t−Tp +1 ]. By using the volatility
filter as defined in (1), we first develop an adaptive filtering based change detection algorithm that utilises a convex combination
of two such distinct volatility filters. We also provide analysis for selecting appropriate filter weights for the respective volatility
filters. The second proposed method is a change point location estimator that is based on the output of the differenced volatility
filter.

A. Adaptive Filtering Based Change Detection Method
This section presents an adaptive filtering and volatility filtering based change detection algorithm. This is achieved by using
an adaptive filter [18] [15] [16] to track changes in the instantaneous standard deviation using a convex combination of two
volatility filters with different window sizes. Where a volatility filter with a small window size is able to effectively track
transient changes in the time series volatility; however, in the steady state the estimate in the volatility of the signal is not
accurate (please see Fig. 2). While a volatility filter with a large window size has stable steady state estimates of volatility,
however, it is not able to track changes effectively (please see Fig. 2). By adaptively estimating the convex weight parameter,
such that if a significant transient change has occurred then volatility filter with a small window size provides a more optimal
(least squared sense) solution than the volatility filter with the large window size. Implying that by tracking the estimate of
the convex weight parameter, it can be determined whether or not a change has occurred. Furthermore, the proposed change
detection method can be extended to data sets drawn from multiple sensors, by using an adaptive filtering architecture that

4

Algorithm 1: Adaptive Filtering Based Change Detection
Require:
- Input signal xt
- Select: µ, Ts , Tf , Td , γ.
- Set: λ(1) ← 1, t ← 1, ρ ← 0.001, Tr ← 1.2Ts .
while t < N do
- σo (t) ← λ(t)σ̂f (t) + (1 − λ(t))σ̂s (t).
- e(t) ← σ̂d (t) − σ̂o (t).
- λ(t + 1) ← λ(t) + µ(|λ(t)| + ρu)e(t)(σ̂f (t) − σ̂s (t))
if λ(t + 1) ≥ γ and Flag(t) = 0 then
- s(t) = 1 in the interval [t, t + Tr ].
µ
in the interval [t, t + Tr ].
- µ ← Var(x
t)
- Flag(t) = 1 in the interval [t, t + Tr ].
end if
- t←t+1
end while

utilises cooperative strategies [17] to enhance both convergence and steady state accuracy of the convex weight parameter.
Given a set of data points x, consider the following volatility filters
σ̂f2 (t) = wf,1 x2t + ... + wf,Tf x2t−Tf +1
= wTf xf
(2)
σ̂s2 (t) = ws,1 x2t + ... + ws,Ts x2t−Ts +1
= wTs xs
where σ̂f (t) and σ̂s (t) correspond1 to the respective outputs of the ‘fast’ and ‘slow’ volatility filters, such that, Tf < Ts . The
proposed estimate of the instantaneous volatility σ̂o (t) is given by a convex combination of the ‘fast’ and ‘slow’ filters, that is
σ̂o (t) = λ(t)σ̂f (t) + (1 − λ(t))σ̂s (t)

(3)

where λ(t) is updated adaptively using the adaptive filter architecture outlined in [15] [19]. The squared error, E{e2 (t)},
required to be minimised is generally given by the following
2
E{e2 (t)} = E{ σ̂d (t) − σ̂o (t) },

(4)

where σ̂d (t) is defined as
σ̂d2 (t) = w1 x2t+1 + ... + wTd x2t−Td

(5)

= wTd xd
where Td is the window length of the desired volatility filter. The desired signal σ̂d (t) provides a coarse estimate of the
instantaneous volatility so as to enable the adaptive filter to track changes in volatility; to this end, a typical window length of
Td ≈ 10 is recommended. It should be noted that λ(t) → 1 during a change in the volatility, that is, a higher weight is placed
1 We now refer to the volatility filters with relatively small and large window sizes as ‘fast’ and ‘slow’ filters respectively, this is to reflect the convergence
speed of the filters when transient changes occur in the signal volatility.

5

Volatility (σ)

Large Window
True Value

Time Index

Volatility (σ)

Small Window
True Value

Time Index

Fig. 2: Figure showing the estimated volatility for time series with instantaneous power changes. (Upper panel) Shows the
output of the ‘slow’ and (lower panel) the output of the ‘fast’ volatility filter.

on the ’fast’ filter. Moreover, in the steady state (variance does not change) a higher weight will be assigned to on the ‘slow’
volatility filter; however, during simulations we found large variations in the weight estimate λ(t) when using the conventional
least mean square (LMS) update. To this end, we propose to use signed sparse-least mean squares (SS-LMS) algorithm [20],
that is the minimised cost function (4) along with a prior distribution on the values taken by λ(t). By using a variation of
the SS-LMS algorithm we ensure that the estimate of λ(t) has lower variability as it converges to the ‘slow’ volatility filter.
Accordingly, the update of, λ(t), is given by2
λ(t + 1) = λ(t) + µ(|λ(t)| + ρu)e(t)(σ̂f (t) − σ̂s (t))

(6)

where µ is the learning rate of the adaptive filter, ǫ is a positive constant for regularization, ρ is a parameter selected so as
to reduce the variability of λ(t) and u corresponds to a Gaussian distributed random variable. The
 updated convex weight

 1, λ(t) ≥ γ
, where
parameter is then mapped to the output of the binary state change sequence s(t), that is, s(t) =

 0, λ(t) < γ
0 ≤ γ ≤ 1, is a threshold and s(t) = 1 corresponds to a change in the volatility of the time series, while s(t) = 0 corresponds to
no change. Algorithm 1 summarizes the proposed method referred to as the adaptive filtering based change detection (AFCD).
It should be noted that after a change has been detected, normalisation of the learning rate is carried out so as to ensure that
the proposed method is independent of the signal scale. Finally, it should be noted that the proposed method assumes that the
time index differences between change point locations are greater than Tr samples, owing to the time period required for the
estimation of the variance as well as the convergence of the filter to the ‘slow’ volatility filter after a change has been detected.

2A

hard threshold is applied so that 0 ≤ λ(t) ≤ 1, that is,
λ(t) =



1,
0,

λ(t) ≥ 1
λ(t) ≤ 0.

6

Algorithm 2: Cooperative Adaptive Filtering Based Change Detection
Require:
- Input signal xc (t)
- Select: µc , Ts,c , Tf,c , Td,c , γc .
- Set: λc (1) ← 1, t ← 1, ρc ← 0.001, Tr,c ← 1.2Ts,c .
PNc
λc (1)
ψ(1) ← c=1
while t < N do
for c = {1, . . . , Nsen } do
- σ̂o,c (t) ← ψ(t)σ̂f,c (t) + (1 − ψs (t))σ̂s,c (t).
- ec (t) ← σ̂d,c (t) − σ̂o,c (t).
- λc (t + 1) ← ψc (t) + µc (|ψ(t)| + ρc uc )ec (t)(σ̂f,c (t) − σ̂s,c (t))
end for
PNc
λc (t + 1)
ψ(t + 1) ← c=1
if ψ(t + 1) ≥ γc and Flag(t) = 0 then
- s(t) ← 1 in the interval [t, t + Tr ].
- µc ← Var(xµcc (t)) for [t, t + Tr ].
- Flag(t) = 1 in the interval [t, t + Tr ].
end if
- t←t+1
end while

Algorithm 2 proposes the multivariate change detection algorithm via an extension of Algorithm 1. That is given a
multichannel signal we seek to propose a cooperative strategy between the data channels for improved detection of transition
points. This is achieved by utilising the combine then adapt (CTA) strategy of the diffusion LMS which was proposed in
[17]; where given Nsen sensors such that each sensor c has an input signal xc (t), desired signal dc (t) and parameters vc (t),
where each sensor also has the ability to send and receive state estimates to one another3, accordingly the CTA strategy for
cooperative adaptive filtering is given by [17]
ψ c (t) =

X

gc,l vl (t − 1)

l∈Nc

vc (t) = ψ c (t) + µc dc −

(7)
ψ Tc (t)xc (t)


xc (t)

where µc is the learning rate of each sensor. The update of the LMS is carried out with respect to the averaged estimates of
the LMS output from the previous time instants, where Nc corresponds to the set of sensors connected to the cth sensor and
gc,l corresponds to the combiner of the sensor estimates. Algorithm 2 provides an overview of the cooperative adaptive filter
based change detection (CAFCD) algorithm for multivariate data.

B. Volatility Filter Weight Selection
In this section we analyse the most appropriate choice of weights for both the ‘fast’ wf and ‘slow’ ws volatility filters, for
detecting changes in time series data. Namely, given a sequence of zero mean data points, x1 , ...xN with a transition time at
t1 , where the variance of the first segment is given by σ12 and the variance of the second segment is given by σ22 . The objective
is to determine, how the choice of volatility filter weights affects the detection performance of the adaptive filtering based
change detection algorithm.
3 The diffusion LMS algorithm considers sensor networks, that is, all the sensors are not fully connected. However, we consider devices where the same
information is being recorded and thus we use full cooperation between all the sensors, in this work.

7

0.07

Filter Weight Value

0.06
0.05
0.04
0.03
0.02
0.01
0
1

5

10

15

20

25

Index
0.35

Filter Weight Value

0.3
0.25
0.2
0.15
0.1
0.05
0
1

2

3

4

5

Index

Fig. 3: Figure showing the triangular impulse response function for ‘slow’ filter (upper panel) and ‘fast’ filter (lower panel).

Note that the detection performance is dependent on the rate at which λ(t) → 1 after the transition time t1 . For rigour, we
calculate the minimum of the expectation of the cost function in (4) to analyse the evolution of λ(t). That is, we determine
the expected value of λ(t) that minimises the cost function at each time index, shown by the following
E{λ(t)} =

E {(σ̂d (t) − σ̂s (t))(σ̂f (t) − σ̂s (t))}
E {(σ̂f (t) − σ̂s (t))2 }

(8)

Given a change point at time t1 , the expression in (8) can be evaluated at three different time intervals: 1) t < t1 and
t > t1 + Ts , 2) t1 ≤ t ≤ t1 + Tf , and 3) t1 + Tf < t ≤ t1 + Ts .
Now the set of time intervals in 1) correspond to the case where both the ‘fast’ and ‘slow’ volatility filters are measuring data
points such that statistics of the data points are identical (strictly stationary). It is straightforward to show that E{λ(t)} ≈

1
2

< 1.

The time intervals in 2) and 3) correspond to the behaviour of λ(t) after a change has occurred. In order to evaluate the
performance of different window parameters on the behaviour of λ(t), we compute (8) over the interval4 2). Moreover, the
interval 2) corresponds to the case where the input data points to both the ‘fast’ and ‘slow’ volatility filters are equivalent
to a mixture of both the first segments statistics with expected variance σ12 and the second segments statistics with expected
variance σ22 . That is, the expected value of the output for the squared volatility filter with weights wf , is given by
E{σ̂f2 (t)} = wf1 (t)σ12 + wf2 (t)σ22
with E{σ̂f (t)} ≈

(9)

q
E{σ̂f2 (t)}, for t1 ≤ t ≤ t1 + Tf , where wf1 (t) is a monotonically increasing function and wf2 (t)

monotonically decreasing function, such that both functions are dependent on the volatility filter weights wf , with sum of the
4A

similar argument can be made for interval 3).

8

two functions equal to one. The corresponding output of the squared volatility filter with weights ws , is given by
E{σ̂s2 (t)} = ws1 (t)σ12 + ws2 (t)σ22
with E{σ̂s (t)} ≈

(10)

p
E{σ̂s2 (t)}, for t1 ≤ t ≤ t1 + Ts , where ws1 (t) and ws2 (t) are monotonically increasing and decreasing

functions respectively, with the sum of both functions equal to one. Accordingly, (8) evaluated over t1 ≤ t ≤ t1 + Tf , is given
by
E{λ(t)} ≈

A(t)
B(t)

(11)

where


A(t) ≈ σ2 E{σ̂f (t)} − E{σ̂s (t)} − Cov σ̂s2 (t), σ̂f2 (t)

(12)


B(t) ≈ E{σ̂f2 (t)} + E{σ̂s2 (t)} − 2Cov σ̂s2 (t), σ̂f2 (t)

(13)

+ E{σ̂s2 (t)} − E{σ̂f (t)}E{σ̂s (t)}

and

− 2E{σ̂f (t)}E{σ̂s (t)}


where E{σ̂d (t)} ≈ σ2 and Cov σ̂s2 (t), σ̂f2 (t) = ǫ, where ǫ → 0 for Ts >> Tf . From (9), (10) and (11) for λ(t) → 1, we

require E{σ̂f2 (t)} ≈ σ1 and E{σ̂s2 (t)} ≈ σ2 over the interval t1 ≤ t ≤ t1 + Tf . The analysis implies that for the ‘fast’ volatility
filter a higher weighting needs to be applied to the most recent samples in wf . While for the ‘slow’ volatility filter the converse
is true, lower weighting needs to be applied to the most recent samples in ws . An example of such a pair of filters with impulse
responses that satisfy the condition above are shown in Fig. 3, namely a pair of triangular weighted filters, that is
x2t + 2x2t−1 + ... + (Ts − 1)x2t−Ts +2 + Ts x2t−Ts +1
1 + 2 + ... + Ts

σ̂s2 (t) =

(14)

for the ’slow’ volatility filter weights, and
σ̂f2 (t) =

Tf x2t + (Tf − 1)x2t−1 + ... + 2x2t−Tf +2 + x2t−Tf +1
1 + 2 + ... + Tf

(15)

for the ’fast’ volatility filter weights.
C. Volatility Change Point Estimator
The second method referred to as the volatility change estimator (VCE), seeks to estimate the location of the change point
transition times. This is achieved by first obtaining the differenced5 estimate σ̂D (t) from the output of the volatility filter σ̂l (t),
as follows
σ̂D (t) = σ̂l (t) − σ̂l (t − Tl ),
where Tl corresponds to the window size and σ̂l (t) is estimated using the square window volatility filter weights wl =

(16)
1
Tl −1 1Tl

(where 1Tl is a vector of length Tl with each element equal assigned the value 1). In order to develop a change point location
estimator, we first observe the following property in the plot of σ̂D (t) with respect to sample indices. Consider a zero mean
5 That

is the difference between the output of the volatility filter at the current time index t and the output of the volatility filter at time index t − Tl .

9

0.6
0.4
0.2
0
-0.2
500

1000 1500 2000 2500 3000 3500 4000 4500 5000 5500 6000
Sample Index

Fig. 4: The output of the volatility filter (of window size T ) (dashed line).

piecewise stationary (with respect to variance) Gaussian distributed sequence of data points with a change in the variance at
the time index 3001. The corresponding plot of σ̂D (t) is shown in Fig. 4. Observe that σ̂D (t) oscillates (with an approximately
constant variation) around the zero until the change point at sample index 3001, where σ̂D (t) increases significantly (with
respect to the prior variation in σ̂D (t)) and at approximately 3000 + Tl samples the time series hits a maximum (shown by
the dashed line in Fig. 4), as a result the following estimator is proposed as a change point location estimate
τ̂ = arg max |σ̂D (t)| − Tl + 1.

(17)

t

While the proposed change point location estimate can be shown to be a reasonable unbiased estimator for a single volatility
change in the time series data (see the following section), a set of rules are required to deal with multiple change point locations.
That is, if a change has been detected by the AFCD algorithm, the estimator in (17) is applied in the interval between the
change point detection time to 2Tl samples after the detection. For multivariate data, either an estimate of the change point
location for each channel is obtained, or a single estimate of the change point location is obtained by averaging across the
individual channel estimates.
D. Change Point Location Estimator Analysis
In this section we verify that the proposed change point estimator (17) is a reasonable unbiased estimator for the transition
point location. Given a sequence of zero mean Gaussian random variables, x1 , . . . , xτ −1 ∈ N (0, σ12 ) and xτ , . . . , xN ∈
N (0, σ22 ), where the change point occurs at the time index τ , the following equality for (17) is satisfied




τ̂ = E arg max |σ̂D (t)| − Tl + 1
t

(18)

in which Tl is the volatility filter window size and arg maxt |σ̂D (t)| is equivalent to showing the following statement
E{|σ̂D (τ + Tl − 1)|} > E{|σ̂D (τ + i)|}

i 6= Tl − 1

(19)

where i is the set of time shifts such that there is no further transitions in volatility within i << N samples of τ . This condition
ensures that the maximum evaluated for arg maxt |σ̂D (t)|, is related to a single instantaneous volatility change. For multiple
instantaneous transitions in volatility, the spacing between any two such transition time points are assumed to be greater than

10

i > 2Tl . If the following volatility change conditions σ2 ≷ σ1 hold, then (19) is equivalent to
E{σ̂D (τ + Tl − 1)} ≷ E{σ̂D (τ + i)}

i 6= T − 1

(20)

For σ2 < σ1 (where an equivalent argument can be made for σ2 > σ1 ),
E{σ̂D (τ + Tl − 1)} = E{σ̂l (τ + Tl − 1)} − E{σ̂l (τ − 1)}
 v

v
u l −1 2
2

u
l −1
uTX
xτ +Tl −1−j  uTX
x
τ
−1−j
−E t
=E t

Tl − 1  
Tl − 1 
j=0

j=0

we next take the first order Taylor expansion

E{σ̂D (τ + Tl − 1)} ≈ E

 12 
 21
l −1
2
 TX

x
τ +T −1−j
τ −1−j
−E
Tl − 1   j=0 Tl − 1 


l −1
TX
x2


j=0

(21)

where it should be noted that the sum of squared standard normal distributed variables are from a chi-squared distribution,
that is,

TX
l −1
j=0

TX
l −1
j=0

x2τ +Tl −j
∽ χ2 (Tl − 1)
σ22
x2τ −j
∽ χ2 (Tl − 1)
σ12

As a result, E{σ̂D (τ + Tl − 1)} ≈ σ2 − σ1 . The corresponding estimates for E{σ̂D (t + Tl + j)} and E{σ̂D (t + Tl − j)} for
k = 2, ..., Tl are given by
E{σ̂D (τ + Tl + k)} ≈ σ2 −

E{σ̂D (τ + Tl − k)} ≈

s

s

(Tl − k)σ12 + (k − 1)σ22
Tl − 1

(Tl − k)σ22 + (k − 1)σ12
− σ1
Tl − 1

where it can be shown that E{σ̂D (τ + Tl − 1)} is less than both E{σ̂D (τ + Tl + k)} and E{σ̂D (τ + Tl − k)}, thus illustrating
that, arg maxt |σ̂D (t)| provides a reasonable unbiased estimate of the change point location.
IV. S IMULATIONS
The proposed change point detection and location estimation methods were analysed on both synthetic and real world data.
In particular, we assessed the performance of the proposed method with the generalized likelihood ratio (GLR) method [13],
owing to the similar computational costs of both algorithms.

A. Synthetic Signal
The first synthetic simulation consists of both univariate/multivariate time series data points x(n) drawn from a univariate/multivariate Gaussian distribution with instantaneous changes in the variance at random points in time; furthermore, each
realization of the multivariate data points had a randomly generated correlation matrix [21]. The instantaneous power changes
for each segment were selected randomly along with the number of samples in the segment as well as the total signal length.

11

Proportion of True Positives

1

0.95

0.9

0.85
Proposed Method - Triangular Weights
Proposed Method - Square Weights
GLR

0.8
1

2

3

4

5

6

7

8

Number of Channels

Proportion of False Positives

10-1

Proposed Method - Triangular Weights
Proposed Method - Square Weights
GLR

10-2

10-3

10

-4

10

-5

1

2

3

4

5

6

7

8

Number of Channels

Detection Delay (Samples Indices)

85

Proposed Method - Triangular Weights
Proposed Method - Square Weights
GLR

80
75
70
65
60
55
50

1

2

3

4

5

6

7

8

Number of Channels

Fig. 5: Figures comparing the proposed method (for both triangular and square filter weights) along with the GLR. (Upper
panel) Proportion of true positives, (middle panel) proportion of false positives and (lower panel) detection latency. Multivariate
data generated using random correlation matrices with low correlations.

The total number of samples was selected uniformly between [5000, 30000] and each segment length was selected with uniform
probability between the interval [300, 700]. Furthermore, the variance for each segment was dependant on the variance of the
prior segment; that is, a scale increase or decrease was selected with equal probability. The magnitude of the scale decrease
was selected with uniform probability between [0.5, 0.85] and similarly for a scale increase [1.2, 1.7].
We evaluate the performance of the proposed method when using different window parameters. Namely, we assess the
performance of the triangular window weights (selected from the analysis in Section III.B) shown in Fig. 3, along with the
standard square window weightings (that is a sliding window estimate of the standard deviation). In order to quantitatively
evaluate the performance of the respective algorithms, we utilised the following measures: the number of false positives and true
positives, the detection latency (the time between the event occurring and the corresponding first detection) and the absolute
error in the change point location estimate. The following parameters were selected for the respective algorithms: 1) GLR:
window size 250, threshold 5, 2) C/AFCD (for both the triangular and square filter parameters): Ts = 250, Tf = 20, γ = 0.8,
Td = 10, furthermore we selected the learning rate µ, such that the proportion of false positives for both the triangular and

12

Proportion of True Positives

1

Low Correlation
High Correlation

0.99
0.98
0.97
0.96
0.95
0.94
0.93
0.92

2

3

4

5

6

7

8

Number of Channels

Proportion of False Positives

0.012

Low Correlation
High Correlation

0.01
0.008
0.006
0.004
0.002
0
2

3

4

5

6

7

8

Number of Sensors

Detection Delay (Sample Indices)

60

Low Correlation
High Correlation

59
58
57
56
55
54
53
52
51

2

3

4

5

6

7

8

Number of Sensors

Fig. 6: Figures comparing the proposed method with triangular filter weights with multivariate data generated using random
correlation matrices with high and low correlations. (Upper panel) Proportion of true positives, (middle panel) proportion of
false positives and (lower panel) detection latency.

square filter weights were equal.
Fig. 5 shows the results of the proposed method (with the two different window parameters) along with the GLR method
when analysing univariate and multivariate data (where for this simulation we utilised randomly generated correlation matrices
with low absolute value in the off-diagonal entries). Observe that the proposed method (for both the triangular and square
window weights) outperformed the GLR with respect to the proportion of false positives and detection latency. Furthermore, it
can be seen that the proportion of true positives for the proposed method with triangular filter weights outperformed the GLR
when the number of sensors was greater than two; while the proposed method with square window weights performed poorly.
Finally, the results indicate that the weighting parameter for the volatility filters identified in section III.B outperformed the
standard square window volatility filter weights, with respect to the proportion of true positives.
We also assessed the performance of the proposed method with randomly generated correlation matrices with high absolute
values in the off-diagonal entries. In particular Fig. 6 compares the randomly generated high and low correlation matrices when
using volatility filters with triangular weights (developed in section III.B). Observe that the detection delay and proportion

13

30

Acceleration (m/s2)

Walk

Run

Walk

20

10

0

-10
0

1000

2000

3000
Time Index

4000

5000

6000

Fig. 7: Figure showing the time domain representation of the accelerometer data, where the three different states of motion in
time have been highlighted.
TABLE I: The performance of the CAFCD and GLR (applied to each channel) algorithms on real world accelerometer data.
False Negative False Negative
Proportion
Latency
CAFCD
GLR
CAFCD<GLR (No. of Samples)
0.927
0.8152
0.64
14.78

of false positives along with the true positives increase when there are large inter-channel correlations. This implies that for
highly correlated data, the proposed method is more susceptible to detecting changes, owing to the proportion of false positives
doubling, while the proportion of true positives does not increase significantly.
Finally, we compared the performance of the proposed change point location estimator with the GLR’s, where we evaluated
the performance of the respective algorithms using synthetic univariate time series data. In particular, we found that the
proposed VCE had an absolute error in the change point location estimate of 4.09 samples, while the GLR’s change point
location estimator had an error of 30.34 samples.

B. Accelerometer Data
The final real world simulation assessed the performance of the proposed algorithm on real world tri-axial accelerometer data
obtained from a smart phone attached to test subjects (where the data was obtained from [22] that consisted of 23 experimental
trials with sampling frequency of 128 Hz). Each experimental trial consists of five different states of motion: ‘no motion’,
‘walking’, ‘running’, ‘walking’ and ‘no motion’; where an example from the output of the accelerometer is shown in Fig. 7.
First differencing of the data was first carried out in order to remove trend components in the data. We then compared the
proposed CAFCD detection method (using the triangular filter weights scheme outlined in Section III.B) with the GLR method
that was applied seperately to each data channel. While the underlying true change point locations are unknown, the four
different changes in the volatility of the accelerometer data can be observed via visual inspection (see Fig. 7). As a result,
we can obtain the following quantitative performance measures: proportion of true positives, the proportion of the CAFCD
algorithm detecting the change in variance before GLR and the relative difference between the detection latencies of GLR and
CAFCD (that is a positive number would indicate CAFCD of having a smaller detection latency from the true change point
location). From Table I, it can be observed that the proposed change detection method had a lower number of false negatives
when compared with the GLR algorithm (applied to each channel). Furthermore, the proposed method was able to detect more

14

rapidly the transitions in 64% percent of the changes in the states of motion.
V. C ONCLUSION
This work developed a class of sequential volatility filtering based change detection and location estimation algorithms for
time series data with instantaneous changes in the volatility. The first method developed a change point detection method based
on adaptive filtering, and a convex combination of two volatility filters with differing window sizes. The adaptive filtering
formulation of the change point problem also enabled a straightforward extension to the analysis of changes in volatility in
multivariate time series data. This work also developed a change point location estimator based on the differenced output of a
volatility filter. The performance of the proposed detection and location estimation methods were evaluated on both synthetic
and real world data sets.
R EFERENCES
[1] M. Basseville and A. Benveniste, “A new statistical approach for the automatic segmentation of continuous speech signals,” IEEE Transactions on
Information Theory, vol. 29, no. 5, pp. 709–724, 1983.
[2] E. Sejdic, C. M. Steele, and T. Chau, “Segmentation of dual-axis swallowing accelerometry signals in healthy subjects with analysis of anthropometric
effects on duration of swallowing activities,” IEEE Transactions on Biomedical Engineering, vol. 56, no. 4, pp. 1090–1097, 2009.
[3] F. Desobry, M. Davy, and C. Doncarli, “An online kernel change detection algorithm,” IEEE Transactions on Signal Processing, vol. 53, no. 8, pp.
2961–2974, 2005.
[4] R. Andre-Obrecht, “A new statistical approach for the automatic segmentation of continuous speech signals,” IEEE Transactions on Acoustics, Speech,
and Signal Processing, vol. 36, no. 1, pp. 29–40, 1988.
[5] C. Inclan and G. C. Tiao, “Use of cumulative sums of squares for retrospective detection of changes of variance,” Journal of the American Statistical
Association, vol. 89, no. 427, pp. 913–923, 1994.
[6] E. Punskaya, C. Andrieu, A. Doucet, and W. Fitzgerald, “Bayesian curve fitting using MCMC with applications to signal segmentation,” IEEE Transactions
on Signal Processing, vol. 50, no. 3, pp. 747–758, 2002.
[7] Z. Wang and P. Willett, “Two algorithms to segment white Gaussian data with piecewise constant variances,” IEEE Transactions on Signal Processing,
vol. 51, no. 2, pp. 373–385, 2003.
[8] M.Lavielle, “Optimal segmentation of random processes,” IEEE Transactions on Signal Processing, vol. 46, no. 5, pp. 1365–1373, 1998.
[9] R. Killick, P. Fearnhead, and I. A. Eckley, “Optimal detection of changepoints with a linear computational cost,” Journal of the American Statistical
Association, vol. 500, no. 107, pp. 1590–1598, 2012.
[10] M. Lavielle, “Using penalized contrasts for the change-point problem,” Signal Processing, vol. 85, no. 8, pp. 1501–1510, 2005.
[11] X. Xuan and K. Murphy, “Modeling changing dependency structure in multivariate time series,” International Conference on Machine Learning, pp.
1055–1062, 2007.
[12] K. Bleakley and J.-P. Vert, “The group fused LASSO for multiple change-point detection,” arXiv:1106.4199, 2011.
[13] A. Brandt, “Detecting and estimating parameter jumps using ladder algorithms and likelihood ratio tests,” IEEE International Conference on Acoustics,
Speech and Signal Processing, pp. 1017–1020, 1983.
[14] C. Alexander, Moving average models for volatility and correlation, and covariance matrices, Handbook of Finance.

Wiley Press, 2008.

[15] Y. Zhang and J. A. Chambers, “Convex combination of adaptive filters for a variable tap-length LMS algorithm,” IEEE Signal Processing Letters, vol. 13,
no. 10, pp. 373–385, 2006.
[16] J. Arenas-Garcia, M. Martinez-Ramon, A. Navia-Vazquez, and A. R. Figueiras-Vidal, “Plant identification via adaptive combination of transversal filters,”
Signal Processing, vol. 86, no. 9, pp. 2430–2438, 2006.
[17] C. G. Lopes and A. H. Sayed, “Diffusion least-mean squares over adaptive networks: Formulation and performance analysis,” IEEE Transactions on
Signal Processing, vol. 56, no. 7, pp. 3122–3136, 2008.

15

[18] B. Jelfs and D. P. Mandic, “Signal modality characterisation using collaborative adaptive filters,” 1st IAPR Workshop on Cognitive Information Process,
2008.
[19] J. Arenas-Garcia, A. Figueiras-Vidal, and A. Sayed, “Mean-square performance of a convex combination of two adaptive filters,” IEEE Transactions on
Signal Processing, vol. 54, no. 3, pp. 1078–1090, 2006.
[20] R. Martin, W. Sethares, R. Williamson, and J. C. Johnson, “Exploiting sparsity in adaptive filters,” IEEE Transactions on Signal Processing, vol. 50,
no. 8, pp. 1883–1894, 2002.
[21] D. Lewandowski, D. Kurowicka, and H. Joe, “Generating random correlation matrices based on vines and extended onion method,” Journal of Multivariate
Analysis, vol. 100, no. 9, pp. 1989–2001, 2009.
[22] A. Brajdic and R. Harle, “Walk detection and step counting on unconstrained smartphones,” ACM International Joint Conference on Pervasive and
Ubiquitous Computing, p. 225234, 2013.

