Введение

library(tidyverse)
library(lubridate)
library(zoo)
library(forecast)
library(tseries)
library(fpp2)

Модель ARIMA сочетает в себе модель авторегрессии порядка p (AR(p)), являющуюся линейной комбинацией p предыдущих значений ряда; и модель скользящего среднего порядка q (MA(q)), являющуюся линейной комбинацией q последних значений шумовой компоненты (ошибки прогноза). В модели ARIMA(p,d,q) параметр d показывает порядок дифференцирования ряда (здесь дифференцирование - это вычитание с шагом 1, т.е. из текущего значения вычитается предыдущее). В модели ARIMA(p,d,q)(P,D,Q) параметры, обозначающиеся заглавными буквами обозначают то же самое, что и маленькими буквами, но для сезонности.

Здесь стоит ввести понятие стационарности временного ряда. Предположение о стационарности включает в себя следующие условия:

  • математические ожидания значений ряда равны, то есть он не имеет выраженного тренда

\[ E(y_1) = E(y_2) = ... =E(y_n) = const\]

  • дисперсия значений ряда постоянна

\[ Var(y_1) = Var(y_2) = ...=Var(y_n)= \gamma_0\]

  • ковариация (сила линейной связи) между соседними значениями ряда постоянна

\[ Cov(y_1,y_2) = ... = Cov(y_{n-1},y_n) = \gamma_1 \]

  • последнее условие можно записать в общем виде для любого шага между значениями - сила линейной связи между значениями отстоящими друг от друга на один и тот же шаг постоянна, то есть эта связь не зависит от времени, а зависит только от “расстояния”" между двумя выбранными точками

\[ Cov(y_t,y{t-k}) = \gamma_k\]

Существует теорема Вольда, которая гласит, что любой стационарный ряд (распределение значений которого не зависит от времени) может быть описан моделью ARMA(p,q). Соответственно, если ряд имеет тренд либо сезонность (либо и то, и другое), он не является стационарным, а следовательно модель ARMA(p,q) не сможет на нем работать.

Рассмотрим несколько временных рядов: 1) данные о ежедневной аренде велосипедов; 2) данные о поголовье рыси по годам; 3) данные о пассажирообороте на американских авиалиниях по месяцам.

bike <- read.csv('data/day.csv') #чтение файлов данных
passml <- read.csv('data/monthly-us-air-passenger-miles-j.csv')
bike$dteday = as.Date(bike$dteday) #приведение даты к нужному формату
passml$date = as.yearmon(passml$Month)
ggplot(bike, aes(dteday, cnt)) + #построение графиков
  geom_line() + 
  labs(title="Динамика количества арендованных велосипедов по дням") +
  ylab("Количество случаев аренды велосипедов") +
  xlab("День")

autoplot(lynx) + 
  labs(title ="Динамика изменения поголовья рыси по годам") +
  ylab("Количество особей рыси") +
  xlab("Год")

ggplot(passml, aes(date, Passenger_miles)) + 
  geom_line() + 
  labs(title = "Динамика пассажирооборота на авиалиниях США по месяцам") +
  ylab("Пассажиро-мили") +
  xlab("Год")

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

Далее будем рассматривать построение моделей на примере временного ряда о пассажирообороте на авиалиниях США. График выглядит довольно “аккуратно”, выбросов на нем не видно. Если бы были выбросы, можно было бы сгладить выбросы с помощью функции tsclean() из библиотеки forecast.

Проверим ряд на наличие пропущенных значений:

sum(is.na(passml$Passenger_miles))
## [1] 0

Компоненты ряда

Временной ряд может включать в себя несколько компонентов (не всегда включает все). Есть 4 основных составляющих временного ряда - тренд, сезонность, циклическая компонента и шумовая компонента. Разложим рассмариваемый ряд на составные части:

pmiles <- ts(na.omit(passml$Passenger_miles), frequency=12)
decomp <- stl(pmiles, s.window="periodic")
deseasonal_cnt <- seasadj(decomp)
autoplot(decomp) +
  labs(title = "Декомпозиция временного ряда",
       x = "Время")

Как сделать временной ряд стационарным

Проверка на стационарность

Чтобы применить модель ARIMA на нашем временном ряде, необходимо привести его к стационарному виду. Как было видно и до, и после деления ряда на составляющие, наш ряд делают нестационарным три основные вещи: 1) тренд, 2) сезонность, 3) рост дисперсии сезонного размаха. К этому заключению мы пришли после визуального анализа, но нужно что-то более конкретное, чтобы убедиться, стационарен наш ряд или нет.

Для этого случая существует статистический критерий Дики-Фуллера. В качестве нулевой гипотезы обычно рассматривается гипотеза о нестационарности ряда, в качестве альтернативы - то, что ряд стационарен. Соответственно, мы сможем отвергнуть гипотезу о нестационарности временного ряда в том случае, если p-value < 0.05 (при 95% уровне значимости).

Критерий Дики-Фуллера:

Временной ряд: \(y^T = y_1,..., y_T\) ;

Нулевая гипотеза: \(H_0:\) ряд нестационарен;

Альтернативная гипотеза: \(H_1:\) ряд стационарен.

В данном примере будем использовать расширенный тест Дики-Фуллера (augmented Dikey-Fuller test), имеющийся в библиотеке tseries. Однако, у этого теста есть особенность в том, что сезонность ряда как-то влияет на его результат (именно в реализации на R). Поскольку у нас период составляет год, то укажем, что в периоде 12 месяцев.

adf.test(pmiles, k=12, alternative = "stationary")
## 
##  Augmented Dickey-Fuller Test
## 
## data:  pmiles
## Dickey-Fuller = -2.3077, Lag order = 12, p-value = 0.4469
## alternative hypothesis: stationary

Итак, тест показал, что ряд нестационарен.

Сглаживание дисперсии сезонных колебаний

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

Преобразование Бокса-Кокса:

\[ y_t' = \begin{cases} ln y_t & \lambda=0 \\ (y^\lambda_t -1)/\lambda & \lambda \neq0 \end{cases} \]

В этом преобразовании необходимо подбирать параметр \(\lambda\).

Логарифмирование:

\[y_t' = ln y_t \]

passml$lnpml <- log(passml$Passenger_miles)
ggplot(passml, aes(date, lnpml)) + 
  geom_line() + 
  labs(title = "Логарифм пассажирооборота на авиалиниях США") +
  ylab("Логарифм пассажиро-миль") +
  xlab("Месяц")

adf.test(passml$lnpml, k = 12, alternative = "stationary")
## 
##  Augmented Dickey-Fuller Test
## 
## data:  passml$lnpml
## Dickey-Fuller = -0.66073, Lag order = 12, p-value = 0.9727
## alternative hypothesis: stationary

На графике видно, что размах колебаний стал ровнее, чем был, но по критерию Дики-Фуллера ряд все равно не стационарен.

Сезонное дифференцирование

Обычно чтобы привести ряд к стационарному виду применяют дифференцирование (сезонное, если есть сезонность и простое, чтобы избавиться от тренда). Начнем с сезонного дифференцирования, может быть, простое уже не понадобится.

deseasonal <- diff(passml$lnpml, lag = 12, differences = 1)

В данном случае используем функцию plot, так как функция autoplot принимает на вход только объекты класса survfit.

plot(passml$date[c(13:length(passml$lnpml))], deseasonal, type = "l",
     main = "Логарифм пассажирооборота без сезонной компоненты", xlab = "Месяц",
     ylab = "Логарифм пассажиро-миль без сезонных изменений") 

adf.test(deseasonal, alternative = "stationary")
## 
##  Augmented Dickey-Fuller Test
## 
## data:  deseasonal
## Dickey-Fuller = -3.5525, Lag order = 5, p-value = 0.03908
## alternative hypothesis: stationary

После сезонного дифференцирования на графике уже не видно “сезонного паттерна”, а p-value критерия Дики-Фуллера значительно уменьшилось. Но попробуем все же еще провести и простое дифференцирование.

Простое дифференцирование

Чтобы в избавиться от тренда временного ряда проведем дифференцирование с лагом 1.

detrended <- diff(deseasonal, lag = 1, differences = 1)
plot(passml$date[c(14:length(passml$lnpml))], detrended, type = "l",
     main = "Логарифм пассажирооборота без сезонной компоненты и тренда",
     xlab = "Месяц",
     ylab = "Логарифм пассажиро-миль без тренда и сезонности") 

adf.test(detrended, alternative = "stationary")
## Warning in adf.test(detrended, alternative = "stationary"): p-value smaller
## than printed p-value
## 
##  Augmented Dickey-Fuller Test
## 
## data:  detrended
## Dickey-Fuller = -6.9202, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationary

Отлично! Теперь на графике виден стационарный ряд, что подкрепляется результатом теста Дики-Фуллера. Сделав сезонное и простое дифференцирование мы подобрали уже два параметра для модели ARIMA(p,d,q)(P,D,Q) (выделенные жирным). То есть большое D - это порядок сезонного дифференцирования (у нас единица, так как вычитали один раз), а маленькое d - это порядок обычного дифференцирования (тоже единица).

Подбор параметров модели

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

Acf(detrended, lag.max = 48, main='Автокорреляция дифференцированного ряда')

Pacf(detrended, lag.max = 48, main='Частичная автокорреляция дифференцированного ряда')

Начальное приближение параметров подбирают следующим образом:

  • Q - номер последнего сезонного лага (по порядку) со значимой автокорреляцией −> Q = 1

  • q - номер последнего несезонного лага со значимой автокорреляцией −> q = 2

  • P - номер последнего сезонного лага со значимой частичной автокорреляцией −> P = 2

  • p - номер последнего незезонного лага со значимой частичной автокорреляцией −> p = 3

Выбор модели

Обратите внимание, что на вход модели мы будем подавать исходный временной ряд, только логарифмированный (без дифференцирования), порядок дифференцирования укажем как входной параметр.

lnpml <- passml$lnpml

Функция auto.arima() в принципе может сама подбирать параметры при заданных ограничениях, однако, мы укажем найденные начальные приближения параметров в качестве максимально возможных, а также зададим исходные параметры (с которых модель начнет обучаться) единицами.

Качество модели оценивается либо информационным критерием Акаике (AIC), либо Байесовским информационным критерием (BIC). Оба эти критерия очень похожи, можно выбирать любой. Чем меньше значение критерия, тем лучше модель.

Обучим модель:

fit1 <- auto.arima(lnpml, d = 1, D = 1, max.p = 3, max.q = 2, max.P = 2,
  max.Q = 1, max.order = 10, seasonal = TRUE, ic = "aic", stepwise = FALSE,
  parallel = FALSE)
fit1
## Series: lnpml 
## ARIMA(2,1,1) with drift 
## 
## Coefficients:
##          ar1     ar2      ma1   drift
##       0.4431  0.0468  -0.8912  0.0094
## s.e.  0.0763  0.0728   0.0341  0.0016
## 
## sigma^2 estimated as 0.01176:  log likelihood=174.22
## AIC=-338.44   AICc=-338.15   BIC=-321.58

Как видно из результатов, все что касается сезонной компоненты было опущено. Однако, как мы видели на графике, сезонность точно присутствовала в модели. Попробуем самостоятельно методом перебора подобрать сезонные параметры. Порядок дифференцирования уже подобран - 1, для двух остальных пройдемся циклом от нуля до максимального значения, найденного по кореллограммам.

for (P in c(0, 1, 2)){
  for (Q in c(0, 1)){
    cat("P=",P,"Q=",Q)
    print(arima(lnpml, order = c(2, 1, 1),
      seasonal = list(order = c(P, 1, Q), period = 12)))
  }
}
## P= 0 Q= 0
## Call:
## arima(x = lnpml, order = c(2, 1, 1), seasonal = list(order = c(P, 1, Q), period = 12))
## 
## Coefficients:
##          ar1     ar2      ma1
##       0.6608  0.0228  -0.9532
## s.e.  0.0778  0.0750   0.0347
## 
## sigma^2 estimated as 0.005673:  log likelihood = 236.47,  aic = -464.93
## P= 0 Q= 1
## Call:
## arima(x = lnpml, order = c(2, 1, 1), seasonal = list(order = c(P, 1, Q), period = 12))
## 
## Coefficients:
##          ar1      ar2      ma1     sma1
##       0.3853  -0.1148  -0.6581  -0.7427
## s.e.  0.1543   0.0893   0.1442   0.0527
## 
## sigma^2 estimated as 0.003488:  log likelihood = 281.33,  aic = -552.65
## P= 1 Q= 0
## Call:
## arima(x = lnpml, order = c(2, 1, 1), seasonal = list(order = c(P, 1, Q), period = 12))
## 
## Coefficients:
##          ar1      ar2      ma1     sar1
##       0.4744  -0.0572  -0.7684  -0.5231
## s.e.  0.1707   0.1028   0.1606   0.0587
## 
## sigma^2 estimated as 0.004138:  log likelihood = 266.83,  aic = -523.67
## P= 1 Q= 1
## Call:
## arima(x = lnpml, order = c(2, 1, 1), seasonal = list(order = c(P, 1, Q), period = 12))
## 
## Coefficients:
##          ar1      ar2      ma1     sar1     sma1
##       0.3837  -0.1169  -0.6590  -0.1155  -0.6812
## s.e.  0.1473   0.0875   0.1363   0.0968   0.0798
## 
## sigma^2 estimated as 0.003467:  log likelihood = 282.02,  aic = -552.04
## P= 2 Q= 0
## Call:
## arima(x = lnpml, order = c(2, 1, 1), seasonal = list(order = c(P, 1, Q), period = 12))
## 
## Coefficients:
##          ar1      ar2      ma1     sar1     sar2
##       0.4157  -0.1124  -0.7044  -0.6857  -0.2966
## s.e.  0.1393   0.0886   0.1271   0.0671   0.0669
## 
## sigma^2 estimated as 0.003739:  log likelihood = 275.98,  aic = -539.95
## P= 2 Q= 1
## Call:
## arima(x = lnpml, order = c(2, 1, 1), seasonal = list(order = c(P, 1, Q), period = 12))
## 
## Coefficients:
##          ar1      ar2      ma1     sar1    sar2     sma1
##       0.3832  -0.1165  -0.6590  -0.1045  0.0126  -0.6919
## s.e.  0.1476   0.0879   0.1369   0.1366  0.1112   0.1218
## 
## sigma^2 estimated as 0.003467:  log likelihood = 282.02,  aic = -550.05

Ориентируясь на критерий информативности Акаике (он минимизируется), найлучшей моделью оказывается следующая: ARIMA(2,1,1)(0,1,1). Сохраним эту модель.

fit2 <- arima(lnpml, order = c(2, 1, 1),
      seasonal = list(order = c(0, 1, 1), period = 12))
fit2
## 
## Call:
## arima(x = lnpml, order = c(2, 1, 1), seasonal = list(order = c(0, 1, 1), period = 12))
## 
## Coefficients:
##          ar1      ar2      ma1     sma1
##       0.3853  -0.1148  -0.6581  -0.7427
## s.e.  0.1543   0.0893   0.1442   0.0527
## 
## sigma^2 estimated as 0.003488:  log likelihood = 281.33,  aic = -552.65

Анализ остатков модели

Чтобы понять, действительно ли хороша модель, посмотрим на ее остатки.

tsdisplay(residuals(fit2), lag.max = 48, main = "Остатки модели ARIMA (2,1,1)(0,1,1)")

Какими должны быть хорошие остатки модели? Они должны быть независимыми (смотрим, чтобы лаги на коррелограммах не выходили за пределы доверительных интервалов для нулевых значений, т.е. были незначимыми) и в среднем равны нулю (что наблюдается на первом графике).

Также посмотрим и на остатки первой модели, где не учитывался сезонный фактор:

tsdisplay(residuals(fit1), lag.max = 48, main = "Остатки модели ARIMA (2,1,1)")

В данном случае на коррелограммах явно видно, что остатки не независимы.

Прогноз

Посмотрим на двухлетний прогноз прологарифмированного ряда (параметр h функции forecast отвечает за число периодов, на которые надо прогнозировать).

lnfcst <- forecast(fit2, h =24)
autoplot(lnfcst, main = "Прогноз для логарифмированного ряда") +
         labs(x = "Номер периода",
         y = "Логарифм пассажиро-миль",
         fill = "Доверительный интервал")

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

\[ y_f = e^{y_f'}\]

fcst <- exp(lnfcst$mean)

Построим на графике исходный ряд и прогноз.

plot(c(passml$Passenger_miles,fcst), type = "l",
     main = "Прогноз пассажирооборота на 2 года (24 периода)",
     xlab = "Номер периода", ylab = "Пассажиро-мили")
lines(fcst, col = "blue")

Похоже, что прогнозные значения довольно неплохо повторяют тенденции исходного временного ряда.