1. AR(p) – модель, запись в виде с оператором сдвига.

Авторегрессионный процесс порядка p определяется следующим образом \[X_t = c + \sum_{i = 1}^p\phi_i X_{t-i} + \epsilon_t,\] где \(\epsilon_t \sim N(0, \sigma^2)\), независимы между собой и с \(\{x_s\}_{s \leqslant t}\). Пусть существуют такие параметры, что процесс \(x_t\) стационарный.

Это выражение можно переписать через оператор сдвига. Пусть \(LX_t = X_{t-1}\) - оператор сдвига. Тогда, \[\epsilon_t = (1 - \sum_{i = 1}^p\phi_i L^i) = \Phi(L)x_t,\] где \(\Phi\) - характеристический полином.

2. AR(p) и модель сигнала в SSA

В SSA рассматривается модель ряда \(x_n = s_n + \epsilon_n\), где \(s_n = \sum_{i = 1}^p\phi_i s_{t-i}\) - сигнал и \(\epsilon_n\) - шум. Модель авторегрессии имеет вид \(X_t = \sum_{i = 1}^p\phi_i X_{t-i} + \epsilon_t\). В случае SSA шум добавляется ко всему сигналу, а в случае AR - на каждом шаге. Несмотря на то, что модели выглядят похожими, SSA рассматривает авторегрессию только как модель для шума.

3. Вид автоковариационной функции acf для AR(p)

Смоделируем данные для авторегрессии порядка 1 с коэффициентом больше нуля и меньше. Построим для них acf.

set.seed(1)
model.sim <- arima.sim(n = 500, list(order = c(1,0,0), ar = 0.7), sd = sqrt(0.4))
acf(model.sim, main = "phi > 0")

set.seed(1)
model.sim <- arima.sim(n = 500, list(order = c(1,0,0), ar = -0.7), sd = sqrt(0.4))
acf(model.sim, main = "phi < 0")

4. Вид pacf для AR(p)

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

set.seed(1)
model.sim <- arima.sim(n = 500, list(order = c(1,0,0), ar = 0.7), sd = sqrt(0.4))
pacf(model.sim, main = "AR(1)")

set.seed(1)
model.sim <- arima.sim(n = 500, list(order = c(2,0,0), ar = c(0.8, -0.4)), sd = sqrt(0.4))
pacf(model.sim, main = "AR(2)")

5. Модель MA(q), вид acf и pacf

Модель скользящего среднего порядка q определяется следующим образом \[X_t=\epsilon_t + \sum_{j=0}^q\theta_j\epsilon_{t-j},\] где \(\epsilon_t\) - белый шум, \(\theta_t\) - параметры модели.

Представление через оператор сдвига \[X_t=(1+\sum_{j=1}^q \theta_j L^j) \epsilon_t =\Theta(L)\epsilon_t.\]

acf и pacf для MA(q) по сравнению с AR(q) “меняются местами”.

Смоделируем данные для MA(1) с коэффициентом больше нуля и меньше. Построим для них pacf.

set.seed(1)
model.sim <- arima.sim(n = 500, list(order = c(0,0,1), ma = 0.7), sd = sqrt(0.4))
pacf(model.sim, main = "theta > 0")

set.seed(1)
model.sim <- arima.sim(n = 500, list(order = c(0,0,1), ma = -0.7), sd = sqrt(0.4))
pacf(model.sim, main = "theta < 0")

Смоделируем теперь данные для модели MA(1) и MA(2). Построим для них acf. По количеству столбцов, значимо отличающихся от нуля, можно определить количество параметров в модели.

set.seed(1)
model.sim <- arima.sim(n = 500, list(order = c(1,0,0), ar = 0.7), sd = sqrt(0.4))
pacf(model.sim, main = "MA(1)")

set.seed(1)
model.sim <- arima.sim(n = 500, list(order = c(2,0,0), ar = c(0.8, -0.4)), sd = sqrt(0.4))
pacf(model.sim, main = "MA(2)")

6. ARMA(p,q)

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

\[X_t = c + \varepsilon_t + \sum_{i=1}^p \phi_i X_{t-i} + \sum_{i=1}^q \theta_i \varepsilon_{t-i}\]

Запись через оператор сдвига \[\Phi_p(L) X_t = \Theta_q(L)\varepsilon_t\]

acf и pacf для модели ARMA(1,1)

set.seed(1)
model.sim <- arima.sim(n = 500, list(order = c(1,0,1), ar = 0.8, ma = -0.2), sd = sqrt(0.4))
acf(model.sim, main = "ARMA(1,1)")

pacf(model.sim, main = "ARMA(1,1)")

7. ARIMA(p,d,q)

ARIMA — интегрированная модель авторегрессии — скользящего среднего. Дифференцирование используются для того, чтобы нестационарный ряд привести с стационарному виду. Модель ARIMA(p,d,q) означает, что разности временного ряда порядка d подчиняются модели ARMA(p,q).

Запись модели через оператор сдвига \[\Phi_p(L)(1 - L)^d X_t = \Theta_q(L)\varepsilon_t\]

Пример модели ARIMA(2,1,1)

set.seed(1)
model.sim <- arima.sim(n = 500, list(ar = c(0.4, -0.6), ma = 0.3, order = c(2,1,1)), sd = sqrt(0.4))
plot(model.sim)

acf(model.sim, main = "ARIMA(2,1,1)")

pacf(model.sim, main = "ARIMA(2,1,1)")

Продифференцируем ряд

model.sim.diff <- diff(model.sim)
plot(model.sim.diff)

acf(model.sim.diff, main = "ARIMA(2,1,1)")

pacf(model.sim.diff, main = "ARIMA(2,1,1)")

8. Seasonal ARIMA(p,d,q)(P,D,Q)

Модель Seasonal ARIMA имеет тройку дополнительных параметров для сезонности (P, D, Q). Запись модели через оператор сдвига

\[\Phi_P^{(s)}(L^s) \Phi_p(L)X_t= \Theta_Q^{(s)}(L^s)\Theta_q(L)\varepsilon_L,\]

где s - заданный период.

Пример Seasonal ARIMA с предсказанием

library("Rssa")
library("zoo")
setwd("E:/Study/Magistracy/10 term/Time series (pr)/Reports")
data_raw <- read.csv("milk.csv")
data_zoo <- zoo(data_raw[,2], as.Date(data_raw$Month, format = '%Y/%m/%d'))
data <- ts(data_zoo)
plot(data_zoo)

summary(auto.arima(data))
## Series: data 
## ARIMA(1,1,4)                    
## 
## Coefficients:
##           ar1     ma1     ma2      ma3      ma4
##       -0.3045  0.2456  0.1500  -0.4257  -0.6493
## s.e.   0.1158  0.0816  0.0545   0.0486   0.0614
## 
## sigma^2 estimated as 1380:  log likelihood=-839.88
## AIC=1691.77   AICc=1692.29   BIC=1710.48
## 
## Training set error measures:
##                   ME     RMSE      MAE       MPE     MAPE      MASE
## Training set 5.83073 36.48179 28.51598 0.5772274 3.782361 0.7310666
##                     ACF1
## Training set -0.01722656
fit.arima <- arima(data, order = c(1,1,4), seasonal = list(order = c(1,1,4), period = 12))

plot(forecast(fit.arima, h = 24))

9. Exponential smoothing, модели тренда, ES и ARIMA

Пусть \(X_t\) - временной ряд. Экспоненциальное сглаживание определяется выражением \[s_t = \alpha x_t + (1 - \alpha) s_{t-1},\] где \(\alpha\) - параметр сглаживания \((0<\alpha <1),\) \(s_t\) - сглаженный ряд.

Существует эквивалентность между ARIMA и некоторыми моделями экспоненциального сглаживания. Например, модель с постоянным трендом и без сезонности (Simple Exponential smoothing) эквивалентна модели ARIMA(0,1,1). Модель с линейным трендом и без сезонности эквивалентна ARIMA(0,2,2).

Применим экспоненциальное сглаживание к нашим данным.

library(smooth)

dataForExp <- ts(data, deltat = 1/12)

fit.exp <- es(dataForExp, model = "AAA", h=24, holdout = TRUE,interval = TRUE)