Как я могу выделить шумные участки во временном ряду?

9

У меня есть много данных временных рядов - уровни воды и скорости против времени. Это результат моделирования гидравлической модели. В качестве части процесса проверки, чтобы подтвердить, что модель работает должным образом, я должен построить каждый временной ряд, чтобы убедиться, что в данных нет «колебаний» (см. Пример незначительного колебания ниже). Использование пользовательского интерфейса программного обеспечения для моделирования - довольно медленный и трудоемкий способ проверки этих данных. Поэтому я написал короткий макрос VBA для импорта различных фрагментов данных из модели, включая результаты, в Excel и построения их всех одновременно. Я надеюсь написать еще один короткий макрос VBA для анализа данных временных рядов и выделения любых подозрительных разделов.

До сих пор я думал только о том, что я могу провести некоторый анализ наклона данных. Везде, где наклон быстро меняется от положительного к отрицательному многократно в данном окне поиска, его можно классифицировать как нестабильное. Я пропускаю какие-то более простые трюки? По сути, «стабильное» моделирование должно обеспечивать очень плавную кривую. Любые внезапные изменения могут быть результатом нестабильности расчетов.

Пример незначительной нестабильности

davehughes87
источник
1
Прочитайте книгу Тьюки EDA для набора простых методов. Например, в начале книги он описывает простые сглаживатели и их использование для получения остатков. Последующее сглаживание абсолютных невязок будет отображать локальную изменчивость ваших кривых, повышаясь там, где у вас есть быстрые, внезапные или отдаленные изменения, и оставаясь на низком уровне. Возможны многие гораздо более сложные методы, но, возможно, этого будет достаточно. Сглаживатели Тьюки относительно легко кодировать в VBA: я сделал это .
whuber
@whuber По сути, это сила скользящего фильтра верхних частот?
амеба
@amoeba Возможно. Мое понимание таких фильтров состоит в том, что они не являются полностью локальными и определенно не устойчивы, тогда как сглаживатели Тьюки обладают обоими этими важными свойствами. (В настоящее время люди используют Loess или GAM для сглаживания, что хорошо, но их гораздо проще реализовать.)
whuber

Ответы:

10

1-αα

фигура

1201αзнак равно0.201

αα0,20α0,20

Детали гладкости не имеют большого значения. В этом примере лесс сглаживать (реализовано в Rкачестве loessс span=0.05локализовать ее) был использован, но даже оконным среднее сделало бы штраф. Чтобы сгладить абсолютные остатки, я запустил оконное среднее значение ширины 17 (около 24 минут), а затем оконное медиану. Эти оконные сглаживания относительно легко реализовать в Excel. Эффективная реализация VBA (для более старых версий Excel, но исходный код должен работать даже в новых версиях) доступна по адресу http://www.quantdec.com/Excel/smoothing.htm .


R Код

#
# 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
источник
1
+1. Вы как-то соскребли данные с сюжета ОП?
амеба
2
@ Amoeba Это было бы слишком много проблем, особенно для волнистых битов через 15 часов. Я наметил дюжину точек на кривой, нанес сплайн, вставил несколько промежуточных точек, чтобы избавиться от странных пиков, которые может вызвать сплайн, и добавил сильно отрицательную гетероскедастическую коррелированную ошибку. Весь процесс занял всего несколько минут и привел к качественному набору данных, подобному показанному в вопросе.
whuber
Интересно, откуда у вас данные с моего участка? Ура! Я попробую.
davehughes87
FWIW, я опубликовал код, который использовал для иллюстрации. Даже если это не VBA, возможно, это прояснит детали. (cc @amoeba)
whuber