# How to fit a robust step function to a time series?

I have a somewhat noisy time series that hovers around different levels.

For example, the following data:

I have the solid line data available, and I would like to obtain an estimate for the dashed line. It should be piecewise constant.

What algorithms are appropriate to try out here?

My ideas so far hover around 0-degree P-splines (but how to find out where to place the knots?) or structural break models. A regression tree is the best idea I currently have, but ideally I would be looking for a method that takes into account the fact that the two levels at y=250 are at equal y-values. If I understand correctly, a regression tree would split these two intervals into two different groups, each with a different mean.

The R code that generated it is this:

set.seed(20181118)
true_fct = stepfun(c(100, 200, 250), c(200, 250, 300, 250))
x = 1:400
y = true_fct(x) + rt(length(x), df=1)
plot(x, y, type="l")
lines(x, true_fct(x), lty=2, lwd=3)

• If your data truly look like the simulated ones, then you could hardly do better than computing a windowed median with a very small window: that would reliably detect all the jumps. Estimate the levels using the medians of the responses within each such interval that is detected. Could you therefore indicate whether the implicit assumptions of the simulation--large jumps, piecewise constant medians, and Student t errors--are precisely the assumptions we should make?
– whuber
Commented Nov 16, 2018 at 14:53
• Thanks for your comment! I have two remarks: (1) How would I get the intervals from the windowed median? (2) The assumptions are piecewise constant medians, and noticeable jumps, but I don't know anything about the error distribution, other than the fact that large outliers can happen. Commented Nov 16, 2018 at 15:33
• Sometimes simple non-parametric methods work when the problem is simple. I would like you to simulate a more challenging/realistic data set where there is an embedded arima structure and perhaps a seasonal pulse or two. Comprehensive approaches to problems like this need to consider and isolate autoregressive structure and anomalies while processing. You might post another question and include the slightly more realistic data set. Commented Nov 16, 2018 at 20:44
• I should also add when the level/step shifts are so large vis-a-vis the error process non-parametric methods may play a useful role and less so as that ratio gets smaller Commented Nov 16, 2018 at 20:53

A simple, robust method to handle such noise is to compute medians.

A rolling median over a short window will detect all but the smallest jumps, while medians of the response within intervals between detected jumps will robustly estimate their levels. (You may replace this latter estimate by any robust estimate that is unaffected by the outliers.)

You should tune this approach with real or simulated data to achieve acceptable error rates. For instance, for the simulation in the question I found it good to use the second and 98th percentiles to set thresholds for detecting the jumps. In other circumstances--such as when many jumps might occur--more central percentiles would work better.

Here is the result showing (a) the three jumps as red dots and (b) the four estimated levels as light blue lines.

The jumps are estimated to occur at indexes 100, 200, 250 (which is exactly where the simulation makes them occur) and the resulting levels are estimated at 199.6, 249.8, 300.0, and 250.2: all within 0.4 of the true underlying values.

This excellent behavior persists with repeated simulations (removing the set.seed command at the beginning).

Here is the R code.

#
# Rolling medians.
#
rollmed <- function(x, k=3) {
n <- length(x)
x.med <- sapply(1:(n-k+10), function(i) median(x[i + 0:(k-1)]))
l <- floor(k/2)
c(rep(NA, l), x.med, rep(NA, k-l))
}
y.med <- rollmed(y, k=5)
#
# Changepoint analysis.
#
dy <- diff(y.med)
fourths <- quantile(dy, c(1,49)/50, na.rm=TRUE)
thresholds <- fourths + diff(fourths)*2.5*c(-1,1)
jumps <- which(dy < thresholds[1] | dy > thresholds[2]) + 1

points(jumps, y.med[jumps], pch=21, bg="Red")
#
# Plotting.
#
limits <- c(1, jumps, length(y)+1)
y.hat <- rep(NA, length(jumps)+1)
for (i in 1:(length(jumps)+1)) {
j0 <- limits[i]
j1 <- limits[i+1]-1
y.hat[i] <- median(y[j0:j1])
lines(x[j0:j1], rep(y.hat[i], j1-j0+1), col="skyblue", lwd=2)
}

• +1, but the "changepoint analysis" part of the code may be not totally clear for some users, so maybe you could comment what is happening there?
– Tim
Commented Nov 16, 2018 at 18:58
• @Tim Thank you for the suggestion. The purpose of the first paragraph is to explain that algorithm. I would like to downplay the details of its implementation because they are unimportant: it should suffice to apply any robust outlier screening method to the residuals.
– whuber
Commented Nov 16, 2018 at 19:19
• You might want to consider zoo::rollmedian or a similar function to simplify your code. Commented Nov 17, 2018 at 1:40
• @usεr11852 Thank you. I am aware of zoo but elected not to use it because I am lazy! It was faster and easier to write rollmed than to review the argument calls to whatever function might already be available. Also, I do like how rollmed clearly illustrates what I am doing rather than hiding the details behind a black box.
– whuber
Commented Nov 17, 2018 at 16:23
• No problem. :) (I was certain you knew of zoo, I was uncertain if you did not use it by choice or by accident. Good answer in any case +1) Commented Nov 17, 2018 at 18:26

If you are still interested in smoothing with L0-penalties I would give a look to the following reference: "Visualization of Genomic Changes by Segmented Smoothing Using an L0 Penalty" - DOI: 10.1371/journal.pone.0038230 (a nice intro to the Whittaker smoother can be found in P. Eilers paper "A perfect smoother" - DOI: 10.1021/ac034173t). Of course, in order to achieve your objective you have to work a bit around the method.

In principle, you need 3 ingredients:

1. The smoother - I would use the Whittaker smoother. Also, I will use matrix augmentation (see Eilers and Marx, 1996 - "Flexible Smoothing with B-splines and Penalties", p.101).
2. Quantile regression - I will use the R package quantreg (rho = 0.5) for laziness :-)
3. L0-penalty - I will follow the mentioned "Visualization of Genomic Changes by Segmented Smoothing Using an L0 Penalty" - DOI: 10.1371/journal.pone.0038230

Of course, you would need also a way to select the optimal amount of smoothing. This is done by my carpenter eyes for this example. You could use the criteria in DOI: 10.1371/journal.pone.0038230 (pg. 5, but I did not try it on your example).

You will find a small code below. I left some comments as guide through it.

# Cross Validated example
rm(list = ls()); graphics.off(); cat("\014")

library(splines)
library(Matrix)
library(quantreg)

# The data
set.seed(20181118)
n = 400
x = 1:n
true_fct = stepfun(c(100, 200, 250), c(200, 250, 300, 250))
y = true_fct(x) + rt(length(x), df = 1)

# Prepare bases - Identity matrix (Whittaker)
# Can be changed for B-splines
B = diag(1, n, n)

# Prepare penalty - lambda parameter fix
nb = ncol(B)
D = diff(diag(1, nb, nb), diff = 1)
lambda = 1e2

# Solve standard Whittaker - for initial values
a = solve(t(B) %*% B + crossprod(D), t(B) %*% y, tol = 1e-50)

# est. loop with L0-Diff penalty as in DOI: 10.1371/journal.pone.0038230
p = 1e-6
nit = 100
beta = 1e-5

for (it in 1:nit) {
ao = a

# Penalty weights
w = (c(D %*% a) ^ 2  + beta ^ 2) ^ ((p - 2)/2)
W = diag(c(w))

# Matrix augmentation
cD = lambda * sqrt(W) %*% D
Bp = rbind(B, cD)
yp =  c(y, 1:nrow(cD)*0)

# Update coefficients - rq.fit from quantreg
a = rq.fit(Bp, yp, tau = 0.5)\$coef

# Check convergence and update
da = max(abs((a - ao)/ao))
cat(it, da, '\n')
if (da < 1e-6) break
}

# Fit
v = B %*% a

# Show results
plot(x, y, pch = 16, cex = 0.5)
lines(x, y, col = 8, lwd = 0.5)
lines(x, v, col = 'blue', lwd = 2)
lines(x, true_fct(x), col = 'red', lty = 2, lwd = 2)
legend("topright", legend = c("True Signal", "Smoothed signal"),
col = c("red", "blue"), lty = c(2, 1))


PS. This is my first answer on Cross Validated. I hope it is useful and clear enough :-)

I would consider using Ruey Tsay's paper Outliers, level shifts, and variance changes in time series Differencing model with AR1 and 21 outliers.

We turned off differencng and the level shifts are specifically called out.

• I wonder whether you overlooked the emphasis on "robust" in the question's title, because any method that identifies 18 spurious parameters (corresponding to outliers introduced in the simulation) in addition to the 3 actual jumps can scarcely be considered robust (or parsimonious, for that matter).
– whuber
Commented Nov 16, 2018 at 19:21
• That is a robust solution. I'm not sure why you are against identifying and adjusting for outliers, but there is a world of research supporting doing it and of course our experiences. Those other variables are outliers. I added a graph that shows the historical data and a cleansed version to contrast the difference. Commented Nov 17, 2018 at 13:40
• Could you be explicit about what your estimate of the step function is?
– whuber
Commented Nov 17, 2018 at 19:43
• There is a flag at period 100(x3), 200(x2), 250(x4) that shows the step. The differencing operator makes it a little tougher to see, but the effect is the same. I added a model without differencing. Commented Nov 17, 2018 at 22:23