Загружаем нужные пакеты:

library("knitr") # создание отчётов
opts_chunk$set(fig.align = 'center') # выравнивание картинок по центру

library("ggplot2") # графики
library("forecast") # прогнозирование временных рядов: ARMA и ETS модели
library("xts") # операции с временными рядами
library("dplyr") # манипуляции с данными
theme_set(theme_bw()) # чёрно-белая тема оформления графиков

Рассмотрим ARMA(1,1) модель

\[ y_t = 0.7y_{t-1} + \varepsilon_t + 0.5\varepsilon_{t-1} \]

Рассчитаем теоретические значения ACF и PACF

acf <- ARMAacf(ar = 0.7, ma = 0.5, lag.max = 10)
pacf <- ARMAacf(ar = 0.7, ma = 0.5, lag.max = 10, pacf = TRUE)

Покажем ACF и PACF в общей табличке. По какой-то необъяснимой логике при расчёте теоретических ACF R начинает с нулевого лага, а при расчёте PACF - с первого. Именно поэтому мы пишем tail(acf, -1), это избавит нас от ACF для нулевого лага, всегда равной 1.

apacf <- data_frame(lag = 1:10, acf = tail(acf, -1), pacf = pacf)
apacf
## Source: local data frame [10 x 3]
## 
##      lag        acf         pacf
##    (int)      (dbl)        (dbl)
## 1      1 0.83076923  0.830769231
## 2      2 0.58153846 -0.350649351
## 3      3 0.40707692  0.168750000
## 4      4 0.28495385 -0.083591331
## 5      5 0.19946769  0.041698842
## 6      6 0.13962738 -0.020837353
## 7      7 0.09773917  0.010417169
## 8      8 0.06841742 -0.005208396
## 9      9 0.04789219  0.002604175
## 10    10 0.03352454 -0.001302084

Одно дело - теоретические ACF/PACF другое дело - выборочные.

Зафиксируем стартовое зерно генератора случайных чисел. Затем сгенерируем интересующий нас процесс

set.seed(49)
y <- arima.sim(n = 120, list(ar = 0.7, ma = 0.5))

Покажем график ряда и рядом оценённые ACF/PACF

ggtsdisplay(y)

Оценим несколько конкурирующих моделей:

Оцениваем ARMA(1,1) или ARIMA(1, 0, 1)

arma11 <- Arima(y, order = c(1, 0, 1))
summary(arma11)
## Series: y 
## ARIMA(1,0,1) with non-zero mean 
## 
## Coefficients:
##          ar1     ma1  intercept
##       0.8361  0.3882     0.1182
## s.e.  0.0537  0.0839     0.7470
## 
## sigma^2 estimated as 1.034:  log likelihood=-171.73
## AIC=351.46   AICc=351.81   BIC=362.61
## 
## Training set error measures:
##                        ME     RMSE       MAE      MPE     MAPE      MASE
## Training set -0.002169651 1.004134 0.7857084 38.42599 111.8123 0.8956911
##                    ACF1
## Training set 0.04530066

То есть оценённое ARMA(1,1) уравнение имеет вид: \[ \begin{cases} z_t = y_t - 0.12 \\ z_t = 0.84\cdot z_{t-1} + \varepsilon_t + 0.39\varepsilon_{t-1} \end{cases} \]

Оценим AR(2) или ARMA(2, 0) или ARIMA(2, 0, 0)

ar2 <- Arima(y, order = c(2, 0, 0)) 
summary(ar2)
## Series: y 
## ARIMA(2,0,0) with non-zero mean 
## 
## Coefficients:
##          ar1      ar2  intercept
##       1.2416  -0.3830     0.0619
## s.e.  0.0842   0.0847     0.6287
## 
## sigma^2 estimated as 1.02:  log likelihood=-170.91
## AIC=349.81   AICc=350.16   BIC=360.96
## 
## Training set error measures:
##                        ME     RMSE       MAE      MPE     MAPE      MASE
## Training set -0.002325698 0.997134 0.7879019 41.46666 117.4694 0.8981917
##                    ACF1
## Training set 0.01947124

То есть оценённое AR(2) уравнение имеет вид: \[ \begin{cases} z_t = y_t - 0.06 \\ z_t = 1.24\cdot z_{t-1} - 0.38\cdot z_{t-2} + \varepsilon_t \end{cases} \]

Отберём наилучшую модель по штрафному критерию AIC:

\[ AIC = -2 \cdot \ln L + 2 \cdot k, \] где \(\ln L\) - логарифм функции правдоподобия, а \(k\) - число параметров модели.

Чем больше параметров, \(k\), тем сложнее модель, тем выше AIC. Чем ниже функция правдоподобия, \(L\), то есть, чем ниже вероятность получить имеющиеся данные при данной модели, тем выше AIC.

AIC(ar2)
## [1] 349.8132
AIC(arma11)
## [1] 351.4617

По критерию AIC лучше оказалась модель AR(2).

Можно также использовать BIC \[ BIC = -2 \cdot \ln L + \ln n \cdot k. \]

BIC(ar2)
## [1] 360.9632
BIC(arma11)
## [1] 362.6117

По критерию BIC также лучше AR(2) моделиь.

Проверим остатки выбранной модели на отсутствие автокорреляции с помощью теста Льюнга-Бокса:

\[ H_0: \rho_1 = \rho_2 = \ldots = \rho_k \]

\(H_a\): хотя бы одна из корреляций не равна нулю

\[ LB = n (n+1) \sum_{k=1}^h \frac{\hat \rho_k^2}{n - k} \]

Если применять статистику \(LB\) к исходному ряду, то при верной \(H_0\) статистика имеет \(\chi^2\)-распределение с \(h\) степенями свободы. Если к остаткам ARMA(p, q) модели, то количество степеней свободы падает до \(h - (p+q)\).

Мы оценивали AR(2) модель, поэтому степени свободы падают на \(p+q = 2\). Указываем значение аргумента fitdf равное \(p+q=2\):

resid_ar2 <- resid(ar2)
Box.test(resid_ar2, lag = 10, type = "Ljung-Box", fitdf = 2)
## 
##  Box-Ljung test
## 
## data:  resid_ar2
## X-squared = 3.6701, df = 8, p-value = 0.8856

У нас \(H_0\) не отвергается, значит можно считать, что модель корректно описывает структуру корреляции.

Для несезонных рядов Rob Hyndman рекомендует брать lag, \(h = 10\), для сезонных \(h = 2m\), где \(m\) — периодичность сезонности, т.е. \(h=24\) для месячных данных.

Для ARMA(p, q) можно посмотреть визуально, где лежат корни AR и MA части:

autoplot(ar2)

По выбранной наилучшей модели можно строить прогнозы:

forecast(ar2, h = 5)
##     Point Forecast      Lo 80    Hi 80      Lo 95    Hi 95
## 121       4.079150  2.7849924 5.373308  2.0999065 6.058394
## 122       3.451918  1.3887296 5.515107  0.2965434 6.607293
## 123       2.732169  0.1817341 5.282604 -1.1683845 6.632723
## 124       2.078785 -0.7597748 4.917344 -2.2624174 6.419987
## 125       1.543238 -1.4573955 4.543871 -3.0458347 6.132310

Или можно построить график прогнозов с предиктивными интервалами:

future <- forecast(ar2, h = 5)
autoplot(future)

Но всё уже придумано до нас

Уже

Уже

И особо ленивые могут воспользоваться готовой функцией

model <- auto.arima(y)
summary(model)
## Series: y 
## ARIMA(1,0,1) with zero mean     
## 
## Coefficients:
##          ar1     ma1
##       0.8357  0.3886
## s.e.  0.0537  0.0839
## 
## sigma^2 estimated as 1.026:  log likelihood=-171.74
## AIC=349.49   AICc=349.69   BIC=357.85
## 
## Training set error measures:
##                      ME     RMSE       MAE      MPE     MAPE      MASE
## Training set 0.01201994 1.004247 0.7856238 37.28652 111.0717 0.8955947
##                    ACF1
## Training set 0.04543413

Пакет forecast автоматически отбирает \(ARMA(1, 1)\)

Из известных мне косяков: пакет forecast может отказываться подключаться, если окошко для графиков в Rstudio очень маленькое :)

Чтиво: