¿Cómo puedo resaltar parches ruidosos en una serie de tiempo?

Tengo muchos datos de series temporales: niveles y velocidades del agua frente al tiempo. Es la salida de una simulación de modelo hidráulico. Como parte del proceso de revisión para confirmar que el modelo está funcionando como se esperaba, tengo que trazar cada serie de tiempo para asegurarme de que no haya "bamboleo" en los datos (ver ejemplo de bamboleo menor a continuación). Usar la interfaz de usuario del software de modelado es una forma bastante lenta y laboriosa de verificar estos datos. Por lo tanto, escribí una macro VBA corta para importar varios bits de datos del modelo, incluidos los resultados en Excel y trazarlos todos a la vez. Espero escribir otra macro VBA corta para analizar los datos de series de tiempo y resaltar cualquier sección que sea sospechosa.

Mi único pensamiento hasta ahora es que podría hacer un análisis sobre la pendiente de los datos. Cualquier lugar donde la pendiente cambie rápidamente de positivo a negativo varias veces dentro de una ventana de búsqueda dada podría clasificarse como inestable. ¿Me estoy perdiendo algún truco más simple? Esencialmente, una simulación "estable" debería proporcionar una curva muy suave. Es probable que cualquier cambio repentino sea el resultado de una inestabilidad en los cálculos.

$1-\alpha$$\alpha$

$1201$$\alpha=0.2$$0$$1$

$\alpha$$\alpha$$0.20$$\alpha$$0.20$

Los detalles de la suavidad no importan mucho. En este ejemplo, un loess suavizar (implementado en Rque loesscon span=0.05localizarla) fue utilizado, pero incluso una ventana media tendría bien hecho. Para suavizar los residuos absolutos ejecuté una media de ventana de ancho 17 (aproximadamente 24 minutos) seguida de una mediana de ventana. Estos suavizados en ventanas son relativamente fáciles de implementar en Excel. Una implementación eficiente de VBA (para versiones anteriores de Excel, pero el código fuente debería funcionar incluso en versiones nuevas) está disponible en http://www.quantdec.com/Excel/smoothing.htm .

R Código

#
# 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))