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)

 A: 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)
}

A: 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:


*

*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).

*Quantile regression - I will use the R package quantreg (rho = 0.5) for laziness :-)

*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 :-)
A: 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.

