Функции разложения / STL R не могут полностью извлечь годовую сезонность из ежедневных временных рядов

Question

Контекст в первую очередь, вопросы внизу.

У меня есть данные по суточным осадкам за 10 лет, которые показывают годовую сезонность, которую я пытаюсь смоделировать с использованием методов ARMA, а затем сделать прогноз. Данные здесь , создание объекта временного ряда ниже.

Мне известно, что обычные пакеты и функции R борются с ежедневными временными рядами. Например, функция arima () в Forecast не принимает частоты выше 350, а ts () не принимает нецелые значения для частот (оба из них были бы полезны, поскольку среднее число дней в году составляет 365,25).

Очевидно, что функция msts () в Forecast может принимать нецелые значения для аргумента season.perdiods, поэтому я создал свой временной ряд следующим образом:

the_ts <- msts(data$PRCP, start=c(2007, 10), end=c(2017, 9),  seasonal.periods=c(365.25))

Time Series:
Start = c(2007, 10) 
End = c(2017, 9) 
Frequency = 365 
   [1] 0.09 0.75 1.63 0.06 0.36 0.63 0.76 0.43 0.13 0.00 0.00 0.02 0.31 1.80 0.03 0.19 0.25 0.01 0.00 0.52 0.01 0.00 0.00 0.00 0.00 ... etc

-

plot(the_ts)

Эта серия стационарна, поэтому ее не нужно дифференцировать.

Затем я хочу разложить серию, извлечь сезонность и тренд, оставив оставшиеся остатки, которые в случае успеха должны приблизиться к белому шуму.

Ниже приведен график разложения. Я также пытался использовать декомпозировать () и играл с многочисленными настройками параметров.

the_ts_decomp = stl(the_ts, s.window = "periodic")
plot(the_ts_decomp)

По-видимому, в остатках есть сезонность некоторого типа, поскольку форма сезонных трендов в данных параллельна остаткам. Давайте удалим идентифицированный сезонный компонент и проверим предположительно несезонные данные:

the_ts_deseasonal <- seasadj(the_ts_decomp)
plot(the_ts_deseasonal)

Все еще выглядит довольно сезонно для меня. ACF и PACF (не показаны) подтверждают наличие автокорреляции.

Acf(the_ts_deseasonal, lag.max=1000)
Pacf(the_ts_deseasonal, lag.max=1000)

Очевидно, что при прогнозировании можно учитывать сезонную составляющую ряда, передавая преобразование Фурье в модель ARMA как экзогенный ковариат, как описано автором пакета Прогноз здесь и здесь , но мне неясно, как именно это связано с необходимостью декомпозировать и десезонализировать ряды до подгонки модели.

Я могу сделать следующий прогноз, используя этот метод, основываясь на приведенных выше недостаточно разложенных данных. Прогноз не похож на домашний пробег, а остатки подтверждают, что что-то не так:

the_ts.fit <- auto.arima(the_ts, seasonal=FALSE, xreg=fourier(the_ts, K=5))
plot(forecast(the_ts.fit, h=365, xreg=fourier(the_ts, K=5, h=365)))
tsdisplay(residuals(the_ts.fit), lag.max=1000)

Я не очень хорошо понимаю exreg = fourier решение Хиндмана, и, возможно, мне не хватает того, как правильно применить его к моим обстоятельствам.

Вопрос 1 : Разве я не должен быть в состоянии разложить данные до без трендового, бессезонного белого шума, чтобы восстановить их для прогноза? А как насчет использования exreg = решения Фурье?

Вопрос 2 : Почему приведенный выше код не может извлечь сезонный компонент ряда, и как я могу это исправить?

Вопрос 3 : Какой пакет, функцию или технику можно использовать для определения годовой сезонности 365,25?

Answer 1

Я постараюсь ответить на вопросы:

Вопрос 1. Разве мне не нужно иметь возможность разлагать данные до без трендового, бессезонного белого шума, чтобы восстановить их для прогноза? А как насчет использования решения exreg = Фурье?

Решение Фурье учитывает детерминированную сезонность, т. Е. Предполагает, что сезонная картина не меняется со временем. В целом, я бы не сказал, что для серии с трендовой и сезонной корректировкой необходим белый шум, поскольку в этой серии могут оставаться краткосрочные паттерны.

Вопрос 2. Почему вышеприведенный код не может извлечь сезонный компонент ряда и как я могу это исправить?

Когда вы указываете s.window = "period" в функции stl, вы в основном предполагаете, что сезонный паттерн не меняется со временем. Это может быть одним из источников проблемы.

Вопрос 3: Какой пакет, функцию или технику можно использовать для определения годовой сезонности 365,25?

Функция dsa в пакете dsa, для которой вам нужно преобразовать временной ряд в формат xts. Например, как это:

data = rnorm(365.25*10, 100, 1)
data_xts <- xts::xts(data, seq.Date(as.Date("2007-01-01"), by="days", length.out = length(data)))

sa = dsa::dsa(data_xts, fourier_number = 24) # the fourier_number is used to model monthly recurring seasonal patterns in the regARIMA part

data_adjusted <- sa$output[,1]