Comment mettre en évidence des patchs bruyants dans une série chronologique?

9

J'ai beaucoup de données chronologiques - niveaux d'eau et vitesses en fonction du temps. C'est la sortie d'une simulation de modèle hydraulique. Dans le cadre du processus d'examen pour confirmer que le modèle fonctionne comme prévu, je dois tracer chaque série chronologique pour m'assurer qu'il n'y a pas de "vacillement" dans les données (voir l'exemple de vacillement mineur ci-dessous). L'utilisation de l'interface utilisateur du logiciel de modélisation est un moyen assez lent et laborieux de vérifier ces données. J'ai donc écrit une courte macro VBA pour importer divers bits de données du modèle, y compris les résultats dans Excel et les tracer tous en même temps. J'espère écrire une autre macro VBA courte pour analyser les données de séries chronologiques et mettre en évidence les sections suspectes.

Jusqu'à présent, ma seule pensée est que je pourrais faire une analyse de la pente des données. Partout où la pente passe rapidement de positive à négative plusieurs fois dans une fenêtre de recherche donnée, elle peut être classée comme instable. Suis-je en train de manquer des astuces plus simples? Essentiellement, une simulation "stable" devrait fournir une courbe très lisse. Tout changement soudain est susceptible d'être le résultat d'une instabilité des calculs.

Exemple d'instabilité mineure

davehughes87
la source
1
Lisez le livre EDA de Tukey pour une suite de méthodes simples. Au début du livre, par exemple, il décrit les lissoirs simples et leur utilisation pour obtenir des résidus. Un lissage de suivi des résidus absolus représenterait la variabilité locale de vos courbes, allant haut où vous avez des changements rapides, soudains ou périphériques, et autrement restant bas. De nombreuses méthodes beaucoup plus sophistiquées sont possibles, mais cela suffirait peut-être. Les lisseurs de Tukey sont relativement faciles à coder en VBA: je l'ai fait .
whuber
@whuber C'est essentiellement la puissance du filtre passe-haut coulissant?
amoeba
@amoeba Peut-être. Ma compréhension de ces filtres est qu'ils ne sont pas entièrement locaux et qu'ils ne sont certainement pas robustes, alors que les lisseurs Tukey ont ces deux propriétés importantes. (De nos jours, les gens utilisent Loess ou GAM pour le lissage, ce qui est bien, mais ceux-ci sont beaucoup moins simples à mettre en œuvre.)
whuber

Réponses:

10

1-αα

Figure

1201α=0,201

αα0.20α0,20

Les détails de la douceur importent peu. Dans cet exemple, un loess lisse (implémenté dans Rcomme loessavec span=0.05pour le localiser) a été utilisé, mais même une moyenne fenêtrée aurait fait l'affaire. Pour lisser les résidus absolus, j'ai utilisé une moyenne fenêtrée de largeur 17 (environ 24 minutes) suivie d'une médiane fenêtrée. Ces lissages fenêtrés sont relativement faciles à implémenter dans Excel. Une implémentation VBA efficace (pour les anciennes versions d'Excel, mais le code source devrait fonctionner même dans les nouvelles versions) est disponible sur 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))
whuber
la source
1
+1. Avez-vous en quelque sorte gratté les données du complot du PO?
amoeba
2
@Amoeba Ce serait trop difficile, surtout pour les mèches après 15 heures. J'ai observé une douzaine de points sur la courbe, tracé une spline, inséré des points intermédiaires pour se débarrasser des pointes étranges qu'une spline peut produire et ajouté une erreur corrélée fortement négativement hétéroscédastique. L'ensemble du processus n'a pris que quelques minutes et a abouti à un ensemble de données qualitativement comme celui montré dans la question.
whuber
Je me demandais comment vous aviez obtenu les données de mon intrigue! À votre santé! Je vais essayer.
davehughes87
FWIW, j'ai posté le code que j'ai utilisé pour faire l'illustration. Même si ce n'est pas VBA, cela clarifiera peut-être les détails. (cc @amoeba)
whuber