Прежде чем перейти к рассмотрению модели ARIMA, познакомимся сначала с двумя другими моделями: скользящего среднего и моделью авторегрессии.
Модель скользящего среднего MA($q$)
Модель скользящего среднего порядка $q$ или просто MA($q$) предполагает следующую зависимость даннных:
$$y_t\ =\ \mu + \varepsilon_t + \theta_1 \varepsilon_{t-1} + ... + \theta_q \varepsilon_{t-q},$$
где $y_t$ — стационарный ряд со средним $\mu$, а $\varepsilon_t$ — гауссовский белый шум, то есть $\varepsilon_t \sim \mathcal{N}(0, \sigma^2)$ и независимы.
По сути наш ряд $y_t$ выражается через сумму некоторого фиксированного среднего $\mu$, значения белого шума в текущий момент времени $\varepsilon_t$ и не более $q$ предыдущих значений белого шума, домноженных на некоторые коэффициенты, которые являются параметрами модели.
Рассмотрим некоторые свойства модели MA($q$). Как уже было упомянуто выше, ряд $y_t$ будет являтьcя стационарным со средним $\mu$. Найдем также $\mathsf{D}y_t$. Воспользовавшись свойством независимости для $\varepsilon_t$, можем заключить, что
$$\mathsf{D}y_t = \left(1 + \theta_1^2 + \dots + \theta_q^2\right)\sigma^2.$$
Посчитаем автоковариационную функцию для ряда $y_t$, то есть найдем значение $cov(y_t, y_{t + \tau})$. Легко понять, что если $\tau > q$, то $cov(y_t, y_{t+\tau})$ = 0, т.к. $\varepsilon_t$ независимы. Если же $\tau \leq q$, то тогда
$$cov(y_t, y_{t+\tau}) = \left(\theta_{\tau} + \theta_1\theta_{\tau + 1} + \dots + \theta_{q-\tau}\theta_{q}\right)\sigma^2.$$
Записав более компактно, можем получить:
$$ cov(y_t, y_{t+\tau}) = \begin{cases} \sigma^2\sum_{j=0}^{q-\tau}\theta_j\theta_{\tau + j}, & \tau \leq q; \ 0 & \tau > q; \end{cases} $$
где $\theta_0 = 1$. Из посчитанных значений для дисперсии и ковариационной функции, можете попробовать получить выражение и для автокорреляционной функции. Ее особенностью будет как раз равенство нулю на лаге, превосходящим $q$.
Посмотрим на визуализацию:
Модель авторегрессии AR($p$)
Модель авторегрессии для временного ряда можно записать следующим образом:
где $y_t$ — стационарный ряд, а $\varepsilon_t$ — гауссовский белый шум, то есть $\varepsilon_t \sim \mathcal{N}(0, \sigma^2)$ и независимы. Отметим, что, вообще говоря, для стационарности нужны некоторые условия на коэффициенты $\varphi_1, ..., \varphi_p$.
По сути наш ряд $y_t$ выражается через сумму некоторого фиксированного числа $\alpha$, значения белого шума в текущий момент времени $\varepsilon_t$ и не более $p$ предыдущих значений этого же ряда, домноженных на некоторые коэффициенты, которые являются параметрами модели.
Другими словами, модель AR($p$) — это модель $\textit{линейной регрессией}$ для которой
- Таргет: $y_t$ — значение ряда в момент времени $t$
- Признаки: $y_{t-1}, ..., y_{t-p}$ — значения ряда в предыдущие моменты времени
Введем $L$ — оператор сдвига, обладающий следующими свойствами:
- применение $L$ к ряду дает предыдущее значение этого же ряда: $Ly_t = y_{t-1}.$
- применение $L$ к белому шуму дает предыдущее значение шума: $L\varepsilon_t = \varepsilon_{t-1}.$
- применение $L$ к константе — это константа: $Lc = c.$
Оператор $L$ иногда называют также лаговым оператором. Можно рассматривать функции от оператора сдвига, например, кратное применение оператора $L$: $L^2 y_t = L(L y_t) = L(y_{t-1}) = y_{t-2}$ или $L^{-1} y_t= y_{t+1}$. Для записей некоторых моделей временных рядов будет удобно использовать лаговый многочлен:
$$\varphi(L) = \sum_{i=1}^p \varphi_i L^{i} $$
Обратным к оператору $\varphi(L)$ называют оператор $\varphi^{-1}(L)$ такой, что:
$$\varphi(L)\varphi^{-1}(L)y_t = \varphi^{-1}(L)\varphi(L)y_t=y_t$$
Так, например, для $\lvert \varphi \rvert < 1$ можно заключить, что: $$\frac{1}{1 - \varphi L} = \left(1-\varphi L \right)^{-1} = \sum_{i=1}^{\infty}\varphi^{i}L^{i}$$
Рассмотрим модель AR($p$):
$$y_t\ =\ \alpha + \varphi_1 y_{t-1} + ... + \varphi_p y_{t-p} + \varepsilon_t$$
С помощью оператора сдвига ее можно представить в следующем виде:
$$a(L) y_t\ =\ \alpha + \varepsilon_t,$$
где $a(z) = 1 - \varphi_1 z - ... - \varphi_p z^p$ — характеристический полином.
Сформулируем пару важных утверждений:
- Любой стационарный (в широком смысле) процесс представим в виде $MA(\infty)$, то есть в виде модели скользящего среднего с неограниченным количеством слагаемых (конечное или бесконечное число). Этот результат так же известен как теорема Волда о декомпозиции временного ряда.
- Модель $AR(p)$ задает стационарный временной ряд $\Longleftrightarrow$ все комплексные корни $a(z)=0$ лежат вне единичного круга.
Приведем пояснение второго утверждения. В самом деле, пусть $z_1, ..., z_p$ — все его комплексные корни (их ровно $p$ с учетом кратности), тогда справедливо представление:
$$a(z) = (z-z_1)...(z-z_p) = z_1...z_p \left(1 - \frac{z}{z_1}\right) ... \left(1 - \frac{z}{z_p}\right)$$
Тогда при представлении временного ряда в виде
$$y_t\ =\ \frac{\alpha + \varepsilon_t}{a(L)}$$
и дальнейшего его разложения на простые дроби возникнут слагаемые вида
$$\frac{\varepsilon_t}{1 - \frac{L}{z_j}}.$$
Если при этом $z_j$ лежит внутри единичного круга или на его границе, то соответствующий ряд будет расходящимся. На самом деле, случай $z_j=1$ мы в дальнейшем учтем.
В качестве примера рассмотрим подробнее модель $AR(1)$. Зависимость имеет вид $y_t = \alpha + \varphi y_{t - 1} + \varepsilon_t$, где $\varepsilon_t \sim \mathcal{N}(0, \sigma^2)$. Для данного ряда можно выписать следующие свойства:
- Уравнение $1 - \varphi z = 0$, имеет корень $\lambda = 1/\varphi$.
- Тем самым, $AR(1)$ стационарен $\Longleftrightarrow$ $\lvert\varphi\rvert < 1$. Кроме того, чем меньше $\varphi$, тем предыдущее значение ряда вносит меньший вклад в текущее значение.
- Если ряд стационарен, то:
- $\mathsf{E} y_t = \frac{\alpha}{1-\varphi}$
- $\mathsf{D} y_t = \dfrac{\sigma^2}{1-\varphi^2}$
- $cov (y_t, y_{t - h}) = \varphi^h \cdot \dfrac{\sigma^2}{1-\varphi^2}$.
Разберем первое равенство, остальные получаются аналогично. Возьмем математическое ожидание в уравнении ряда
$$\mathsf{E} y_t = \alpha + \varphi \mathsf{E} y_{t - 1} +\mathsf{E} \varepsilon_t$$
Поскольку ряд стационарен, то его математическое ожидание не меняется во времени, а для белого шума математическое ожидание равно нулю. Тем самым мы получаем уравнение на $m = \mathsf{E} y_t$, откуда следует доказываемая формула.
Таким образом, в зависимости от значения $\varphi$ мы можем получить следующие результаты:
- Если $\lvert\varphi\rvert < 1$, то $y_t = \mu + \sum\limits_{j = 0}^{\infty}\varphi^j\varepsilon_{t - j}$ — представление ряда в виде MA($\infty$).
- Если $\lvert\varphi\rvert = 1$, то $AR(1)$ — это случайное блуждание.
- Если $\lvert\varphi\rvert > 1$, то $AR(1)$ — экспоненциально растущий процесс.
Посмотрим на визуализацию.
В первом случае мы имеем модель $y_t = - 0.5 y_{t - 1} + \varepsilon_t$, отрицательный коэффициент является следствием больших колебаний ряда.
Во втором случае модель $y_t = 0.9 y_{t - 1} + \varepsilon_t$, большой положительный коэффициент делает ряд менее шумным.
В третьем случае показано несколько рядов вида случайного блуждания $y_t = y_{t - 1} + \varepsilon_t$, что соответствует случаю $\varphi=1$.
В четвертом случае показан экспоненциальный процесс $y_t = 1.1 y_{t - 1} + \varepsilon_t$, на графике шум уже не заметен из-за масштаба.
На немного вернемся к модели MA($q$). Чуть выше мы выяснили, что при некоторых условиях на коэффиценты $\varphi$ временной ряд модели AR($p$) будет стационарным, а значит имеет представление в виде MA($\infty$). На самом деле, модель скользящего среднего порядка $q$ тоже можно представить с помощью оператора $L$ следующим образом:
$$y_t\ =\ \mu + \varepsilon_t + \theta_1 \varepsilon_{t-1} + ... + \theta_q \varepsilon_{t-q}\ =\ \mu + b(L)\varepsilon_t$$
где $b(z) = 1 + \theta_1 z + ... + \theta_q z^q$ — характеристический многочлен. Для простоты изложения пусть $\mu = 0$. Важным при такой записи оказывается понятие обратимости, то есть представления в виде
$$ \varepsilon_t = b^{-1}(L)y_t, $$
которое означает, что ряд можно представить в виде бесконечной авторегрессионной модели. Здесь, как и в рассуждениях выше, можно заключить, что временной ряд $y_t$ обратим, если все комплексные корни $b(z) = 0$ лежат вне единичного круга.
Модель ARMA($p, q$)
Модель ARMA($p, q$) по сути является суммой моделей $AR(p)$ и $MA(q)$, иначе говоря, модель есть сумма нескольких предыдущих значений ряда и нескольких предыдущих значений белого шума с некоторым коэффициентами.
$$y_t\ =\ \alpha + \varphi_1 y_{t-1} + ... + \varphi_p y_{t-p} + \varepsilon_t + \theta_1 \varepsilon_{t-1} + ... + \theta_q \varepsilon_{t-q}$$
Эквивалентную запись ряда в терминах оператора сдвига можно получить, рассмотрев два многочлена
$$a(L) y_t = \alpha + b(L) \varepsilon_t$$
или
$$y_t = \mu + \frac{b(L)}{a(L)} \varepsilon_t,$$
где $a(z) = 1 - \varphi_1 z - ... - \varphi_p z^p,$ и $b(z) = 1 + \theta_1 z + ... + \theta_q z^q$. Заметим, что во втором представлении константа $\alpha$ заменена на $\mu = \mathsf{E} y_t$. На самом деле, стационарность такого ряда будет определяться только его AR($p$) компонентой, то есть значениями коэффициентов $\varphi$, так ряд в модели MA($q$) всегда является стационарным.
Модель ARIMA($p, d, q$)
Модель ARIMA($p, d, q$) — это расширение моделей типа ARMA на нестационарные временные ряды, которые однако могут стать стационарным после применениея процедуры дифференцирования ряда. Модель ARIMA($p, d, q$) для ряда $y_t$ определяется как модель ARMA($p, q$) для ряда разностей порядка $d$ ряда $y_t$.
- Разность порядка 1: $y_t - y_{t-1} = (1 - L) y_t$.
- Разность порядка 2: $(1 - L)^2 y_t = (1 - L) (y_t - y_{t-1}) = (y_t - y_{t-1}) - (y_{t-1} - y_{t-2}) = y_t - 2y_{t-1} + y_{t-2}$.
Получаем формулу модели ARIMA:
$$a(L) (1 - L)^d y_t = \alpha + b(L) \varepsilon_t$$
или
$$(1 - L)^d y_t = \mu + \frac{b(L)}{a(L)} \varepsilon_t.$$
То есть многочлен $\widetilde{a}(z) = a(z) (1-z)^d$ имеет $d$ единичных корней. Тем самым такая модель позволяет учесть нестационарности, в частности, тренд.
В качестве примера рассмотрим процесс случайного блуждания:
$$y_t = y_{t-1} + \varepsilon_t, $$
где $\varepsilon_t$ — белый шум. Как уже упомяналось ранее, такой ряд не является стационарным. Однако, если мы применим операцию дифференцирования, то можем перейти к новому, уже стационарном ряду $ y't = y_t - y{t-1} $, который можно записать в виде:
$$ y'_t = \varepsilon_t$$
Частичная автокорреляция
Для модели скользящего среднего порядка $q$ мы выяснили, что значения автокорреляционной функции для такого ряда оказывается равной нулю после лага $q$. Эта особенность позволяет использовать автокорреляционную функцию для определения порядка модели скользящего среднего. Возникает разумный вопрос, как оценить порядок $p$ для модели AR($p$)? Здесь оказывается полезным понятие частичной (частной) автокорреляционной функции.
Частичная автокорреляция (PACF) — корреляция ряда с собой после снятия линеной зависимости от промежуточных значений ряда. Иначе говоря, мы хотим как-то учесть опосредованного влияние промежуточных значений ряда и оценить непосредственное влияние $y_{t - \tau}$ на $y_t$. Чуть более формально частичную автокорреляцию можно записать следующим образом:
где $y_t^{\tau-1}$ — линейная регрессия на $y_{t-1}, y_{t-2}, ..., y_{t - (\tau-1)}$:
- $y_t^{\tau-1} = \varphi_1 y_{t-1} + \varphi_2 y_{t-2} + ... + \varphi_{\tau-1} y_{t - (\tau-1)}$
- $y_{t+\tau}^{\tau-1} = \varphi_1 y_{t+\tau-1} + \varphi_2 y_{t+\tau-2} + ... + \varphi_{\tau-1} y_{t+1}$
Пример для $\tau=2$:
$$\gamma_{2} = corr\left(y_{t+2} - \varphi_1 y_{t+1}, y_t - \varphi_1 y_{t-1}\right)$$
где $\varphi_1$ — МНК-оценка в модели $y_t = \varphi y_{t-1}$.
Можно показать, что значение частиной автокорреляции для модели авторегресии AR($p$) будет ненулевой для лагов $\tau \leq p$ и равняться нулю для лагов $\tau > p$. Имеет место быть полная аналогия с автокорреляционной функцией и моделью MA($q$). Таким образом, исследование поведения автокорреляционной и частичной автокорреляционной функции может быть использовано для определения порядка $q$ модели скользящего среднего и порядка $p$ модели авторегрессии соответсвтенно.
Оценка коэффициентов в ARIMA
Пусть гиперпараметры $p, d, q$ фиксированы. В предположении, что $\varepsilon_t$ — гауссовский белый шум, в нашей модели мы можем выписать функцию правдоподобия $L_y(\theta, \varphi, \alpha) = p_{\theta, \varphi, \alpha}(y_1, ..., y_T),$ где $ p_{\theta, \varphi, \alpha}(a_1, ..., a_T)$ — соместная плотность. Из-за того, что $\varepsilon_t$ имеют нормальное распределение, она будет иметь разумный вид. Соответственно, в качестве оценок параметров берется оценка максимального правдоподобия.
Для поиска начальных приближение для параметров $p$ и $q$ воспользуемся автокорреляционной и частичной автокорреляционной функцией.
- Начальное приближение $p$: последний значимый пик у PACF.
- Начальное приближение $q$: последний значимый пик у ACF.
Далее обычно используется поиск по сетке вокруг подобранных значений, минимизируя информационный критерий:
$AIC = -2\ell^* + 2(p+q+1)$ — критерий Акаике;
$AIC_c = -2\ell^* + \frac{2(p+q+1)(p+q+2)}{T-p-q-2}$ — критерий Акаике (короткие ряды);
$BIC = -2\ell^* + (\log T - 2)(p+q+1)$ — Байесовский информационный критерий или критерий Шварца,
где $\ell^* = \ln L_y\left(\widehat{\theta}, \widehat{\varphi}, \widehat{\alpha}\right)$ — логарифм функции правдоподобия, $T$ — длина временного ряда.
Приведем некоторый план при применению модели ARIMA для прогнозирования временных рядов.
Анализ выбросов: замена нерелевантых выбросов на
NA
или усреднение по соседним элементам.Стабилизация дисперсии (преобразования).
Дифференцирование, если ряд не стационарен.
Выбор пилотных $p$ и $q$ по PACF и ACF.
Вокруг этих параметров подбираем оптим. модель по $AIC$/$AIC_c$.
Пошаговое построение прогноза:
— для $t \leqslant T$: $\varepsilon_t \Longrightarrow \widehat{\varepsilon}_t = y_t - \widehat{y}_t$;
— для $t > T$: $\varepsilon_t \Longrightarrow 0$;
— для $t > T$: $y_t \Longrightarrow \widehat{y}_t$.Построение предсказательного интервала:
— если остатки модели нормальны и гомоскедастичны (дисперсия постоянна), то строится теоретический предсказательный интервал
$$\widehat \sigma^2(h) = \widehat \sigma^2 \left(1 + \sum\limits_{i = 1}^{h - 1} \widehat{\psi}_i^2\right)$$
где $h$ — горизонт прогнозирования, $\widehat\sigma^2$ — оценка на дисперсию шума $\varepsilon_t$, $\widehat{\psi}_i$ — коэф. для ряда при его представлении в виде бесконечного процесса скользящего среднего. И $\widehat\sigma^2$, и $\widehat{\psi}_i$ могут быть выражены через оценки на параметры $\varphi$ и $\theta$.
— иначе интервалы строятся с помощью бутстрепа.
Модели SARIMA и ARIMAX
Рассмотрим некоторые расширение модели ARIMA. Обобщение модели ARIMA на ряды с наличием сезонной составляющей назвается SARIMA. Пусть $s$ — известная сезонность ряда. Добавим в модель ARIMA($p, d, q$) компоненты, отвечающие за значения в предыдущие сезоны. Тогда модель SARIMA $(p, d, q)\times (P, D, Q)_s$ может быть записана следующим образом:
$$(1 - L)^d (1 - L^s)^D y_t = \mu + \frac{b(L) B(L^s)}{a(L) A(L^s)} \varepsilon_t,$$
где
$$a(z) = 1 - \varphi_1 z - ... - \varphi_p z^p,$$
$$b(z) = 1 + \theta_1 z + ... + \theta_q z^q,$$
$$A(z) = 1 - \varphi_1^s z - \dots - \varphi_P z^P,$$
$$B(z) = 1 + \theta^s_1 z + \dots + \theta_Q^s z^Q.$$ Параметр сезонного дифференцирования $D$, а также параметры $P, Q$ подбираются из тех же соображений, что и для $p, d, q$, но только с поправкой, что делается это с учетом сезонности $s$. ARIMAX — обобщение модели ARIMA, которая учитывает некоторые экзогенные факторы. Пусть $x_t \in \mathbb{R}^n$ — ряд регрессоров, известный до начала прогноза.
Простой вариант:
$$(1 - L)^d y_t = \mu + \sum_{i=1}^n \frac{\beta_i}{a(L)} x_t^i + \frac{b(L)}{a(L)} \varepsilon_t$$
Общий случай: $$(1 - L)^d y_t = \mu + \sum_{i=1}^n \frac{u_i(L)}{v_i(L)} x_t^i + \frac{b(L)}{a(L)} \varepsilon_t$$
Пример: $x_t = I {\text{в момент времени t праздник}}$
Вышеукзанные модели можно объединить и получить SARIMAX $(p, d, q)\times (P, D, Q)_s $:
$$(1 - L)^d (1 - L^s)^D y_t = \mu + \sum_{i=1}^n \frac{u_i(L)}{v_i(L)} x_t^i + \frac{b(L) B(L^s)}{a(L) A(L^s)} \varepsilon_t$$
Проведем аналогию с линейной регрессией. Это линейная по признакам модель, в которой
- Отклик: $y_t$ — значение ряда в моменты времени $t$
- Признаки:
- $y_{t-1}, ..., y_{t-p}$ — значения ряда в предыдущие моменты времени
- Значение ряда за предыдущие сезоны
- Значения признаков в предыдущие моменты времени
- Значения признаков в предыдущие сезоны
- Ошибка: сумма шума за предыдущие моменты времени и предыдущие сезоны.