Загрузка данных

Ряд данных (IM_T_M) содержит месячные показатели импорта товаров в млрд долларов.

dygraph(IM, main = "Импорт товаров",
        xlab = "Год", 
        ylab = "млрд. долл.") %>%
  dyLegend(show = "follow") %>%
  dyRangeSelector()

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

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

sum(is.na(IM))
## [1] 0

Разложим рассмариваемый ряд на составные части:

autoplot(stl(IM, s.window="periodic")) +
  labs(title = "Декомпозиция временного ряда",
       x = "Время")

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

Преобразование в стационарный ряд

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

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

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

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

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

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

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

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

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

lnpml <- log(IM)
autoplot(lnpml) + 
  labs(title = "Логарифм импорта") +
  xlab("Год")

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

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

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

Чтобы привести ряд к стационарному виду, начнем с применения сезонного дифференцирования.

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

autoplot(deseasonal,
     main = "Логарифм импорта без сезонной компоненты", xlab = "Год") 

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

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

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

detrended <- diff(deseasonal, lag = 1, differences = 1)

autoplot(detrended,
     main = "Логарифм импорта без сезонной компоненты", xlab = "Год") 

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 = -5.5273, Lag order = 6, p-value = 0.01
## alternative hypothesis: stationary

На графике все же есть выброс, остносящийся к кризису 1998 года, которые построенной моделью плохо описывается.

На графике виден стационарный ряд, что подкрепляется результатом теста Дики-Фуллера. По итогу сезонного и простого дифференцирование были подобраны два параметра для модели ARIMA: d и D.

D - это порядок сезонного дифференцирования = 1;

d - это порядок обычного дифференцирования = 1. Можно рассматривать и d = 0, так как уже после сезонного дифференцирования ряд был стационарным.

ACF и PACF

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

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

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

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

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

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

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

Выбор модели прогнозирования

Подбор коэффициентов

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

fit1 <- auto.arima(lnpml, d = 0, D = 1, max.p = 3, max.q = 4, max.P = 2,
  max.Q = 1, max.order = 10, seasonal = TRUE, ic = "aic", stepwise = FALSE,
  parallel = FALSE)
fit1
## Series: lnpml 
## ARIMA(3,0,2)(2,1,1)[12] 
## 
## Coefficients:
##           ar1     ar2     ar3     ma1     ma2    sar1     sar2     sma1
##       -0.1310  0.1406  0.9501  1.0234  0.8984  0.1258  -0.0497  -0.7899
## s.e.   0.0475  0.0286  0.0327  0.0781  0.0370  0.0954   0.0845   0.0770
## 
## sigma^2 estimated as 0.005609:  log likelihood=332.7
## AIC=-647.4   AICc=-646.75   BIC=-614.5

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

for (P in c(0, 1, 2)){
  for (Q in c(0, 1)){
    cat("P=",P,"Q=",Q)
    print(arima(lnpml, order = c(3, 0, 2),
      seasonal = list(order = c(P, 1, Q), period = 12)))
  }
}
## P= 0 Q= 0
## Call:
## arima(x = lnpml, order = c(3, 0, 2), seasonal = list(order = c(P, 1, Q), period = 12))
## 
## Coefficients:
##           ar1     ar2     ar3     ma1     ma2
##       -0.2804  0.3758  0.7553  1.1475  0.6458
## s.e.   0.1400  0.1298  0.0818  0.1809  0.1022
## 
## sigma^2 estimated as 0.008013:  log likelihood = 283.21,  aic = -554.42
## P= 0 Q= 1
## Call:
## arima(x = lnpml, order = c(3, 0, 2), seasonal = list(order = c(P, 1, Q), period = 12))
## 
## Coefficients:
##           ar1     ar2     ar3     ma1     ma2     sma1
##       -0.1367  0.1391  0.9547  1.0242  0.8979  -0.7443
## s.e.   0.0396  0.0246  0.0278  0.0709  0.0340   0.0528
## 
## sigma^2 estimated as 0.005534:  log likelihood = 330.8,  aic = -647.59
## P= 1 Q= 0
## Call:
## arima(x = lnpml, order = c(3, 0, 2), seasonal = list(order = c(P, 1, Q), period = 12))
## 
## Coefficients:
##           ar1     ar2     ar3     ma1     ma2     sar1
##       -0.1818  0.1189  0.9512  1.0708  0.9421  -0.3861
## s.e.   0.0252  0.0210  0.0206  0.0465  0.0238   0.0551
## 
## sigma^2 estimated as 0.006553:  log likelihood = 310.07,  aic = -606.15
## P= 1 Q= 1
## Call:
## arima(x = lnpml, order = c(3, 0, 2), seasonal = list(order = c(P, 1, Q), period = 12))
## 
## Coefficients:
##           ar1     ar2     ar3     ma1     ma2    sar1     sma1
##       -0.1767  0.1470  0.9861  1.1566  0.9966  0.1679  -0.8295
## s.e.   0.0094  0.0108  0.0092     NaN     NaN  0.0819   0.0578
## 
## sigma^2 estimated as 0.00535:  log likelihood = 333.33,  aic = -650.65
## P= 2 Q= 0
## Call:
## arima(x = lnpml, order = c(3, 0, 2), seasonal = list(order = c(P, 1, Q), period = 12))
## 
## Coefficients:
##           ar1     ar2     ar3     ma1     ma2     sar1     sar2
##       -0.1725  0.1259  0.9572  1.0641  0.9289  -0.4966  -0.2563
## s.e.   0.0275  0.0201  0.0212  0.0599  0.0293   0.0590   0.0590
## 
## sigma^2 estimated as 0.006122:  log likelihood = 319.03,  aic = -622.06
## P= 2 Q= 1
## Call:
## arima(x = lnpml, order = c(3, 0, 2), seasonal = list(order = c(P, 1, Q), period = 12))
## 
## Coefficients:
##           ar1     ar2     ar3     ma1     ma2    sar1     sar2     sma1
##       -0.1310  0.1406  0.9501  1.0234  0.8984  0.1258  -0.0497  -0.7899
## s.e.   0.0475  0.0286  0.0327  0.0781  0.0370  0.0954   0.0845   0.0770
## 
## sigma^2 estimated as 0.005452:  log likelihood = 332.7,  aic = -647.4

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

fit2 <- arima(lnpml, order = c(3, 0, 2),
      seasonal = list(order = c(1, 1, 1), period = 12))
fit2
## 
## Call:
## arima(x = lnpml, order = c(3, 0, 2), seasonal = list(order = c(1, 1, 1), period = 12))
## 
## Coefficients:
##           ar1     ar2     ar3     ma1     ma2    sar1     sma1
##       -0.1767  0.1470  0.9861  1.1566  0.9966  0.1679  -0.8295
## s.e.   0.0094  0.0108  0.0092     NaN     NaN  0.0819   0.0578
## 
## sigma^2 estimated as 0.00535:  log likelihood = 333.33,  aic = -650.65

Построение моделейна тестовом периоде (несколько вариантов с разными коэффициентами)

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

IM_ob <- window(lnpml, end = c(2017, 6))
fit2_ob <- arima(IM_ob, order = c(3, 0, 2),
      seasonal = list(order = c(1, 1, 1), period = 12))
IM_test <- window(lnpml, start = c(2017, 7))

autoplot(forecast(fit2_ob, h = 12), PI = F) + autolayer(IM_test)

Сравним выбранные на предыдущем этапе значения: ARIMA(3,0,2)(1,1,1); и те значения, которые были предложены автоматически: ARIMA(3,0,2)(2,1,1). На графике изобразим только тестовый период для максимальной наглядности.

fit2_ob_1 <- arima(IM_ob, order = c(3, 0, 2),
      seasonal = list(order = c(1, 1, 1), period = 12))
fit2_ob_2 <- arima(IM_ob, order = c(3, 0, 2),
      seasonal = list(order = c(2, 1, 1), period = 12))

IM_1 <- forecast(fit2_ob_1, h = 12)
IM_2 <- forecast(fit2_ob_2, h = 12)

autoplot(IM_test, series = "test") + 
  autolayer(IM_1, PI = F, series = "ARIMA(3,0,2)*(1,1,1)") +
  autolayer(IM_2, PI = F, series = "ARIMA(3,0,2)*(2,1,1)")

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

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

Суди по визуальному анализу остатков и по сравнению значения ошибок, модель ARIMA(3,0,2)*(2,1,1) незначительно лучше подходит для постоения прогноза.

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

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

rbind(accuracy(fit2_ob_1),
      accuracy(fit2_ob_2)) %>%
  as_tibble() %>%
  round(4) %>%
  mutate(`Метод` = c('ARIMA(3,0,2)*(1,1,1)', 
                     'ARIMA(3,0,2)*(2,1,1)')) %>%
  select(`Метод`, ME, MAPE, MPE) 
## # A tibble: 2 x 4
##   Метод                    ME  MAPE   MPE
##   <chr>                 <dbl> <dbl> <dbl>
## 1 ARIMA(3,0,2)*(1,1,1) 0.0042  2.48 0.130
## 2 ARIMA(3,0,2)*(2,1,1) 0.0038  2.55 0.116

Прогноз на 5 лет

fit3 <- arima(lnpml, order = c(3, 0, 2),
      seasonal = list(order = c(2, 1, 1), period = 12))

lnfcst <- forecast(fit3, h = 60)
autoplot(lnfcst, main = "Прогноз для логарифмированного ряда") +
         labs(x = "Год",
         fill = "Доверительный интервал")

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

fcst <- exp(lnfcst$mean)
autoplot(IM, series = 'Оригинал') +
  autolayer(fcst, series = 'Прогноз') +
  scale_color_manual(values = c("#CC99FF", "#CC0000")) +
     labs(y = 'млрд.$', x = NULL, title = 'Прогноз импорта на 5 лет (60 периодов)', color = NULL)

Источники