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

Question

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

Вот пример, в котором десятое наблюдение отсутствует в 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

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

РЕДАКТИРОВАТЬ: немного изменил код, чтобы прояснить мою проблему. Существуют 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)

Answer 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. Возможно, вы захотите снова изменить подгонянные значения на zoo.

Answer 2

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

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

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

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

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

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

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

Answer 3

В этом случае, я думаю, вы хотите исправить сезонность в модели 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) действительно близка, но это просто удача.С более длинными временными рядами вы сможете откалибровать сезонную коррекцию, чтобы дать вам хорошие прогнозы.Было бы полезно иметь предварительные сведения о том, какие механизмы как для сигнала, так и для сезонной коррекции используются для улучшения результатов выборки.