¿Cómo puedo resaltar las zonas ruidosas de una serie temporal?

Question

Tengo muchos datos de series temporales: niveles de agua y velocidades frente al tiempo. Es el resultado de la simulación de un modelo hidráulico. Como parte del proceso de revisión para confirmar que el modelo funciona como se espera, tengo que trazar cada serie temporal para asegurarme de que no hay "bamboleos" en los datos (véase el ejemplo de bamboleo menor más abajo). Utilizar la interfaz de usuario del software de modelización es una forma bastante lenta y laboriosa de comprobar estos datos. Por lo tanto, he escrito una breve macro VBA para importar varios datos del modelo, incluidos los resultados, a Excel y trazarlos todos a la vez. Espero escribir otra macro VBA corta para analizar los datos de la serie temporal y resaltar cualquier sección que sea sospechosa.

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

Answer 1

Para simplificar, sugeriría analizar los tamaños (valores absolutos) de los residuos en relación con un suavizado robusto de los datos. Para la detección automática, considere la posibilidad de sustituir esos tamaños por un indicador: 1 cuando superen algún cuantil alto, por ejemplo en el nivel $1-\alpha$ y 0 en caso contrario. Suavizar este indicador y resaltar cualquier valor suavizado que exceda $\alpha$ .

El gráfico de la izquierda muestra $1201$ puntos de datos en azul junto con un suave local robusto en negro. El gráfico de la derecha muestra los tamaños de los residuos de ese suavizado. La línea negra de puntos es su percentil 80 (correspondiente a $\alpha=0.2$ ). La curva roja está construida como se ha descrito anteriormente, pero se ha escalado (a partir de valores de $0$ y $1$ ) al rango medio de los residuos absolutos para su trazado.

Variando $\alpha$ permite controlar la precisión. En este caso, la configuración de $\alpha$ menos de $0.20$ identifica una breve brecha en el ruido alrededor de las 22 horas, mientras que el establecimiento $\alpha$ mayor que $0.20$ también recoge el rápido cambio cerca de las 0 horas.

Los detalles del liso no importan mucho. En este ejemplo se ha utilizado un suavizado de loess (implementado en R como loess con span=0.05 para localizarlo) se utilizó, pero incluso una media con ventana habría funcionado bien. Para suavizar los residuos absolutos, utilicé una media con ventana de 17 (unos 24 minutos), seguida de una mediana con ventana. Estos suavizados de ventana son relativamente fáciles de implementar en Excel. Una eficiente implementación de VBA (para versiones antiguas de Excel, pero el código fuente debería funcionar incluso en las nuevas versiones) 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))