How can I highlight noisy patches in a time series?

Question

I have lots of time series data - water levels and velocities vs time. It is the output from a hydraulic model simulation. As part of the review process to confirm that the model is performing as expected, I have to plot each time series to ensure that there are no "wobbles" in the data (see example minor wobble below). Using the modelling software's UI is a pretty slow and laborious way to check this data. I've therefore written a short VBA macro to import various bits of data from the model including results into Excel and plot them all at once. I'm hoping to write another short VBA macro to analyse the time series data and highlight any sections that are suspect.

My only thought so far is that I could do some analysis on the slope of the data. Anywhere that the slope rapidly changes from being positive to negative multiple times within a given search window could be classed as unstable. Am I missing any simpler tricks? Essentially, a "stable" simulation should provide a very smooth curve. Any sudden changes are likely to be the result of an instability in the calculations.

Answer 1

For simplicity, I would suggest analyzing the sizes (absolute values) of the residuals relative to a robust smooth of the data. For automated detection, consider replacing those sizes by an indicator: 1 when they exceed some high quantile, say at level $1-\alpha$, and 0 otherwise. Smooth this indicator and highlight any smoothed values that exceed $\alpha$.

The graphic at left plots $1201$ data points in blue along with a robust, local smooth in black. The graphic at right shows the sizes of the residuals of that smooth. The black dotted line is their 80th percentile (corresponding to $\alpha=0.2$). The red curve is constructed as described above, but has been scaled (from values of $0$ and $1$) to the midrange of the absolute residuals for plotting.

Varying $\alpha$ allows control over the precision. In this instance, setting $\alpha$ less than $0.20$ identifies a short gap in the noise around 22 hours, while setting $\alpha$ greater than $0.20$ also picks up the rapid change near 0 hours.

The details of the smooth don't matter much. In this example a loess smooth (implemented in R as loess with span=0.05 to localize it) was used, but even a windowed mean would have done fine. To smooth the absolute residuals I ran a windowed mean of width 17 (about 24 minutes) followed by a windowed median. These windowed smooths are relatively easy to implement in Excel. An efficient VBA implementation (for older versions of Excel, but the source code ought to work even in new versions) is available at http://www.quantdec.com/Excel/smoothing.htm.

---

R Code

#
# Emulate the data in the plot.
#
xy <- matrix(c(0, 96.35,  0.3, 96.6, 0.7, 96.7, 1, 96.73, 1.5, 96.74, 2.5, 96.75, 
               4, 96.9, 5, 97.05, 7, 97.5, 10, 98.5, 12, 99.3, 12.5, 99.35, 
               13, 99.355, 13.5, 99.36, 14.5, 99.365, 15, 99.37, 15.5, 99.375, 
               15.6, 99.4, 15.7, 99.41, 20, 99.5, 25, 99.4, 27, 99.37),
             ncol=2, byrow=TRUE)
n <- 401
set.seed(17)
noise.x <- cumsum(rexp(n, n/max(xy[,1])))
noise.y <- rep(c(-1,1), ceiling(n/2))[1:n]
noise.amp <- runif(n, 0.8, 1.2) * 0.04
noise.amp <- noise.amp * ifelse(noise.x < 16 | noise.x > 24.5, 0.05, 1)
noise.y <- noise.y * noise.amp

g <- approxfun(noise.x, noise.y)
f <- splinefun(xy[,1], xy[,2])
x <- seq(0, max(xy[,1]), length.out=1201)
y <- f(x) + g(x)
#
# Plot the data and a smooth.
#
par(mfrow=c(1,2))
plot(range(xy[,1]), range(xy[,2]), type="n", main="Data", sub="With Smooth",
     xlab="Time (hours)", ylab="Water Level")
abline(h=seq(96, 100, by=0.5), col="#e0e0e0")
abline(v=seq(0, 30, by=5), col="#e0e0e0")
#curve(f(x) + g(x), xlim=range(xy[,1]), col="#2070c0", lwd=2, add=TRUE, n=1201)
lines(x,y, type="l", col="#2070c0", lwd=2)

span <- 0.05
fit <- loess(y ~ x, span=span)
y.hat <- predict(fit)
lines(fit$x, y.hat)
#
# Plot the absolute residuals to the smooth.
#
r <-  abs(resid(fit))
plot(fit$x, r, type="l", col="#808080",
     main="Absolute Residuals", sub="With Smooth and a Threshold",
     xlab="Time hours", ylab="Residual Water Level")
#
# Smooth plot an indicator of the smoothed residuals.
#
library(zoo)
smooth <- function(x, window=17) {
  x.1 <- rollapply(ts(x), window, mean)
  x.2 <- rollapply(x.1, window, median)
  return(as.vector(x.2))
}
alpha <- 0.2
threshold <- quantile(r, 1-alpha)
abline(h=threshold, lwd=2, lty=3)
r.hat <- smooth(r >threshold)
x.hat <- smooth(fit$x)
z <- max(r)/2 * (r.hat > alpha)
lines(x.hat, z, lwd=2, col="#c02020")
par(mfrow=c(1,1))