Como destacar patches barulhentos em uma série temporal?

9

Eu tenho muitos dados de séries temporais - níveis e velocidades da água versus tempo. É a saída de uma simulação de modelo hidráulico. Como parte do processo de revisão para confirmar que o modelo está executando como esperado, tenho que plotar cada série temporal para garantir que não haja "oscilações" nos dados (veja o exemplo de oscilação menor abaixo). Usar a interface do usuário do software de modelagem é uma maneira bastante lenta e trabalhosa de verificar esses dados. Portanto, escrevi uma macro VBA curta para importar vários bits de dados do modelo, incluindo resultados para o Excel e plotar todos de uma vez. Espero escrever outra macro VBA curta para analisar os dados da série temporal e destacar as seções suspeitas.

Meu único pensamento até agora é que eu poderia fazer algumas análises na inclinação dos dados. Qualquer local em que a inclinação mude rapidamente de positiva para negativa várias vezes em uma determinada janela de pesquisa pode ser classificado como instável. Estou perdendo alguns truques mais simples? Essencialmente, uma simulação "estável" deve fornecer uma curva muito suave. Qualquer mudança repentina provavelmente resultará de uma instabilidade nos cálculos.

Exemplo de instabilidade menor

davehughes87
fonte
11
Leia o livro EDA de Tukey para obter um conjunto de métodos simples. No início do livro, por exemplo, ele descreve smoothers simples e seu uso para obter resíduos. Uma subseqüente suavização dos resíduos absolutos traçaria a variabilidade local de suas curvas, subindo alto onde ocorreram alterações rápidas, repentinas ou periféricas e permanecendo baixo. Muitos métodos muito mais sofisticados são possíveis, mas talvez isso seja suficiente. Os smoothers de Tukey são relativamente fáceis de codificar no VBA: eu fiz isso .
whuber
@whuber Isso é essencialmente a força do filtro passa-alta deslizante?
Ameba
@amoeba Talvez. Minha compreensão desses filtros é que eles não são totalmente locais e definitivamente não são robustos, enquanto as esmolas de Tukey têm essas duas propriedades importantes. (Atualmente, as pessoas usam Loess ou GAMs para suavizar, o que é bom, mas elas são muito menos simples de implementar.)
whuber

Respostas:

10

Para simplificar, sugiro analisar os tamanhos (valores absolutos) dos resíduos em relação a uma suavização robusta dos dados. Para detecção automatizada, considere substituir esses tamanhos por um indicador: 1 quando excederem algum quantil alto, digamos no nível e 0 caso contrário. Suavize esse indicador e realce quaisquer valores suavizados que excedam .1αα

Figura

O gráfico à esquerda representa pontos de dados em azul, juntamente com um robusto local suave em preto. O gráfico à direita mostra os tamanhos dos resíduos desse liso. A linha pontilhada preta é o percentil 80 (correspondendo a ). A curva vermelha é construída como descrito acima, mas foi escalada (dos valores de e ) para a faixa média dos resíduos absolutos para plotagem.1201α=0.201

Variação permite controle sobre a precisão. Nesse caso, a configuração menor que identifica uma pequena lacuna no ruído em torno de 22 horas, enquanto a configuração maior que também capta a mudança rápida perto de 0 horas.αα0.20α0.20

Os detalhes do liso não importam muito. Neste exemplo, um loess suavizar (implementado em Rcomo loesscom span=0.05a localizá-lo) foi usado, mas até mesmo uma janela médio teria feito bem. Para suavizar os resíduos absolutos, corri uma média com janelas da largura 17 (cerca de 24 minutos), seguida por uma mediana com janelas. Esses suaves janelas são relativamente fáceis de implementar no Excel. Uma implementação eficiente do VBA (para versões mais antigas do Excel, mas o código-fonte deve funcionar mesmo em novas versões) está disponível em 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))
whuber
fonte
11
+1. De alguma forma, você raspou os dados do gráfico do OP?
Ameba
2
@ Ammoeba Isso seria muito problemático, especialmente para os pedaços perfeitos depois de 15 horas. Olhei uma dúzia de pontos na curva, tracei um spline, inseri alguns pontos intermediários para livrar-nos dos picos estranhos que um spline pode produzir e acrescentei um erro correlacionado fortemente heteroscedástico negativamente. Todo o processo levou apenas alguns minutos e resultou em um conjunto de dados qualitativamente como o mostrado na pergunta.
whuber
Gostaria de saber como você conseguiu os dados do meu enredo! Felicidades! Vou tentar.
davehughes87
FWIW, publiquei o código que usei para fazer a ilustração. Mesmo que não seja VBA, talvez ele esclareça os detalhes. (cc @amoeba)
whuber