9
$\begingroup$

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.

Example minor instability

$\endgroup$
3
  • 1
    $\begingroup$ Read Tukey's book EDA for a suite of simple methods. Early in the book, for instance, he describes simple smoothers and their use to obtain residuals. A follow-on smooth of the absolute residuals would chart the local variability of your curves, going high where you have rapid, sudden, or outlying changes, and otherwise staying low. Many much more sophisticated methods are possible, but perhaps this would suffice. Tukey's smoothers are relatively easy to code in VBA: I have done it. $\endgroup$
    – whuber
    Feb 10, 2017 at 14:17
  • $\begingroup$ @whuber This is essentially power of the sliding high-pass filter? $\endgroup$
    – amoeba
    Feb 10, 2017 at 14:51
  • $\begingroup$ @amoeba Maybe. My understanding of such filters is that they are not entirely local and they're definitely not robust, whereas Tukey's smoothers have both of these important properties. (Nowadays people use Loess or GAMs for smoothing, which is fine, but those are much less simple to implement.) $\endgroup$
    – whuber
    Feb 10, 2017 at 14:57

1 Answer 1

11
$\begingroup$

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$.

Figure

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))
$\endgroup$
4
  • 1
    $\begingroup$ +1. Did you somehow scrape the data from the OP's plot? $\endgroup$
    – amoeba
    Feb 10, 2017 at 21:16
  • 3
    $\begingroup$ @Amoeba That would be too much trouble, especially for the wiggly bits after 15 hours. I eyeballed a dozen points on the curve, plotted a spline, inserted some intermediate points to get rid of the strange spikes a spline can produce, and added strongly negatively heteroscedastic correlated error. The whole process took just a few minutes and resulted in a dataset qualitatively like the one shown in the question. $\endgroup$
    – whuber
    Feb 10, 2017 at 21:51
  • $\begingroup$ I wondered how you'd got the data from my plot! Cheers! I'll give it a go. $\endgroup$ Feb 11, 2017 at 10:07
  • $\begingroup$ FWIW, I posted the code I used to make the illustration. Even though it's not VBA, maybe it will clarify the details. (cc @amoeba) $\endgroup$
    – whuber
    Feb 11, 2017 at 16:10

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.