Detecting Abrupt Changes in the Presence of Local Fluctuations and Autocorrelated Noise

Abstract

While there are a plethora of algorithms for detecting changes in mean in univariate time-series, almost all struggle in real applications where there is autocorrelated noise or where the mean fluctuates locally between the abrupt changes that one wishes to detect. In these cases, default implementations, which are often based on assumptions of a constant mean between changes and independent noise, can lead to substantial over-estimation of the number of changes. We propose a principled approach to detect such abrupt changes that models local fluctuations as a random walk process and autocorrelated noise via an AR(1) process. We then estimate the number and location of changepoints by minimizing a penalized cost based on this model. We develop a novel and efficient dynamic programming algorithm, DeCAFS, that can solve this minimization problem; despite the additional challenge of dependence across segments, due to the autocorrelated noise, which makes existing algorithms inapplicable. Theory and empirical results show that our approach has greater power at detecting abrupt changes than existing approaches. We apply our method to measuring gene expression levels in bacteria. Supplementary materials for this article are available online.

Keywords:

• Breakpoints
• Changepoints
• Dynamic programming
• FPOP
• Optimal partitioning
• Structural breaks

1 Introduction

Detecting changes in data streams is a ubiquitous challenge across many modern applications of statistics. It is important in such diverse areas as bioinformatics (Olshen et al. 2004; Futschik et al. 2014), ion channels (Hotz et al. 2013), climate records (Reeves et al. 2007), oceanographic data (Killick et al. 2010), and finance (Kim, Morley, and Nelson 2005). The most common and important change detection problem is that of detecting changes in mean, and there have been a large number of different approaches to this problem that have been proposed (e.g., Olshen et al. 2004; Killick, Fearnhead, and Eckley 2012; Fryzlewicz 2014; Frick, Munk, and Sieling 2014; Maidstone et al. 2017; Eichinger and Kirch 2018; Fearnhead and Rigaill 2019; Fryzlewicz 2018 b, amongst many others). Almost all of these methods are based on modeling the data as having a constant mean between changes and the noise in the data being independent. Furthermore, all change-point methods require specifying some threshold or penalty that affects the amount of evidence that there needs to be for a change before an additional changepoint is detected. In general, the methods have default choices of these thresholds or penalties that have good theoretical properties under strong modeling assumptions.

Although these methods perform well when analyzing simulated data where the assumptions of the method hold, they can be less reliable in real applications, particularly if the default threshold or penalties are used. Reasons for this include the noise in the data being autocorrelated, or the underlying mean fluctuating slightly between the abrupt changes that one wishes to detect. To see this, consider change detection for the well-log data (taken from Fearnhead and Liu 2011; Ruanaidh and Fitzgerald 2012) shown in Figure 1. These data come from lowering a probe into a bore-hole, and taking measurements of the rock structure as the probe is lowered. The data we plot has had outliers removed. As the probe moves from one rock strata to another, we expect to see an abrupt change in the signal from the measurements, and it is these changes that an analyst would wish to detect. Previous analyses of this data have shown that, marginally, the noise in the data is very well approximated by a Gaussian distribution; but by eye we can see local fluctuations in the data that suggest either autocorrelation in the measurement error, or structure in the mean between the abrupt changes.

Fig. 1 Segmentations of well-log data: wild binary segmentation using the strengthened Schwarz information criteria (top); segmentation under square error loss with penalty inflated to account for autocorrelation in measurement error (middle); optimal segmentation from DeCAFS with default penalty (bottom). Each plot shows the data (black line) the estimated mean (red line) and changepoint location (vertical blue dashed lines).

The top plot shows an analysis of the well-log data that uses wild binary segmentation (Fryzlewicz 2014) with the standard cusum test for a change in mean, and then estimates the number of changepoints based on a strengthened Schwarz information criteria. Both the cusum test and the strengthened Schwarz information criteria are based on modeling assumptions of a constant mean between changepoints and independent, identically distributed (IID) Gaussian noise, and are known to consistently estimate the number and location of the changepoints if these assumptions are correct. However, in this case, we can see that it massively overfits the number of changepoints. Similar results are obtained for standard implementation of other algorithms for detecting changes in mean, see Figure 13 in the supplementary material.

Lavielle and Moulines (2000) and Bardwell et al. (2019) suggested that if we estimate changepoints by minimizing the squared error loss of our fit with a penalty for each change, then we can correct for potential autocorrelation in the noise by inflating the penalty used for adding a changepoint. The middle plot of Figure 1 shows results for such an approach (Bardwell et al. 2019); this gives an improved result but it still noticeably overfits.

By comparison, the method we propose models both autocorrelation in the noise and local fluctuations in the mean between changepoints—and analysis of the data using default settings produces a much more reasonable segmentation of the data (see bottom plot of Figure 1). This method is model-based, and assumes that the local fluctuations in the mean are realizations of a random walk and that the noise process is an AR(1) process. We then segment the data by minimizing a penalized cost that is based on the log-likelihood of our model together with a BIC penalty for adding a changepoint.

The key algorithmic challenge with our approach is minimizing the penalized cost. In particular, many existing dynamic programming approaches (e.g., Jackson et al. 2005; Killick, Fearnhead, and Eckley 2012) do not work for our problem due to the dependence across segments caused by the autocorrelated noise. We introduce a novel extension of the functional pruned optimal partitioning algorithm of Maidstone et al. (2017), and we call the resulting algorithm DeCAFS, for Detecting Changes in Autocorrelated and Fluctuating Signals. It is both computationally efficient (analysis of the approx 4000 data points in the well-log data taking a fraction of a second on a standard laptop) and guaranteed to find the best segmentation under our criteria.

While we are unaware of any previous method that tries to model both autocorrelation and local fluctuations, Chakar et al. (2017) introduced AR1Seg which aims to detect changes in mean in the presence of autocorrelation. Their approach is similar to ours if we remove the random walk component, as they aim to minimize a penalized cost where the cost is the negative of the log-likelihood under a model with an AR(1) noise process. However they were unable to minimize this penalized cost, and instead minimized an approximation that removes the dependence across segments. One consequence of using this approximation is that it often estimates two consecutive changes at each changepoint, and AR1Seg uses a further post-processing step to try and correct this. Moreover, our simulation results show that using the approximation leads to a loss of power, particularly when the autocorrelation in the noise is high.

Particularly if interest lies in estimating how the underlying mean of the data varies over time, natural alternatives to DeCAFS are trend filtering methods (Kim et al. 2009; Tibshirani et al. 2014) which estimate the mean function under an L₁ penalty on a suitably chosen discrete derivative of the mean. Depending on the order of the derivative, these methods can fit a piecewise constant (in this case trend filtering is equivalent to the fused lasso of Tibshirani et al. 2005), linear, or quadratic etc. function to the mean. One advantage of trend filtering is its flexibility: depending on the application one can fit mean functions with different degrees of smoothness. However, it does not allow one to jointly estimate a smoothly changing mean and abrupt changes—the L₁ penalty must be chosen to capture one or other of these effects. This is particularly an issue if the main interest is in detecting abrupt changes rather than estimating the mean. These issues are investigated empirically in Appendix F.3.

Distinguishing between local fluctuations and abrupt changes is possible by methods that model the mean as a sum of functions (Jalali, Ravikumar, and Sanghavi 2013)—for example one smoothly varying and one piecewise constant. And trend-filtering, or other regularized estimators can then be used to estimate each component. The LAVA approach of Chernozhukov, Hansen, and Liao (2017) is one such approach, fitting the piecewise constant function and the smoothly varying function using, respectively, an L₁ and L₂ penalty on the function’s first discrete derivative. The main difference between LAVA and DeCAFS is thus that LAVA has an L₁ penalty for abrupt changes, and DeCAFS has an L₀ penalty. If interest is primarily in detecting changes, previous work has shown the use of an L₀ penalty to be preferable to an L₁ penalty – as the latter can often over-fit the number of changes (see, e.g., the empirical evidence in Fearnhead, Maidstone, and Letchford 2018; Jewell et al. 2020), and a post-processing step is often needed to correct for this (Lin et al. 2017; Safikhani and Shojaie 2020).

One important feature of modeling the random fluctuations in the mean via a random walk, is that the information that a data point y_s has about a change at time t decays to 0 as $|t-s|$ gets large. This is most easily seen if we consider applying DeCAFS to detect a single abrupt change. We show in Section 5 that whether DeCAFS detects a change at a time t depends on some contrast of the data before and after t. This contrast compares a weighted mean of the data before t to a weighted mean of the data after t, with the weights decaying (essentially) geometrically with the length of time before/after t. This is appropriate for the random walk model which enables, for example, the mean of $y_{t+h}$ to be very different to the mean at $y_{t+1}$ if h > 0 is large, and thus $y_{t+h}$ contains little to no information about whether there has been an abrupt change in the mean of $y_{t+1}$ compared to the mean of y_t . By contrast, standard CUSUM methods for detecting a change in mean quantify the evidence for a change at t by a contrast of the unweighted mean of the data before and after t – thus each data point, $y_{t+h}$, has the same amount of information about the change regardless of the value of h. In situations where there are local fluctuations in the mean but through some stationary process, such as a mean-reverting random walk, data far from a change would still have a non-negligble amount of information about it. In such situations, particularly if the segments are long, DeCAFS could have lower power than CUSUM or other methods that ignore any local fluctuations in the data. We investigate this empirically in Appendix F.4.

The outline of the article is as follows. In the next section, we introduce our model-based approach and the associated penalized cost. In Section 3 we present DeCAFS, a novel dynamic programming algorithm that can exactly minimize the penalized cost. To implement our method we need estimates of the model parameters, and we present a simple way of preprocessing the data to obtain these in Section 4. We then look at the theoretical properties of the method. These justify the use of the BIC penalty, show that our method has more power at detecting changes when our model assumptions are correct than standard approaches, and also that we have some robustness to model error—in that we can still consistently estimate the number and location of the changepoints in such cases by adapting the penalty for adding a changepoint. Sections 6 and 7 evaluate the new method on simulated and real data; and the article ends with a discussion.

Code implementing the new algorithm is available in the R package DeCAFS on CRAN. The package and full code from our simulation study is also available at github.com/gtromano/DeCAFS.

2 Modeling and Detecting Abrupt Changes

2.1 Model

Let $y_{1:n}=(y_1,\dots ,y_n)\in {\mathbb{R}}^n$ be a sequence of n observations, and assume we wish to detect abrupt changes in the mean of this data in the presence of local fluctuations and autocorrelated noise. We take a model-based approach where the signal vector is a realization of a random walk process with abrupt changes, and we super-impose an AR(1) noise process.

So for $t=1,\dots ,n$,(1) $$y_t={\mu }_t+{ϵ}_t,$$(1) where for $t=2,\dots ,n$ (2) $${\mu }_t={\mu }_{t-1}+{\eta }_t+{\delta }_t,\hspace{1em}\text{with}{\eta }_t\underset{\text{iid}}{\sim }\mathcal{N}(0,{\sigma }_{\eta }^2),{\delta }_t\in \mathbb{R},$$(2) and ${\delta }_t=0$ except at time points immediately after a set of m changepoints, $0<{\tau }_1<\cdots <{\tau }_m<n$. That is ${\delta }_t=0$ unless $t={\tau }_j+1$ for some j. This model is unidentifiable at changepoints. If τ is a changepoint, then whilst the data is informative about ${\mu }_{\tau +1}$ and ${\mu }_{\tau }$, we have no further information about the specific value of ${\delta }_{\tau +1}$ relative to ${\eta }_{\tau +1}$. We thus take the convention that ${\delta }_{\tau +1}={\mu }_{\tau +1}-{\mu }_{\tau }$ and ${\eta }_{\tau +1}=0$, which is the most likely value of ${\eta }_{\tau +1}$ under our model, and consistent with maximizing the likelihood criteria we introduce below. The noise process, ϵ_t is a stationary AR(1) process with, for $t=2,\dots ,n$,(3) $${ϵ}_t=\varphi {ϵ}_{t-1}+{\nu }_t\hspace{1em}\text{with}{\nu }_t\underset{\text{iid}}{\sim }\mathcal{N}(0,{\sigma }_{\nu }^2),$$(3) for some autocorrelation parameter, $\varphi$ with $\left|\varphi \right|<1$ and ${ϵ}_1\sim \mathcal{N}(0,{\sigma }_{\nu }^2/(1-{\varphi }^2)).$ Special cases of our model occur when $\varphi =0$ or when ${\sigma }_{\eta }^2=0$. When $\varphi =0$ our noise process ϵ_t is then IID, and the model is equivalent to a random walk plus noise with abrupt changes. When ${\sigma }_{\eta }^2=0$ we are detecting changes in mean with an AR(1) noise process, resulting in a formulation equivalent to the one of Chakar et al. (2017).

2.2 Penalized Maximum Likelihood Approach

In the following, we will assume that $\varphi ,{\sigma }_{\eta }^2$ and ${\sigma }_{\nu }^2$ are known; we consider robust approaches to estimate these parameters from the data in Section 4. We can then write down a type of likelihood for our model, defined as the joint density of the observations, $y_{1:n}$, and the local fluctuations in the mean, ${\eta }_{2:n}$. We will express this as a function of ${\mu }_{1:n}$ and ${\delta }_{2:n}$. Writing $f(·|·)$ for a generic conditional density, we have that this is$$\begin{array}{c}\mathcal{L}(y_{1:n};{\mu }_{1:n},{\delta }_{2:n})=\left(\prod _{t=2}^{n}f({\mu }_t|{\mu }_{t-1},{\delta }_t)\right)f(y_1|{\mu }_1)\\ \times \left(\prod _{t=2}^{n}f(y_t|y_{t-1},{\mu }_{t-1},{\mu }_t)\right)\\ \propto \left(\prod _{t=2}^n\hspace{0.17em}\text{exp}\hspace{0.17em}\left\{-\frac{{({\mu }_t-{\mu }_{t-1}-{\delta }_t)}^2}{2{\sigma }_{\eta }^2}\right\}\right)\hspace{0.17em}\text{exp}\hspace{0.17em}\left\{-\frac{{(y_1-{\mu }_1)}^2}{2{\sigma }_{\nu }^2/(1-{\varphi }^2)}\right\}\\ \times \left(\prod _{t=2}^n\hspace{0.17em}\text{exp}\hspace{0.17em}\left\{-\frac{{((y_t-{\mu }_t)-\varphi (y_{t-1}-{\mu }_{t-1}))}^2}{2{\sigma }_{\nu }^2}\right\}\right) .\end{array}$$

We have used the specific Gaussian densities of our model, and dropped multiplicative constants, to get the second expression.

If we knew the number of changepoints, then we could estimate their position by maximizing this likelihood subject to the constraints on the number of nonzero entries of ${\delta }_{2:n}$. However, as we need to also estimate the number of changepoints we proceed by maximizing a penalized version of the log of the likelihood where we introduce a penalty for each changepoint—this is a common approach to changepoint detection, see, for example, Maidstone et al. (2017). It is customary to restate this as minimizing a penalized cost, rather than maximizing a penalized likelihood, where the cost is minus twice the log-likelihood. That is, we estimate the number and location of the changepoints by solving the following minimization problem:(4) $$\begin{array}{c}{\mathcal{F}}_n={\min }_{\stackrel{{\mu }_{1:n}}{{\delta }_{2:n}}}\{(1-{\varphi }^2)\gamma {(y_1-{\mu }_1)}^2 +\\ \sum _{t=2}^n[\lambda {({\mu }_t-{\mu }_{t-1}-{\delta }_t)}^2+\gamma {((y_t-{\mu }_t)-\varphi (y_{t-1}-{\mu }_{t-1}))}^2\\ +\beta 1_{{\delta }_t\ne 0}\left]\right\},\end{array}$$(4) where $\beta >0$ is the penalty for adding a changepoint, $\lambda =1/{\sigma }_{\eta }^2,\gamma =1/{\sigma }_{\nu }^2$, and $1\in \{0,1\}$ is an indicator function. For the special case of a constant mean between changepoints, corresponding to ${\sigma }_{\eta }^2=0$, we require ${\mu }_t={\mu }_{t-1}+{\delta }_t\forall t=2,\dots ,n$ and simply drop the first term in the sum.

2.3 Dynamic Programming Recursion

We will use dynamic programming to minimize the penalized cost Equation (4)(4) $$\begin{array}{c}{\mathcal{F}}_n={\min }_{\stackrel{{\mu }_{1:n}}{{\delta }_{2:n}}}\{(1-{\varphi }^2)\gamma {(y_1-{\mu }_1)}^2 +\\ \sum _{t=2}^n[\lambda {({\mu }_t-{\mu }_{t-1}-{\delta }_t)}^2+\gamma {((y_t-{\mu }_t)-\varphi (y_{t-1}-{\mu }_{t-1}))}^2\\ +\beta 1_{{\delta }_t\ne 0}\left]\right\},\end{array}$$(4) . The challenge here is to deal with the dependence across changepoints due to the AR(1) noise process which means that some standard dynamic approaches for change-point detection, such as optimal partitioning (Jackson et al. 2005) and PELT (Killick, Fearnhead, and Eckley 2012), cannot be used. To overcome this, as in Rigaill (2015) or Maidstone et al. (2017), we define the function $\mu \mapsto Q_t(\mu )$ to be the minimum penalized cost for data $y_{1:t}$ conditional on ${\mu }_t=\mu$,$$\begin{array}{c}{Q}_t(\mu )={\min }_{\stackrel{{\mu }_{1:t}}{{\delta }_{2:t},{\mu }_t=\mu }}\{(1-{\varphi }^2)\gamma {(y_1-{\mu }_1)}^2 +\\ \sum _{i=2}^t[\lambda {({\mu }_i-{\mu }_{i-1}-{\delta }_i)}^2\\ +\gamma {((y_i-{\mu }_i)-\varphi (y_{i-1}-{\mu }_{i-1}))}^2+\beta 1_{{\delta }_i\ne 0}\left]\right\}.\end{array}$$

So ${\mathcal{F}}_n={\min }_{\mu \in \mathbb{R}}{Q}_n(\mu )$; and the following proposition gives a recursion for $Q_t(\mu )$.

Proposition 1.

The set of functions $\{\mu \mapsto Q_t(\mu ) , t=1,\dots ,n\}$ satisfies$$Q_1(\mu )=(1-{\varphi }^2)\gamma {(y_1-\mu )}^2\text{and},\mathrm{ }\text{for}\mathrm{ }t=2,\dots ,n,$$(5) $$\begin{array}{c}{Q}_t(\mu )={\min }_{u\in \mathbb{R}}\{Q_{t-1}(u)+\min \{\lambda {(\mu -u)}^2,\beta \}\\ \hspace{1em}+\gamma {((y_t-\mu )-\varphi (y_{t-1}-u))}^2\} .\end{array}$$(5)

The intuition behind the recursion is that we first condition on ${\mu }_{t-1}=u$, with the term in braces being the minimum penalized cost for $y_{1:t}$ given u and ${\mu }_t=\mu$, and then minimize over u. The cost in braces is the sum of three terms: (i) the minimum penalized cost for $y_{1:t-1}$ given u; (ii) the cost for the change in mean from u to μ; and (iii) the cost of fitting data point y_t with μ_t . The cost for the change in mean, (ii), is just the minimum of the constant cost for adding a change and the quadratic cost for a change due to the random walk. The recursion applies to the special case of a constant mean between changepoints, where $\lambda =\infty$, if we replace $\min \{\lambda {(\mu -u)}^2,\beta \}$ with its limit as $\lambda \to \infty$, which is $\beta 1_{\mu \ne u}$.

Although recursions of this form have been considered in earlier changepoint algorithms (e.g., Maidstone et al. 2017; Hocking et al. 2020), the dependence between the current mean u and the previous mean μ that appears in terms (ii) and (iii) makes our recursion more challenging to solve. We next show one efficient way of solving by combining existing functional pruning dynamic programming ideas with properties of infimal convolutions.

3 Computationally Efficient Algorithm

3.1 The DeCAFS Algorithm

Algorithm 1 gives pseudo code for solving the dynamic programming recursion introduced in Proposition 1. The key to implementing this algorithm is performing the calculations in line 5, and how this can be done efficiently will be described below. Throughout we give the algorithm for the case where there is a random walk component, that is, $\lambda <\infty$, though it is trivial to adapt the algorithm to the $\lambda =\infty$ case.

As well as solving the recursion for $Q_t(\mu )$, Algorithm 1 shows how we can also obtain the estimate of the mean, through a standard back-tracking step. The idea is that our estimate of μ_n , ${\widehat{\mu }}_n$, is just the value of μ that maximizes $Q_n(\mu )$. We then loop backwards through the data, and our estimate of μ_t is the value that minimizes the penalized cost for the data $y_{1:t}$ conditional on ${\mu }_{t+1}={\widehat{\mu }}_{t+1}$, which can be calculated as $B_t(\mu )$ in line 11.

Finally, as we obtain the estimates of the mean, we can also directly obtain the estimated changepoint locations. It is straightforward to see, by examining the form of the penalized cost, that the optimal solution for ${\delta }_{2:n}$ has ${\delta }_{t+1}\ne 0$ (and hence, t is a changepoint) if and only if $\lambda {({\widehat{\mu }}_{t+1}-{\widehat{\mu }}_t)}^2>\beta$.

Algorithm 1: DeCAFS

Data: $\mathbf{y}=y_{1:n}$ a time series of length n

Input: $\beta >0,\lambda >0,\gamma >0$ and $0\le \varphi <1$.

1 begin Initialization

2  $Q_1(\mu )\stackrel{}{\leftarrow }(1-{\varphi }^2)\gamma {(y_1-\mu )}^2$

3 end

4 for t = 2 to n do

5  $Q_t(\mu )\stackrel{}{\leftarrow }\underset{u}{\min }\{Q_{t-1}(u)+\min \{\lambda {(\mu -u)}^2,\beta \}+$

$\gamma {((y_t-\mu )-\varphi (y_{t-1}-u))}^2\}$

6 end

7 begin Backtracking

8  ${\widehat{\mu }}_n\stackrel{}{\leftarrow }\text{argmin}{Q}_n(\mu )$

9  $\widehat{\tau }\stackrel{}{\leftarrow }n$

10  for $t=n-1$ to 1 do

11   $B_t(\mu )\stackrel{}{\leftarrow }{Q}_t(\mu )+\min \{\lambda {(\mu -{\widehat{\mu }}_{t+1})}^2,\beta \}+$

 $\gamma {((y_{t+1}-{\widehat{\mu }}_{t+1})-\varphi (y_t-\mu ))}^2$

12   ${\widehat{\mu }}_t\stackrel{}{\leftarrow }\text{argmin}{B}_t(\mu )$

13   if ${({\widehat{\mu }}_t-{\widehat{\mu }}_{t+1})}^2>\beta /\lambda$ then

14    $\widehat{\tau }\stackrel{}{\leftarrow }(t,\widehat{\tau })$

15   end

16  end

17 end

18 Return ${\widehat{\mu }}_{1:n},\widehat{\tau }$

3.2 The Infimal Convolution

The main challenge with Algorithm 1 is implementing line 5. First, this needs a compact way of characterizing $Q_t(\mu )$. This is possible as $Q_1(\mu )$ is a quadratic function; and the recursion maps piecewise quadratic functions to piecewise quadratic functions. Hence $Q_t(\mu )$ will be piecewise quadratic and can be defined by storing a partition of the real-line together with the coefficients of the quadratics for each interval in this partition.

Next we can simplify line 5 of Algorithm 1. As written line 5 involves minimizing a two-dimensional function, in $(u,\mu )\in {\mathbb{R}}^2$, over the variable u. We can recast this operation into a one-dimensional problem by introducing the concept of an infimal convolution (see Chapter 12 of Bauschke and Combettes 2011).

Definition 1.

Let f be a real-valued function defined on $\mathbb{R}$ and ω a nonnegative scalar. We define $IN{F}_{f,\infty }(\theta )=f(\theta )$ and for $\omega >0$,(6) $$IN{F}_{f,\omega }(\theta )={\min }_{u\in \mathbb{R}}\left(f(u)+\omega {(u-\theta )}^2\right),$$(6) as the infimal convolution of f with a quadratic term.

The following proposition presents a reformulation of the update-rule into a minimization involving infimal convolutions, for the case $\varphi \ge 0$. The proof is in Appendix B, together with details of equivalent results when $\varphi <0$.

Proposition 2.

Assume $\varphi \ge 0$. The functions $\{Q_t(\mu ) , t=2,\dots ,n\}$ can be written as$$Q_t(\mu )=\min \left\{Q_t^{=}(\mu ),Q_t^{\ne }(\mu )\right\},$$where$$\begin{array}{c}{Q}_t^{=}(\mu )=IN{F}_{{\mathbb{Q}}_{t-1},\gamma \varphi +\lambda }(\mu )\\ \hspace{1em}+\frac{\gamma }{1-\varphi }{(y_t-\varphi y_{t-1}-(1-\varphi )\mu )}^2 ,\\ Q_t^{\ne }(\mu )=IN{F}_{{\mathbb{Q}}_{t-1},\gamma \varphi }(\mu )\\ \hspace{1em}+\frac{\gamma }{1-\varphi }{(y_t-\varphi y_{t-1}-(1-\varphi )\mu )}^2+\beta ,\end{array}$$ and$${\mathbb{Q}}_{t-1}(u)=Q_{t-1}(u)-\gamma \varphi (1-\varphi ){\left(u-\frac{y_t-\varphi y_{t-1}}{1-\varphi }\right)}^2.$$

3.3 Fast Infimal Convolution Computation

As noted above we can represent Q_t by $Q_t=(q_t^1,\dots ,q_t^s)$ where each $q_t^i$ is a quadratic defined on some interval $[d_i,d_{i+1}[$ with $d_1=-\infty$ and $d_{s+1}=+\infty$. It is this representation of Q_t that we update at each time step. Some operations involved in solving the recursion, such as adding a quadratic to a piecewise quadratic, or calculating the pointwise minimum of two piecewise quadratics are easy to perform with a computational cost that is linear in the number of intervals (see, e.g., Rigaill 2015). The following theorem shows that a fast update for the infimal convolution of a piecewise quadratic is also possible, and is important for developing a fast algorithm for solving the dynamic programming recursions.

Theorem 1.

Let $Q_t=(q_t^1,\dots ,q_t^s)$ be the representation of the functional cost Q_t . For all $\omega \ge 0$, the representation returned by the infimal convolution ${\text{INF}}_{Q_t,\omega }$ has the following order-preserving form:$${\text{INF}}_{Q_t,\omega }=(\text{INF}{q}_t^{u_1},\text{INF}{q}_t^{u_2},\dots ,\text{INF}{q}_t^{u_{s^{*}-1}},\text{INF}{q}_t^{u_{s^{*}}}) ,$$with $1=u_1<u_2<\dots <u_{s^{*}-1}<u_{s^{*}}=s$ and $s^{*}\le s$.

The key part of this result is that the order of the quadratics is not changed when we apply the infimal convolution, and thus we can calculate ${\text{INF}}_{Q_t,\omega }$ using a linear scan over the real-line. The proof of this theorem is given in Appendix C, and an example algorithm for calculating ${\text{INF}}_{Q_t,\omega }$ with complexity that is linear in s is shown in Appendix D.

4 Robust Parameter Estimation

Our optimization problem Equation (4)(4) $$\begin{array}{c}{\mathcal{F}}_n={\min }_{\stackrel{{\mu }_{1:n}}{{\delta }_{2:n}}}\{(1-{\varphi }^2)\gamma {(y_1-{\mu }_1)}^2 +\\ \sum _{t=2}^n[\lambda {({\mu }_t-{\mu }_{t-1}-{\delta }_t)}^2+\gamma {((y_t-{\mu }_t)-\varphi (y_{t-1}-{\mu }_{t-1}))}^2\\ +\beta 1_{{\delta }_t\ne 0}\left]\right\},\end{array}$$(4) depends on three unknown parameters: ${\sigma }_{\eta }^2,{\sigma }_{\nu }^2$ and $\varphi$. We estimate these parameters by fitting to robust estimates of the variance of the k-lag differenced data, $z_t^k=y_{t+k}-y_t$Z$\mathrm{(}\text{HTML}\mathrm{ }\text{translation}\mathrm{ }\text{failed}\mathrm{)}$, for $k\ge 1$.

Proposition 3.

With the model defined by Equations (1)–(3),$$z_t^k\sim \mathcal{N}(\sum _{i=t+1}^{t+k}{\delta }_i,k{\sigma }_{\eta }^2+2\frac{1-{\varphi }^k}{1-{\varphi }^2}{\sigma }_{\nu }^2) ,\hspace{1em}t=1,\dots ,n-k.$$

Providing k is small relative to the average length of a segment, the mean of $z_t^k$ will be zero for most t: if there are m changes then at most km of the $z_t^k$’s will overlap a change and have a nonzero mean. This suggests that we can estimate the variance of $z_t^k$ using a robust estimator, such as the median absolute difference from the median, or MAD, estimator.

Fix K, and let v_k be the MAD estimator of the variance of $z_t^k$ for $k=1,\dots ,K$. We estimate the parameters by minimizing the least-square fit to these estimates,$${\mathcal{S}}_{\varphi }({\sigma }_{\eta }^2,{\sigma }_{\nu }^2)=\sum _{k=1}^K(k{\sigma }_{\eta }^2+2\frac{1-{\varphi }^k}{1-{\varphi }^2}{\sigma }_{\nu }^2-v_k{)}^2 .$$

In practice, we can minimize these criteria by using a grid of values for $\varphi$ and then for each $\varphi$ value analytically minimize with respect to ${\sigma }_{\eta }^2\ge 0$ and ${\sigma }_{\nu }^2\ge 0$. Obviously, if we are fitting a model without the random walk component, then we can set ${\sigma }_{\eta }^2=0$, or if we wish to have uncorrelated noise, then we set $\varphi =0$.

An empirical evaluation of this method for estimating the parameters is shown in Section F of the supplementary material. These include an investigation of the accuracy in situations where we have changepoints. In our simulation study, we use K = 10, though similar results were obtained as we varied K.

5 Theoretical Properties

We can reformulate our model as a linear-regression. To do this it is helpful to introduce new variables, ${\tilde{\eta }}_{1:n}$, that give the cumulative effect of the random-walk fluctuations. To simplify exposition, it is further helpful to define this process so it has an invertible covariance matrix. So we will let ${\tilde{\eta }}_1\sim \mathcal{N}(0,{\sigma }_{\eta }^2)$ and ${\tilde{\eta }}_t={\tilde{\eta }}_{t-1}+{\eta }_t$ for $t=2,\dots ,n$. For a set of m changepoints ${\tau }_{1:m}$, and defining ${\tau }_0=0$, we can introduce a $n\times (m+1)$ matrix $X_{{\tau }_{0:m}}$ where the ith column is a column of ${\tau }_{i-1}$ zeros followed by $n-{\tau }_{i-1}$ ones. Our model is then(7) $$y_{1:n}=X_{{\tau }_{0:m}}\Delta +{\zeta }_{1:n},$$(7) where ${\zeta }_{1:n}$ is a vector of Gaussian random variables with$$\text{var}({\zeta }_{1:n})=\text{var}({ϵ}_{1:n})+\text{var}({\tilde{\eta }}_{1:n}):={\Sigma }_{\text{AR}}+{\Sigma }_{\text{RW}}$$ the sum of the variance matrices for the AR component of the model, ${ϵ}_{1:n}$, and the random walk component of the model, ${\tilde{\eta }}_{1:n}$; and Δ is a $(m+1)\times 1$ vector whose first entry is ${\mu }_1-{\tilde{\eta }}_1$ and whose ith entry is ${\delta }_{{\tau }_{i-1}+1}$ the change at the $(i-1)$th changepoint. This formulation allows us to consider the impact of model error on DeCAFS. Later, when we consider its asymptotic properties, we will allow for the data-generating process to be Equation (7)(7) $$y_{1:n}=X_{{\tau }_{0:m}}\Delta +{\zeta }_{1:n},$$(7) but with $\text{var}({\zeta }_{1:n})$ different from that assumed by DeCAFS.

As shown in Section E of the supplementary material, the unpenalized version of the cost that we minimize, conditional on a specific set of changepoints, can be written as$$\begin{array}{c}\mathcal{C}({\tau }_{1:m})={\min }_{\Delta ,{\tilde{\eta }}_{1:n},{\tilde{\eta }}_1=0}[{(y_{1:n}-X_{{\tau }_{0:m}}\Delta -{\tilde{\eta }}_{1:n})}^T\\ \hspace{1em}\times {\Sigma }_{\text{AR}}^{-1}(y_{1:n}-X_{{\tau }_{0:m}}\Delta -{\tilde{\eta }}_{1:n})+{\tilde{\eta }}_{1:n}^T{\Sigma }_{\text{RW}}^{-1}{\tilde{\eta }}_{1:n}],\end{array}$$where ${\tilde{\eta }}_{1:n}$ is assumed to be a column vector. Thus, the penalized cost Equation (4)(4) $$\begin{array}{c}{\mathcal{F}}_n={\min }_{\stackrel{{\mu }_{1:n}}{{\delta }_{2:n}}}\{(1-{\varphi }^2)\gamma {(y_1-{\mu }_1)}^2 +\\ \sum _{t=2}^n[\lambda {({\mu }_t-{\mu }_{t-1}-{\delta }_t)}^2+\gamma {((y_t-{\mu }_t)-\varphi (y_{t-1}-{\mu }_{t-1}))}^2\\ +\beta 1_{{\delta }_t\ne 0}\left]\right\},\end{array}$$(4) is ${\mathcal{F}}_n={\min }_{m,{\tau }_{1:m}}\left[\mathcal{C}({\tau }_{1:m})+m\beta \right]$. In the remainder of this section, we will call $\mathcal{C}({\tau }_{1:m})$ the cost, and $\mathcal{C}({\tau }_{1:m})+m\beta$ the penalized cost.

While our cost is obtained by minimizing over ${\eta }_{2:n}$, the following result shows that it is equal to the weighted residual sum of squares from fitting the linear model defined in Equation (7)(7) $$y_{1:n}=X_{{\tau }_{0:m}}\Delta +{\zeta }_{1:n},$$(7) .

Proposition 4.

The cost for fitting a model with changepoints, ${\tau }_{1:m}$ is(8) $$\mathcal{C}({\tau }_{1:m})={\min }_{\Delta }{(y_{1:n}-X_{{\tau }_{0:m}}\Delta )}^T{\left({\Sigma }_{\text{AR}}+{\Sigma }_{\text{RW}}\right)}^{-1}(y_{1:n}-X_{{\tau }_{0:m}}\Delta )$$(8)

Let ${\mathcal{C}}_0$ denote the cost if we fit a model with no changepoints. The following corollary, which follows from standard arguments, gives the behavior of the cost under a null model of no changepoints. This includes a bound on the impact of mis-specifying the covariance matrix, for example due to mis-estimating the parameters of the AR(1) or random walk components of the model, or if our model for the residuals is incorrect.

Corollary 1.

Assume that data is generated from the model defined in Equation (7)(7) $$y_{1:n}=X_{{\tau }_{0:m}}\Delta +{\zeta }_{1:n},$$(7) , with m = 0 but with ${\zeta }_{1:n}$ a mean-zero Gaussian vector with $\text{var}({\zeta }_{1:n})=\Sigma$. Let ${\alpha }_n^{+}$ be the largest eigenvalue of ${({\Sigma }_{\text{AR}}+{\Sigma }_{\text{RW}})}^{-1}\Sigma$. If $\Sigma ={\Sigma }_{\text{AR}}+{\Sigma }_{\text{RW}}$ then ${\mathcal{C}}_0-\mathcal{C}({\tau }_{1:d})\sim {\chi }_d^2$. Otherwise, for any x $$\mathrm{Pr}({\mathcal{C}}_0-\mathcal{C}({\tau }_{1:d})>x)\le \mathrm{Pr}({\chi }_d^2>x/{\alpha }_n^{+}).$$

Furthermore, if we estimate the number of changepoints using the penalized cost Equation (4)(4) $$\begin{array}{c}{\mathcal{F}}_n={\min }_{\stackrel{{\mu }_{1:n}}{{\delta }_{2:n}}}\{(1-{\varphi }^2)\gamma {(y_1-{\mu }_1)}^2 +\\ \sum _{t=2}^n[\lambda {({\mu }_t-{\mu }_{t-1}-{\delta }_t)}^2+\gamma {((y_t-{\mu }_t)-\varphi (y_{t-1}-{\mu }_{t-1}))}^2\\ +\beta 1_{{\delta }_t\ne 0}\left]\right\},\end{array}$$(4) with penalty $\beta =C{\alpha }_n^{+}\hspace{0.17em}\text{log}\hspace{0.17em}n$ for any C > 2, then the estimated number of changepoints, $\widehat{m}$, satisfies $\mathrm{Pr}(\widehat{m}=0)\to 1$ as $n\to \infty$.

To gain insight into the behavior of the procedure in the presence of changepoints, and how it differs from standard standard change-in-mean procedures, it is helpful to consider the reduction in cost if we add a single changepoint.

Proposition 5.

Given a fixed changepoint location τ₁:

1. The reduction in cost for adding a single changepoint at τ₁ can be written as ${\mathcal{C}}_0-\mathcal{C}({\tau }_1)={(v^T{y}_{1:n})}^2$, for some vector v defined as$$\begin{array}{c}v=\frac{1}{\sqrt{c_{{\tau }_1}-c_{0,{\tau }_1}^2/c_0}}\{{({\Sigma }_{\text{AR}}+{\Sigma }_{\text{RW}})}^{-1}{u}_{{\tau }_1}\\ \hspace{1em}-\frac{c_{0,{\tau }_1}}{c_0}{({\Sigma }_{\text{AR}}+{\Sigma }_{\text{RW}})}^{-1}{u}_0\},\end{array}$$

   where u₀ is a column vector of n ones, $u_{{\tau }_1}$ is a column vector of τ₁ zeroes followed by $n-{\tau }_0$ ones, and$$\begin{array}{c}{c}_0\mathrm{}=\mathrm{}{u}_0^T{({\Sigma }_{\text{AR}}\mathrm{}+\mathrm{}{\Sigma }_{\text{RW}})}^{-1}{u}_0,c_{0,{\tau }_1}\mathrm{}=\mathrm{}{u}_0^T{({\Sigma }_{\text{AR}}\mathrm{}+\mathrm{}{\Sigma }_{\text{RW}})}^{-1}{u}_{{\tau }_1}, \\ c_{{\tau }_1}=u_{{\tau }_1}^T{({\Sigma }_{\text{AR}}+{\Sigma }_{\text{RW}})}^{-1}{u}_{{\tau }_1}.\end{array}$$

2. The vector v in (i) satisfies ${\sum }_{i=1}^n{v}_i=0$ and $v^T({\Sigma }_{\text{AR}}+{\Sigma }_{\text{RW}})v=1$.

3. For any vector w that satisfies ${\sum }_{i=1}^n{w}_i=0$ and $w^T({\Sigma }_{\text{AR}}+{\Sigma }_{\text{RW}})w=1$,$${\left(\sum _{i={\tau }_1+1}^n{w}_i\right)}^2\le {\left(\sum _{i={\tau }_1+1}^n{v}_i\right)}^2.$$

The vector v in part (i) of this proposition defines a projection of the data that is used to determine whether to add a changepoint at τ₁. The properties in part (ii) mean that this projection is invariant to shifts of the data, and that the distribution of the reduction in cost if our model is correct and there are no changes will be ${\chi }_1^2$. The statistic $v^T{y}_{1:n}$ can be viewed as analogous to the cusum statistic (Hinkley 1971) that is often used for a standard change-in-mean problem, and in fact if we set $\varphi =0$ and ${\sigma }_{\eta }=0$ so as to remove the auto regressive and random-walk aspects of the model, $\left|v^T{y}_{1:n}\right|$ is just the standard cusum statistic. The power of our method to detect a change at τ₁ will be governed by the distribution of this projection applied to the data in the segments immediately before and after τ₁. For a single changepoint where the mean changes by δ this distribution is a noncentral chi-squared with 1 degree of freedom and noncentrality parameter ${\delta }^2{({\sum }_{i={\tau }_1+1}^n{v}_i)}^2$. Thus, part (iii) shows that v is the best linear projection, in terms of maximizing the non-centrality parameter, over all projections that are invariant to shifts in the data and that are scaled so that the null distribution is ${\delta }^2{({\sum }_{i={\tau }_1+1}^n{v}_i)}^2$.

To gain insight into how the auto regressive and random-walk parts of the model affect the information in the data about a change we have plotted different projections v for different model scenarios in the top row of Figure 2. The top-left plot shows the projections if we have $\varphi =0$ for different, values of the random walk variance. The projection, naturally, places more weight to data near the putative changepoint, and the weight decays essentially geometrically as we move away from the putative changepoint. In the top-right plot we show the impact of increasing the autocorrelation of the AR(1) process, with the absolute value of the weight given to data points immediately before and after the putative change increasing with $\varphi$.

Fig. 2 Top row: projections of data v for detecting a change in the middle of n = 100 data-points. Random walk model (top-left) for varying ${\sigma }_{\eta }^2$ of 0.03 (black), 0.02 (red) and 0.01 (green); AR(1) plus random walk model (top-right) for ${\sigma }_{\eta }^2=0.01$ and varying $\varphi$ of 0.4 (black), 0.2 (red) and 0.1 (green). In both plots the blue line shows the standard cusum projection. Bottom row: noncentrality parameter for a ${\chi }_1^2$ test of a change using the optimal projection (solid line) and the cusum projection (dashed line) for a change of size 1 in the middle of the data as we vary n. Out-fill asymptotics (bottom-left) where $({\sigma }_{\eta }^2,\varphi )$ is (0.0025,0) (black), (0.01,0) (red), (0.0025,0.5) (green) and (0.01,0.5) (blue); In-fill asymptotics (bottom-right) where for n = 50 $({\sigma }_{\eta }^2,\varphi )$ is (0.0025,0) (black), (0.01,0) (red), (0.0025,0.5) (green) and (0.01,0.5) (blue).

A key feature of the random walk model is that for any fixed ${\sigma }_{\eta }^2>0$ the amount of information about a change will be bounded as we increase the segment lengths either side of the change. This is shown in the bottom-left plot of Figure 2 where we show the noncentrality parameter for detecting a change in the middle of the data as we vary n. For comparison we also show the non-centrality parameter of a test based on the cusum statistic (scaled so that it also has a ${\chi }_1^2$ distribution under the null of no change). We can see that ignoring local fluctuations in the mean, if they exist and come from a random walk model, by using the cusum statistic leads to a reduction of power as segment lengths increase. For comparison in the bottom right we show an equivalent comparison where we consider an infill asymptotic regime, so that as n increases we let the random walk variance decay at a rate proportional to $1/n$ and we increase the lag-1 autocorrelation appropriately. In this case using the optimal projection gives a noncentrality parameter that increases with n, whereas the cusum statistic has power that can be shown to be bounded as we increase n.

We now turn to the property of our method at detecting multiple changes. Based on the above discussion, we will consider in-fill asymptotics as $n\to \infty$.

• (C1) Let $y_1,\dots ,y_n$ be generated as a finite sample from a Gaussian process on $[0,1]$; that is $y_i=z(i/n)$ where, for $t\in [0,1]$ $z(t)=\mu (t)+\zeta (t),\mu (t)$ is a piecewise constant with m⁰ changepoints at locations $r_1,\dots ,r_{m^0}$, and $\zeta (t)$ is a mean zero Gaussian process. For a given n define the true changepoint locations as ${\tau }_i^0=\lfloor n{r}_i^0\rfloor$. The change in mean at each changepoint is fixed and non-zero.

• (C2) Assume there exists strictly positive constants $c_{\eta },c_{\nu }$ and $c_{\varphi }$, such that we implement DeCAFS with ${\sigma }_{\eta }^2=c_{\eta }/n$ and either (i) $\varphi =0$ and ${\sigma }_{\nu }^2=c_{\nu }$; or (ii) $\varphi =\hspace{0.17em}\text{exp}\hspace{0.17em}\{-c_{\varphi }/n\}$ and ${\sigma }_{\nu }^2=c_{\nu }(1-\hspace{0.17em}\text{exp}\hspace{0.17em}\{-2c_{\varphi }/n\})$.

• (C3) There exists an α such that for any large enough n if ${\Sigma }_n^0$ is the covariance of the noise in the data generating model (C1), and ${\Sigma }_{\text{AR}}^{(n)}+{\Sigma }_{\text{RW}}^{(n)}$ is the covariance assumed by DeCAFS in (C2) then the largest eigenvalue of ${({\Sigma }_{\text{AR}}^{(n)}+{\Sigma }_{\text{RW}}^{(n)})}^{-1}{\Sigma }_n^0$ is less than α.

The two regimes covered by condition C2 are due to the different limiting behavior of an AR(1) model under in-fill asymptotics, depending on whether the AR(1) noise is independent, case (i), or there is autocorrelation, case (ii). The form of ${\sigma }_{\nu }^2$ in each case ensures that the AR(1) process has fixed marginal variance, $c_{\nu }$, for all values of n.

The key condition here is (C3) which governs how accurate the model assumed by DeCAFS is to the true data-generating procedure. Clearly, if the model is correct then (C3) holds with α = 1. The following proposition gives upper bound on α in the case where the covariance of the data-generating model is that of a random walk plus AR(1) process, but with different parameter values to those assumed by DeCAFS in (C2), for example, due to misestimation of these parameters.

Proposition 6.

Assume the noise process $\zeta (t)$ of the data-generating process (C1) is equal to a random walk plus an AR(1) process.

1. If $\text{cov}(\zeta (t),\zeta (s))=c_{\eta }^0\min (t,s)$ for $t\ne s$ and $\text{var}(\zeta (t))=c_{\eta }^{0}t+c_{\nu }$, and DeCAFS is implemented as in (C2)(i), then (C3) holds with $\alpha =\max \{c_{\nu }^0/c_{\nu },c_{\eta }^0/c_{\eta }\}$.

2. If $\text{cov}(\zeta (t),\zeta (s))=c_{\eta }^0\min (t,s)+c_{\nu }^0\hspace{0.17em}\text{exp}\hspace{0.17em}\{-c_{\varphi }^0|t-s\left|\right\}$ and DeCAFS is implemented as in (C2)(ii), then for any $ϵ>0$ (C3) holds with$$\alpha =\max \left\{\frac{c_{\nu }^0}{c_{\nu }}\frac{c_{\varphi }^0}{c_{\varphi }}(1+ϵ),\frac{c_{\nu }^0}{c_{\nu }}\left(1+\frac{c_{\varphi }}{c_{\varphi }^0}\right)(1+ϵ),\frac{c_{\eta }^0}{c_{\eta }}\right\}.$$

The following result shows that we can consistently estimate the number of changepoints and gives a bound on the error in the estimate of changepoint locations, if we use DeCAFS under an assumption of a maximum number of changepoints (the assumption of a maximum number changes is for technical convenience, though is common in similar results, e.g., Yao 1988).

Theorem 2.

Assume data, $y_{1:n}$, is generated as described in (C1), and let $\widehat{m}$ and ${\widehat{\tau }}_{1:\widehat{m}}$ be the estimated number and location of the changepoints from DeCAFS implemented with parameters given by (C2), penalty $\beta =C\alpha \hspace{0.17em}\text{log}\hspace{0.17em}n$ for some C > 2, and a maximum number of changes $m_{\max }\ge m^0$. Then as $n\to \infty$: if $\varphi >0$ $$\mathrm{Pr}\left(\widehat{m}=m^0,{\max }_{i=1,\dots ,m^0}\left|{\widehat{\tau }}_i-{\tau }_i^0\right|=0\right)\to 1;$$and if $\varphi =0$ $$\mathrm{Pr}\left(\widehat{m}=m^0,{\max }_{i=1,\dots ,m^0}\left|{\widehat{\tau }}_i-{\tau }_i^0\right|\le {(\hspace{0.17em}\text{log}\hspace{0.17em}n)}^2\right)\to 1.$$

The most striking part of this result is the very different behavior between $\varphi =0$ and $\varphi >0$. In the latter case, asymptotically we detect the position of the changepoints without error. This is because the positive autocorrelation in the noise across the changepoint helps us detect it. In fact, as $n\to \infty$ the signal for a change at t comes just from the lag-1 difference, $y_{t+1}-y_t$. The variance of $(y_{t+1}-y_t)$ is $O(1/n)$, and its mean is 0 except at changepoints, where it takes a fixed non-zero value. A simple rule based on detecting a change at t if and only if ${(y_{t+1}-y_t)}^2$ is above some threshold, $c_1(\hspace{0.17em}\text{log}\hspace{0.17em}n)/n$ for some suitably large constant c₁, would consistently detect the changes. For the infill asymptotics, we consider, empirically DeCAFS converges to such an approach as $n\to \infty$.

The theorem also gives insight into the choice of penalty β. It is natural to choose this to be the smallest value that ensures consistency, as larger values will mean loss of power for detecting changes. Assuming the DeCAFS model is correct and we have the true hyper parameters this suggests using $\beta =2\hspace{0.17em}\text{log}\hspace{0.17em}n$, the infimum of the penalties that are valid according to the theorem. This is the value that we use within our simulation study—though slightly inflating the penalty may be beneficial to account for error in the estimated hyper parameters or if we want to account for substantial model error.

6 Simulation Study

6.1 Comparison with Changepoint Methods

We now asses the performance of our algorithm in a simulation study on four different change scenarios, illustrated in Figure 3. In all cases we run DeCAFS with $\beta =2\hspace{0.17em}\text{log}\hspace{0.17em}n$, and estimate the parameters $\varphi ,{\sigma }_{\eta },{\sigma }_{\nu }$ as described in Section 4.

Fig. 3 Four different change scenarios. Top-left, no change present, top-right, change pattern with 19 different changes, bottom-left up changes only, bottom-right, up-down changes of the same magnitude. In this particular example data were generated from an AR model with $\varphi =0.7,{\sigma }_{\nu }=2$.

Simulations were performed over a range of evenly spaced values of $\varphi ,{\sigma }_{\eta },{\sigma }_{\nu }$. There are no current algorithms that directly model local fluctuations in the mean, so we compare with two approaches the assume a constant mean between changes: FPOP (Maidstone et al. 2017) which also assumes IID noise, and AR1Seg (Chakar et al. 2017) that models the noise as an AR(1) process. We compare default implementation of each method, which involves robust estimates of the assumed model parameters. We also compare an implementation of FPOP with an inflated penalty (Bardwell et al. 2019) to account for the autocorrelated noise. To see the impact of possible misestimation of the model parameters, we also implement DeCAFS and AR1Seg using the true parameters when this is possible.

We focus on the accuracy of these methods at detecting the changepoints. We deem a predict change as correct if it is within ±2 observations of a true changepoint. As a measure of accuracy we use the F1 score, which is defined as the harmonic mean of the precision (the proportion of detected changes that are correct) and the recall (the proportion of true changes that are detected). The F1 score ranges from 0 to 1, where 1 corresponds to a perfect segmentation. Separate figures for precision and recall can be found in Section G of the supplementary material. Results reported are based over 100 replications of each simulation experiment, with each simulated data having n = 5000.

In Figure 4A, we report performances of the various algorithms as we vary $\varphi$ for fixed values of ${\sigma }_{\nu }=2$ and ${\sigma }_{\eta }=0$. In Figure 4B, we additionally fix $\varphi =0.85$, but we vary the size of changes. In these cases, there is no random walk component and the model assumed by AR1Seg is correct.

Fig. 4 F1 Scores on the 4 different scenarios. In A a pure AR(1) over a range of values of $\varphi$, for fixed values of ${\sigma }_{\nu }=2,{\sigma }_{\eta }=0$ and a change of magnitude 10. In B a pure AR(1) process with fixed $\varphi =0.85$ and changes in the signal of various magnitudes. In C the full model with $\varphi =0.85$ for a range of values of ${\sigma }_{\eta }$. The gray line represents the cross-section between parameters values in A, B, and C. AR1Seg est. and DeCAFS est. refer to the segmentation of the relative algorithms with estimated parameters. Note, in B the results from DeCAFS and DeCAFS est overlap so only one line is visible. Other algorithms use the true parameter values.

There are a number of conclusions to draw from these results. First, we see that the impact of estimating the parameters on the performance of DeCAFS and AR1Seg is small. Second, we see that using a method which ignores autocorrelation but just inflates the penalty for a change does surprisingly well unless the autocorrelation is large, $\varphi >0.5$, this is inline with results on the robustness of using a square error cost for detecting changes in mean (Lavielle and Moulines 2000). For high values of $\varphi$, DeCAFS is the most accurate algorithm. The one exception are the simulations where there are no changes: the default penalty choice for AR1Seg is such that it rarely introduces a false positive.

In Figure 4C, we explore the effect of local fluctuations in the mean by varying ${\sigma }_{\eta }$. We see a quick drop off in performance for all methods as ${\sigma }_{\eta }$ increases, consistent with the fact that it is harder to detect abrupt changes when the local fluctuations of the mean are greater. Across all experiments, DeCAFS was the most accurate algorithm.

One word of caution when fitting the full DeCAFS model, is that when ${\sigma }_{\eta }$ is large it can be difficult to estimate the parameters, as a model with a very high random walk variance produces data similar to that of a model with constant mean but high autocorrelation. Whilst the impact on detecting changes of any errors when estimating the parameters is small, it can lead to larger errors in the estimate of the signal, μ_t : as different parameter estimates mean that the fluctuations in the data are viewed as either fluctuations in the noise process or in the signal. An example of this is shown in Section F of the supplementary material.

6.2 Robustness to Model Mis-specification

We now investigate the performance of DeCAFS when its model is incorrect. First we follow Chakar et al. (2017) and simulate data with a constant mean between changes but with the noise process being AR(2), i.e. ${ϵ}_t={\varphi }_1{ϵ}_{t-1}+{\varphi }_2{ϵ}_{t-2}+{\nu }_t$. In Figure 5 we report F1 Scores for DeCAFS and AR1Seg as we vary range ${\varphi }_2$. Obviously as $\left|{\varphi }_2\right|$ increases, all algorithms perform worse, but the segmentations returned from DeCAFS are the more reliable as we increase the level of model error.

Fig. 5 F1 score on different scenarios with AR(2) noise as we vary ${\varphi }_2$. Data simulated fixing ${\sigma }_{\nu }=2,{\sigma }_{\eta }=0$ and ${\varphi }_1=0.3$ over a change of size 20.

Second, we consider local fluctuations in the mean that are generated by a sinusoidal process rather than the random walk model, see Figure 6B. In Figure 6A, we compare performance of DeCAFS and AR1Seg as we vary the frequency of the sinusoidal process. Again we see that DeCAFS gives more reliable segmentations in these cases. In the three change scenarios performance decrease as we increase the frequency of the process. In these cases, it becomes significantly harder to detect any changepoints, however, DeCAFS still has higher scores than AR1Seg since it is more robust and returns fewer false positives.

Fig. 6 In A the F1Score on the 4 scenarios for the Sinusoidal Model for fixed amplitude of 15, changes of size 5 and IID Gaussian noise with a variance of 4, as we vary the frequency of the sinusoidal process. In B an example of a realization for the updown scenario, vertical segments refer to estimated changepoint locations of DeCAFS (in light green) and AR1Seg (in blue).

Additional simulation results, showing the robustness of DeCAFS to the mean fluctuations being from an Ornstein-Uhlenbeck process or to the noise being AR(1) within a segment but independent across segments are shown in Sections F and G in the Supplementary Material.

6.3 Comparison to LAVA

The LAVA method of Chernozhukov, Hansen, and Liao (2017) can be applied to model a signal as the sum of a piecewise constant function and a locally fluctuating function: and is thus a natural alternative to DeCAFS. If we let X denote the n × n matrix whose ith column has i – 1 zeroes followed by $n-i+1$ ones then LAVA can estimate the mean ${\mu }_{1:n}$ as $X(f_{1:n}+g_{1:n})$ where the vectors f and g minimize$${(y-X(f_{1:n}+g_{1:n}))}^T(y-X(f_{1:n}+g_{1:n}))+{\lambda }_1\left|\right|f_{1:n}|{|}_1+{\lambda }_2|\left|g_{1:n}\right|{|}_2^2,$$with $\left|\right|·|{|}_1$ and $\left|\right|·|{|}_2$ denoting, respectively, the L₁ and L₂ norms.

The interpretation of this is that ${(f+g)}_i$ is the change in the mean of the data from time i – 1 to time i. The penalties are such that $f_{1:n}$ is sparse, and thus is modeling the abrupt changes in the mean, while $g_{1:n}$ is dense and is accounting for the local fluctuations. It is possible to show that DeCAFS is equivalent to LAVA if we have independent noise and replace the L₁ norm by an L₀ norm.

We implemented LAVA using the lavash package in R. This implementation of LAVA is substantially slower than DeCAFS, and empirically the computational cost appears to increase with the cube of the number of data points. As a result we compare DeCAFS with LAVA on simulated data of length n = 1000 (for which LAVA takes, on average, 132 seconds and DeCAFS 0.02 seconds per dataset).

LAVA tunes λ₁ by cross-validation. We used a plug-in value for λ₂ based on the oracle choice suggested in Chernozhukov, Hansen, and Liao (2017). This choice depends on the variance of the local fluctuations, i.e. the variance of the random walk component in our model, and we implemented LAVA using both the true variance (denoted LAVA), and using the same estimate of the variance as we use for DeCAFS (denoted LAVA_est). The plug-in approach seemed to give better results, and was substantially faster than using cross-validation to tune λ₂. It also makes a comparison with DeCAFS easier, as in situations where neither method estimates any abrupt changes, the two methods then give essentially identical estimates of the mean.

Results from analyzing data simulated under a pure random-walk model (so $\varphi =0$) are shown in Figure 7. While both methods perform similarly at estimating the underlying mean, we see that LAVA is less reliable at estimating the changepoint locations—and often substantially over-estimates the number of changepoints. As summarized in the introduction, this is not unexpected as it is common for methods that use ${\mathcal{l}}_1$ penalties on the size of an abrupt change to overestimate the number of changes (see, e.g., Fearnhead, Maidstone, and Letchford 2018; Jewell et al. 2020).

Fig. 7 On top: comparison of the F1 Score in A1, Precision in A2 and MSE in A3, for DeCAFS (in green) and LAVA (in red) with oracle initial parameters and the relative results with estimated initial parameters (in lighter colours), on the updown scenario for a random walk signal over a range of values of ${\sigma }_{\eta }$. On the bottom the first 250 observations of two realizations of the experiment with, in B1, ${\sigma }_{\eta }$ equal to 0.5 and in B2 ${\sigma }_{\eta }$ equal to 2. Again, the continuous lines over the data points represent the signal estimates of DeCAFS and LAVA; and the vertical lines below show their estimated changepoint locations.

7 Gene Expression in Bacilus subtilis

We now evaluate DeCAFS on estimating the expression of cells in the bacteria Bacilus subtilis. Specifically, we analyze data from Nicolas et al. (2009), which is data from tiling arrays with a resolution of less than 25 base pairs. Each array contains several hundred thousand probes which are ordered according to their position on the bacterial chromosome. For a probe, labeled t say, we get an RNA expression measure, Y_t . Figure 8 shows data from 2000 probes. Code and data used in our analyses, presented below, are available on forgemia: https://forgemia.inra.fr/guillem.rigaill/decafsrna.

Fig. 8 Data on 2000 bp of the plus-strand of the Bacilus subtilis chromosome. Gray dots show the original data. The plain red line represents the estimated signal of DeCAFS with a penalty of $10\hspace{0.17em}\text{log}\hspace{0.17em}(n)$. The dashed black line represents the estimated signal of hmmTiling.

The underlying expression level is believed to undergo two types of transitions, large changes which Nicolas et al. (2009) call shifts and small changes which they call drifts. Thus, it naturally fits our modeling framework of abrupt changes, the shifts, between which there are local fluctuations caused by the drifts. To evaluate the performance of DeCAFS at estimating how the gene expression levels vary across the genome we will compare to the hmmTiling method of Nicolas et al. (2009). This method fits a discrete state hidden Markov model to the data, with the states being the gene expression level, and the dynamics of the hidden Markov model corresponding to either drifts or shifts. As a comparison of computational cost of the two methods, DeCAFS takes about 7 min to analyze data from one of the strands, each of which contains around 192,000 data points. Nicolas et al. (2009) reported a runtime of 5 h and 36 min to analyze both strands.

A comparison of the estimated gene expression level from DeCAFS and from hmmTiling, for a 2000 base pair region of the genome, is shown in Figure 8. We see a close agreement in the estimated level for most of the region, except for a couple of regions where hmmTiling estimates abrupt changes in gene expression level that DeCAFS does not.

To evaluate which of DeCAFS and hmmTiling is more accurate, we follow Nicolas et al. (2009) and see how well the estimated gene expression levels align with bioinformatically predicted promoters and terminators. A promoter roughly corresponds to the start of a gene, and a terminator the end, and we expect gene expression to increase around a promoter and decrease around a terminator.

For promoters, consider all probe locations t from the tiling chip and consider a threshold parameter δ. We can count the number of probe locations with a predicted difference ${\widehat{d}}_t={\widehat{\mu }}_{t+1}-{\widehat{\mu }}_t$ strictly greater than δ. We call this $R(\delta ).$ Among those probes, we can count how many have a promoter nearby (within 22 base pairs). We call this $M(\delta )$. By symmetry we can define an equivalent measure for terminators. A method is better than another if for the same $R(\delta )$ it achieves a larger $M(\delta ).$

We used a data-driven approach to choose the penalty, β for DeCAFS, benefitting from having separate data from the plus and minus strand of the chromosome. For Figure 9 the penalty was learned on the minus strand data and tested on the plus strand data. More specifically we ran DeCAFS on the minus strand for $\beta =\{2\hspace{0.17em}\text{log}\hspace{0.17em}(n),2.5\hspace{0.17em}\text{log}\hspace{0.17em}(n)\dots 30\hspace{0.17em}\text{log}\hspace{0.17em}(n)\}.$ For each β we computed $M(\delta )$ for a fixed $R(\delta )=750$ and took the β maximizing $M(\delta )$: $8\hspace{0.17em}\text{log}\hspace{0.17em}(n)$ for promoters and $5\hspace{0.17em}\text{log}\hspace{0.17em}(n)$ for terminators.

Fig. 9 Benchmark comparisons. The number of promoters (left) and terminators (right) correctly predicted on the plus strand, $M(\delta )$ using a 22 bp distance cutoff, as a function of the number of predicted breakpoints, $R(\delta )$. Plain black lines are the results of hmmTiling (as reported in Figure 4 of Nicolas et al. 2009)). Dotted black lines are the results of hmmTiling when considering all probes rather than only those called transitions. Plain red lines are the results of DeCAFS using $\beta =8\hspace{0.17em}\text{log}\hspace{0.17em}(n)$ for promoters and $5\hspace{0.17em}\text{log}\hspace{0.17em}(n)$ for terminators. These values were learned on the minus strand using a data-driven approach. The thin dark-green leaning line represent $y=x.$

Figure 9 plots $M(\delta )$ against $R(\delta )$ for the plus strand as we vary δ for DeCAFS and two different estimates from hmmTiling. The first, hmmTiling.ori, are the prediction presented in Nicolas et al. (2009). The second, hmmtTiling.all, are those obtained when using all probes rather than only those called transitions by hmmTiling.

In the case of promoters, the prediction of hmmTiling is slightly better than DeCAFS for lower thresholds but noticeably worse for higher thresholds. In the case of terminators, the prediction of DeCAFS are clearly better than those of hmmTiling. Given that DeCAFS was not developed to analyze such data we believe that its relatively good performances for promoters and better performances for terminators is a sign of its versatility.

Unsurprisingly, we get different results and plots if we vary the value for $R(\delta )$ used to select the penalty or if we use the plus strand as a training dataset and the minus strand as a testing dataset. However, Figure 16 in Appendix F.5 shows that $M(\delta )$ is robust to these choices and that the range of penalty for which DeCAFS is close to or better than hmmTiling is quite large, particularly for terminators.

8 Discussion

There are various ways of developing the DeCAFS algorithm, that build on other extensions of the functional pruning version of optimal partitioning. For example, to make the method robust to outliers, we can use robust losses, such as the bi-weight loss, instead of square error loss to measure our fit to the data (Fearnhead and Rigaill 2019). Alternatively, we can incorporate additional constraints on the underlying mean such as monotonicity (Hocking et al. 2020) or that the mean decays geometrically between changes (Jewell and Witten 2018; Jewell et al. 2020). Also it can be extended to allow for the noise to be independent across segments, but AR(1) within a segment. Finally, the algorithm is inherently sequential and thus should be straightforward to adapt to an online analysis of a data stream.

We do not claim that the method we present in Section 4 for estimating the parameters in our model is best. It is likely that more efficient or more robust methods are possible, for example using different robust estimates of the variances of the k-lag difference data (Rousseeuw and Croux 1993); or using iterative procedures where we estimate the changepoints, and then conditional on these changepoints re-estimate the parameters. Using better estimates should lead to further improvement on the statistical performance we observed in Section 6. Our theoretical results suggest that for estimating changes, mis-estimation of the parameters, or errors in our model for the noise or local fluctuations, can be corrected by inflating the penalty for adding a changepoint. As such, in applications we would suggest implementing the method for a range of penalty values, for example using the CROPS algorithm (Haynes, Eckley, and Fearnhead 2017), and then choosing the number of penalties using criteria that consider how the fit to data improves as we add more changes (e.g., Arlot et al. 2016; Fryzlewicz 2018 a; Arlot 2019); or use training data to evaluate performance for different values of β, as we did in Section 7.

Acknowledgments

We thank Pierre Nicolas for providing the output of hmmTiling on the Bacillus subtilis data and his R code allowing us to generate Figure 9, which closely resembles Figure 4 of Nicolas et al. (2009).

Additional information

Funding

This work was supported by EPSRC grant EP/N031938/1, and an ATIGE grant from Genopole. The IPS2 benefits from the support of the LabEx Saclay Plant Sciences-SPS.