Bagaimana saya bisa menyoroti tambalan bising dalam rangkaian waktu?

Saya punya banyak data deret waktu - ketinggian dan kecepatan air vs waktu. Ini adalah output dari simulasi model hidrolik. Sebagai bagian dari proses peninjauan untuk mengonfirmasi bahwa model berkinerja seperti yang diharapkan, saya harus merencanakan setiap rangkaian waktu untuk memastikan bahwa tidak ada "goyangan" dalam data (lihat contoh goyangan kecil di bawah). Menggunakan UI perangkat lunak pemodelan adalah cara yang cukup lambat dan melelahkan untuk memeriksa data ini. Karena itu saya telah menulis makro VBA singkat untuk mengimpor berbagai bit data dari model termasuk hasil ke Excel dan plot semuanya sekaligus. Saya berharap dapat menulis makro VBA pendek lain untuk menganalisis data deret waktu dan menyoroti bagian mana pun yang dicurigai.

Satu-satunya pemikiran saya sejauh ini adalah bahwa saya bisa melakukan beberapa analisis pada kemiringan data. Di mana pun lereng cepat berubah dari positif ke negatif beberapa kali dalam jendela pencarian yang diberikan dapat digolongkan tidak stabil. Apakah saya melewatkan trik yang lebih sederhana? Pada dasarnya, simulasi "stabil" harus memberikan kurva yang sangat halus. Perubahan mendadak apa pun kemungkinan merupakan hasil dari ketidakstabilan dalam perhitungan.

$1-\alpha$$\alpha$

$1201$$\alpha=0.2$$0$$1$

$\alpha$$\alpha$$0.20$$\alpha$$0.20$

Detail dari smooth tidak terlalu penting. Dalam contoh ini loess sebuah halus (dilaksanakan Rsebagai loessdengan span=0.05melokalisasi itu) digunakan, tetapi bahkan berjendela berarti harus dilakukan dengan baik. Untuk menghaluskan residu absolut, saya menjalankan rata-rata berjendela 17 (sekitar 24 menit) diikuti oleh median berjendela. Smoothing berjendela ini relatif mudah diimplementasikan di Excel. Implementasi VBA yang efisien (untuk versi Excel yang lebih lama, tetapi kode sumber harus bekerja bahkan dalam versi baru) tersedia di http://www.quantdec.com/Excel/smoothing.htm .

R Kode

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