Интерполировать отсутствующие значения в временных рядах с сезонным циклом

У меня есть временной ряд, для которого я хочу разумно интерполировать отсутствующие значения. Значение в определенное время зависит от многодневной тенденции, а также от ее позиции в ежедневном цикле.

Вот пример, в котором десятое наблюдение отсутствует в myzoo

start <- as.POSIXct("2010-01-01") 
freq <- as.difftime(6, units = "hours") 
dayvals <- (1:4)*10 
timevals <- c(3, 1, 2, 4) 
index <- seq(from = start, by = freq, length.out = 16)
obs <- (rep(dayvals, each = 4) + rep(timevals, times = 4))
myzoo <- zoo(obs, index)
myzoo[10] <- NA

Если бы мне пришлось реализовать это, я бы использовал какое-то средневзвешенное значение близких времен в близлежащие дни или добавлял значение для дня к функциональной строке, привязанной к большей тенденции, но я надеюсь, что там уже есть некоторые пакет или функции, которые применяются к этой ситуации?

EDIT: немного изменил код, чтобы прояснить мою проблему. Существуют методы na.*, которые интерполируют из ближайших соседей, но в этом случае они не признают, что недостающее значение находится в момент, когда это самое низкое значение дня. Возможно, решение состоит в том, чтобы преобразовать данные в широкий формат и затем интерполировать, но я не хотел бы полностью игнорировать смежные значения с того же дня. Стоит отметить, что diff(myzoo, lag = 4) возвращает вектор из 10. Решение может быть связано с некоторой комбинацией reshape, na.spline и diff.inv, но я просто не могу понять это.

Вот три подхода, которые не работают:

EDIT2. Изображение, созданное с использованием следующего кода.

myzoo <- zoo(obs, index)
myzoo[10] <- NA # knock out the missing point
plot(myzoo, type="o", pch=16) # plot solid line
points(na.approx(myzoo)[10], col = "red")
points(na.locf(myzoo)[10], col = "blue")
points(na.spline(myzoo)[10], col = "green")
myzoo[10] <- 31 # replace the missing point
lines(myzoo, type = "o", lty=3, pch=16) # dashed line over the gap
legend(x = "topleft", 
       legend = c("na.spline", "na.locf", "na.approx"), 
       col=c("green","blue","red"), pch = 1)

Ответы

Ответ 1

Попробуйте следующее:

x <- ts(myzoo,f=4)
fit <- ts(rowSums(tsSmooth(StructTS(x))[,-2]))
tsp(fit) <- tsp(x)
plot(x)
lines(fit,col=2)

Идея состоит в том, чтобы использовать базовую структурную модель для временного ряда, которая обрабатывает недостающее значение с использованием фильтра Калмана. Тогда гладкая Калмана используется для оценки каждой точки в временном ряду, включая любые пропущенные.

Мне пришлось преобразовать объект вашего зоопарка в объект ts с частотой 4, чтобы использовать StructTS. Возможно, вы захотите снова изменить установленные значения обратно в зоопарк.

Ответ 2

В этом случае, я думаю, вы хотите корректировать сезонность в модели ARIMA. Там не хватает даты здесь, чтобы соответствовать сезонной модели, но это должно помочь вам начать.

library(zoo)
start <- as.POSIXct("2010-01-01") 
freq <- as.difftime(6, units = "hours") 
dayvals <- (1:4)*10 
timevals <- c(3, 1, 2, 4) 
index <- seq(from = start, by = freq, length.out = 16)
obs <- (rep(dayvals, each = 4) + rep(timevals, times = 4))
myzoo <- myzoo.orig <- zoo(obs, index)
myzoo[10] <- NA

myzoo.fixed <- na.locf(myzoo)

myarima.resid <- arima(myzoo.fixed, order = c(3, 0, 3), seasonal = list(order = c(0, 0, 0), period = 4))$residuals
myzoo.reallyfixed <- myzoo.fixed
myzoo.reallyfixed[10] <- myzoo.fixed[10] + myarima.resid[10]

plot(myzoo.reallyfixed)
points(myzoo.orig)

В моих тестах ARMA (3, 3) действительно близко, но это просто удача. С более длинными временными рядами вы сможете откалибровать сезонные поправки, чтобы дать вам хорошие прогнозы. Было бы полезно иметь хорошую предварительную информацию о том, какие базовые механизмы как для сигнала, так и для сезонной коррекции будут лучше соответствовать показателям производительности.

Ответ 3

forecast::na.interp - хороший подход. Из документации

Использует линейную интерполяцию для несезонных серий и периодическое stl-декомпозицию с сезонными рядами для замены отсутствующих значений.

library(forecast)
fit <- na.interp(myzoo)
fit[10]  # 32.5, vs. 31.0 actual and 32.0 from Rob Hyndman answer

В этой статье оценивается несколько методов интерполяции по сравнению с рядами реального времени и находит, что na.interp является точной и эффективной:

Из реализаций R, протестированных в этой статье, na.interp из пакета прогноза и na.StructTS из пакета zoo показал наилучшие общие результаты.

Функция na.interp также не намного медленнее, чем na.approx [самый быстрый метод], поэтому разложение лёсса кажется не очень требовательным с точки зрения вычислительного времени.

Также стоит отметить, что Роб Хиндман написал пакет forecast и включил na.interp после предоставления ответа на этот вопрос. Вероятно, что na.interp является улучшением этого подхода, хотя в этом случае он оказался хуже (вероятно, из-за указания периода в StructTS, где na.interp показывает его).