Jak mogę wyróżnić hałaśliwe łatki w szeregu czasowym?

Mam wiele danych szeregów czasowych - poziomy wody i prędkości w funkcji czasu. Jest to wynik symulacji modelu hydraulicznego. W ramach procesu przeglądu, aby potwierdzić, że model działa zgodnie z oczekiwaniami, muszę wykreślić każdy szereg czasowy, aby upewnić się, że w danych nie ma żadnych „wahnięć” (patrz przykład mniejszego wahania poniżej). Korzystanie z interfejsu użytkownika oprogramowania do modelowania jest dość powolnym i pracochłonnym sposobem sprawdzania tych danych. Dlatego napisałem krótkie makro VBA, aby zaimportować różne bity danych z modelu, w tym wyniki do Excela i wykreślić je wszystkie naraz. Mam nadzieję napisać kolejne krótkie makro VBA, aby przeanalizować dane szeregów czasowych i wyróżnić wszelkie podejrzane sekcje.

Do tej pory myślałem tylko, że mógłbym przeprowadzić analizę nachylenia danych. Wszędzie tam, gdzie nachylenie gwałtownie zmienia się od wielokrotnego dodatniego do ujemnego w danym oknie wyszukiwania, można je zaklasyfikować jako niestabilne. Czy brakuje mi prostszych sztuczek? Zasadniczo „stabilna” symulacja powinna zapewniać bardzo płynną krzywą. Wszelkie nagłe zmiany mogą wynikać z niestabilności obliczeń.

Dla uproszczenia sugerowałbym przeanalizowanie wielkości (wartości bezwzględnych) reszt w stosunku do solidnego wygładzenia danych. W przypadku automatycznego wykrywania rozważ zastąpienie tych rozmiarów wskaźnikiem: 1, gdy przekroczą jakiś wysoki kwantyl, powiedzmy na poziomie , a 0 w przeciwnym razie. Wygładź ten wskaźnik i wyróżnij wygładzone wartości przekraczające .$1-\alpha$$\alpha$

Grafika po lewej stronie przedstawia punktów danych w kolorze niebieskim wraz z solidnym, lokalnym wygładzeniem w kolorze czarnym. Grafika po prawej stronie pokazuje rozmiary resztek tej gładkości. Czarna linia przerywana to ich 80 percentyl (odpowiadający ). Czerwona krzywa jest skonstruowana jak opisano powyżej, ale została przeskalowana (od wartości i ) do średnicy bezwzględnych reszt do wykreślenia.$1201$$\alpha=0.2$$0$$1$

Zmienianie pozwala kontrolować precyzję. W tym przypadku ustawienie poniżej identyfikuje krótką przerwę w hałasie około 22 godzin, podczas gdy ustawienie większe niż powoduje także szybką zmianę w pobliżu 0 godzin.$\alpha$$\alpha$$0.20$$\alpha$$0.20$

Szczegóły gładkości nie mają większego znaczenia. W tym przykładzie less gładka (realizowane w Rjak loessze span=0.05aby go zlokalizować) użyto, ale nawet z okienkiem średnią zrobiłby grzywny. Aby wygładzić resztki bezwzględne, przeprowadziłem okienkową średnią szerokość 17 (około 24 minut), a następnie środkową okienkową. Te wygładzone okna są stosunkowo łatwe do wdrożenia w programie Excel. Wydajna implementacja VBA (dla starszych wersji programu Excel, ale kod źródłowy powinien działać nawet w nowych wersjach) jest dostępna pod adresem http://www.quantdec.com/Excel/smoothing.htm .

R Kod

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