Guess The Correlation.
Пусть мы посчитали значение автокорреляции. А как понять, значение 0.1 – это много или мало? На этот вопрос может ответить статистический критерий Льюнга-Бокса (Ljung–Box), который проверяет значимость отклонения $r_{\tau}$ от нуля. Основное правило, которое нужно здесь понять – если значение p-value критерия не превосходит $0.05$ (или другого заранее фиксированного порога значимости), то автокорреляция с лагом $\tau$ значима. Это число вычисляется с помощью стандартных статистических пакетов (например, statsmodels в Питоне).
Рассмотрим временной ряд дорожно-транспортных происшествий за 14 лет с дискретностью 1 месяц, то есть с 12 измерениями в год. На графике мы видим явную сезонность. На нижнем графике изображена коррелограмма – график, визуализирующий автокорреляционную функцию. Точками на графике показаны значения автокорреляционной функции. Ее значение в нуле всегда равно 1, так как это автокорреляция ряда с собой же. Также мы видим, что значение $r_{12}$ является локальным максимумом, что означает высокую положительную корреляцию значения ряда за текущий месяц с аналогичным значением год назад. Иными словами, ряд со сдвигом на год ведёт себя «похожим образом», и это подкрепляет наше наблюдение про наличие годичной сезонности. Наоборот, значение $r_{6}$ минимально в своей окрестност и, что означает высокую отрицательную корреляцию значения ряда со значением полгода назад. Закрашенная область визуализирует границу незначимой автокорреляции, то есть тех значений автокорреляции, для которых не выявлена статистически значимое отличие от 0 (иначе говоря, доверительный интервал пересекает 0). Так мы видим, что последняя значимая сезонная автокорреляция это – это $r_{24}$. Кроме нее также значима сезонная корреляция $r_{12}$ и корреляция за половину сезона $r_6$. Из несезонных автокорреляций значимы оказались $r_{1}$ и $r_{2}$. Отсюда можно сделать вывод, что для построения прогноза значения $y_t$ имеет смысл рассматривать признаки $y_{t-1}, y_{t-2}, y_{t-6}, y_{t-12}, y_{t-24}$.
Рассмотрим временной ряд потребления электричества в Австралии с дискретностью 30 минут. По графику мы можем заметить две разных сезонности – суточную и недельную. Кроме того, при наличии большего количества данных мы смогли бы увидеть еще и годовую сезонность. На данном графике видно повышенное потребление электричества в январе-феврале, когда в Австралии жарко и работает много кондиционеров. По графику автокорреляций мы видим, что наибольшую корреляцию имеют соседние измерения, а также значения сутки назад, двое суток назад, и т.д. Наоборот, значения 12, 36, ... часов назад имеют отрицательную корреляцию.
Временной ряд $y_t$ называется стационарным
Пример. Рассмотрим временной ряд $y_t = \xi_1 \cos t + \xi_2 \sin t$, где $\xi_1, \xi_2$ и независимы и одинаково распределены, причем $\mathsf{P}(\xi_1=1) = \mathsf{P}(\xi_1=-1) = 1/2$.
Можем заметить, что $\mathbb{E} y_1 = 0$ и $cov(y_t, y_s) = \cos t\cos s\ \mathsf{D} \xi_1 + \sin t\sin s\ \mathsf{D} \xi_2 = \cos (t-s)$. Тем самым имеем стационарность в широком смысле.
Но при $t=0$ получаем $y_0 = \xi_1 \in {-1, 1}$, а при $t = \pi/4$ получаем
$$y_{\pi/4} = \frac{\xi_1 + \xi_2}{\sqrt{2}} \in \left\{-\sqrt{2}, 0, \sqrt{2}\right\}.$$
Мы получили разные распределения, поэтому нет стационарности в узком смысле.
Некоторые примеры нестационарных временных рядов:
$$\mathbb{E} y_t=\begin{cases} -1, \text{ при }t = -\pi/2 + 2\pi k; \ 1, \text{ при }t = t = \pi/2 + 2\pi k. \end{cases}$$
Стационарный ряд визуально не имеет предсказуемых закономерностей. Если посмотреть на график такого ряда издалека, то он будет горизонтален.
Ряд можно проверить строго на станционарность с помощью различных статистических критериев. Наиболее популярны следующие критерии:
Рассмотрим несколько примеров временных рядов:
Когда вы определяете, чем вызвано изменение данных – трендом или шумом – стоит учитывать природу данных. Например, колебания значений ряда f теоретически можно было бы объяснить шумом в данных, но по временной оси мы видим, что данные представлены за 15 лет, соответственно, понимаем данные колебания как изменения тренда. Аналогично, для ряда d мы говорим о наличии меняющегося тренда помимо годовой сезонности.
Зачем?
Данные методы рекомендуется использовать, если задача требует некоторой аналитики временного ряда. Если же требуется только построить точечный прогноз на будущее без построения предсказательных интервалов, то стабилизация дисперсии не является необходимой процедурой. Если же нас интересует предсказательный интервал, то многие методы лучше обрабатывают именно стационарные ряды, поэтому имеет смысл стабилизировать дисперсию.
Преобразования:
Класс преобразований Бокса-Кокса с параметром $\lambda$:
$$z_t = \begin{cases} \ln y_t, & \lambda = 0 \ (y_t^\lambda - 1) / \lambda, & \lambda \not=0 \end{cases}.$$
Если есть предположения о зависимости $\mathsf{D} y_t$ от $t$, то можно рассмотреть ряд $z_t = y_t \left/ \sqrt{\mathsf{D} y_t} \right.$.
После построения прогноза для преобразованного ряда нужно сделать обратное преобразование.
Преобразования:
Дифференцирование ряда, то есть переход к ряду $(y'_t, t \in \{2, ..., T\})$, где $y'_t = y_t - y_{t-1}$. Данное преобразование используется для снятие тренда.
Сезонное дифференцирование ряда, то есть переход к ряду $(y'_t, t \in \{s+1, ..., T\})$, где $y'_t = y_t - y_{t-s}$, $s$ – длина сезона.
Преобразования можно применять несколько раз. Обычно сначала применяют сезонное дифференцирование.
Посмотрим на пример. В критерии KPSS для исходного ряда $pvalue < 0.01$, то есть ряд можно считать нестационарным. После логарифмирования ряда $pvalue < 0.01$, а после ещё и сезонного дифференцированная $pvalue > 0.1$, тем самым полученный ряд мы уже не можем отличить от стационарного.
Не редко временной ряд выглядит довольно шумным, что может достаточно плохо сказаться на работе других моделей и подходов к анализу этого временного ряда. В таком случае можно попытаться сгладить значения ряда. Далее мы рассмотрим несколько моделей сглаживания ряда, в том числе при наличии тренда и сезонности ряда. Помимо сглаживания истории ряда, с помощью данных методов можно также осуществлять простое прогнозирование ряда.
Пусть имеется временной ряд $y_t$. В результате экспоненциального сглаживания получается новый временной ряд $\widehat{y}$ по правилу
$$\widehat{y}_{t+1\vert t} = \alpha y_t + (1 - \alpha) \widehat{y}_{t\vert t - 1},$$
где $\widehat{y}_{t+h\vert t}$ – прогноз значения $y_{t+h}$ в момент времени $t$, а $\alpha$ – *параметр сглаживания*.
Смысл преобразования следующий – сглаженное значение в момент времени $t+1$ есть взвешенная комбинация предыдущего значения ряда $y_t$ и предыдущего сглаженного значения ряда $\widehat{y}_{t\vert t - 1}$.
Свойства:
$$\sum\limits_{t=t_0}^T \left(\widehat{y}_{t}(\alpha) - y_t\right)^2 \to \min_{\alpha}.$$
Существуют следующие эмпирические правила:
Примеры: на каждом из графиков изображен исходный ряд (синий) и сглаженный ряд (оранжевый) для разных значений параметра сглаживания. Если имеется тренд или сезонность, то при большом сглаживании полученный ряд начинает «запаздывать» за исходным рядом.
Откуда взялась эта формула экспоненциального сглаживания?
Покажем, что сглаженное значение соответствует прогнозу величины $x$ в момент времени $T+1$, подбираемому по правилу
$$\sum\limits_{t = 0}^T\beta^{T - t}(y_{t} - x)^2 \rightarrow \min\limits_x.$$
Иначе говоря, для прогнозирования мы берем взвешенный MSE с экспоненциально убывающими по времени весами.
Приравняем производную к нулю:
$$2\sum\limits_{t = 0}^T\beta^{T - t}(y_{t} - x) = 0$$
Отсюда выразим $x$ и воспользуемся разложением функции $\frac{1}{1-x}$ в ряд Тейлора
$$x = \dfrac{\sum\limits_{t = 0}^{T}\beta^{T - t}y_t}{\sum\limits_{t = 0}^{T}\beta^t} \approx \dfrac{\sum\limits_{t = 0}^{T}\beta^{T - t}y_t}{1/(1 - \beta)} =(1 - \beta) \sum\limits_{t = 0}^T\beta^{T - t}y_t =$$
$$ = (1 - \beta) y_T + (1 - \beta) \beta \sum\limits_{t=0}^{T - 1}\beta^{T - 1 - t}y_t = (1 - \beta)y_T + \beta \widehat{y}_{T|T - 1}$$
Тем самым мы получили модель экспоненциального сглаживания для $\beta = 1 - \alpha$.
Аддитивный линейный тренд:
Прогноз на $d$ шагов вперед выражается с помощью линейной функции от числа шагов, где коэффициенты меняются по формулам, аналогичным экспоненциальному сглаживанию
$$\widehat{y}_{t+d|t} = a_t + b_t \cdot d,$$
$$a_t = \alpha y_t + (1-\alpha) (a_{t-1} + b_{t-1})$$
$$b_t = \beta (a_t - a_{t-1}) + (1-\beta) b_{t-1}$$
Модель для мультипликативного линейного тренда выглядит аналогично
$$\widehat{y}_{t+d|t} = a_t b_t^d,$$
$$a_t = \alpha y_t + (1-\alpha) (a_{t-1} b_{t-1})$$
$$b_t = \beta \frac{a_t}{a_{t-1}} + (1-\beta) b_{t-1}$$
Модель Хольта: пояснение формул
Поясним на примере аддитивного тренда, почему формулы для $a_t$ и $b_t$ получаются именно такими.
Заметим, что для $d=0$ и $d=1$ и момента времени $t+1$ получаем $\widehat{y}_{t + 1\vert t + 1} = a_{t + 1}$, $\widehat{y}_{t + 1\vert t} = a_{t} + b_{t}$. Хотелось бы, чтобы эти прогнозы примерно совпадали, то есть чтобы имело место $a_{t + 1} - a_{t} \approx b_{t}$.
Рассмотрим ряд разностей $\Delta y_t = a_{t + 1} - a_{t}$ и задачу константного прогноза для него (то есть $\Delta \widehat{y}_{t \vert t - 1} = b$) методом простого экспоненциального сглаживания
$$\sum\limits_{i = 0}^{t}(1 - \beta)^{t - i}(\Delta y_i - b)^2 \rightarrow \min_b.$$
Ее решение мы уже получили ранее:
$$b_t = \beta \Delta y_{t - 1} + (1 - \beta)b_{t - 1} = \beta (a_t - a_{t-1}) + (1-\beta) b_{t-1}.$$
Получилась формулу для $b_t$ в модели аддитивного тренда.
Далее, мы хотим, чтобы имело место $a_{t + 1} \approx b_{t} + a_{t}$. Рассматривая для $y_t$ экспоненциальное сглаживание, в котором в качестве предыдущего значения сглаженного ряда берется $b_t+a_t$, а в качестве нового – $a_{t + 1}$, получаем
$$a_t = \alpha y_t + (1-\alpha) (a_{t-1} + b_{t - 1}).$$
Аддитивная сезонность с трендом:
$$\widehat{y}_{t+d|t} = a_t + db_t + s_{t-m+(d \text{ mod } m)},$$
$$a_t = \alpha (y_t - s_{t-m}) + (1-\alpha) (a_{t-1} + b_{t-1})$$
$$b_t = \beta (a_t - a_{t-1}) + (1-\beta) b_{t-1}$$
$$s_t = \gamma(y_t - a_t) + (1-\gamma)s_{t-m}$$
где $m$ – длина сезона
Мультипликативная сезонность
Без тренда
$$\widehat{y}_{t+d|t} = a_t \cdot s_{t-m+(d \text{ mod } m)},$$
$$a_t = \alpha (y_t / s_{t-m}) + (1-\alpha) a_{t-1}$$
$$s_t = \gamma(y_t / a_t) + (1-\gamma)s_{t-m}$$
С линейным трендом
$$\widehat{y}_{t+d|t} = (a_t + db_t) s_{t-m+(d \text{ mod } m)},$$
$$a_t = \alpha \frac{y_t}{s_{t-m}} + (1-\alpha) (a_{t-1} + b_{t-1})$$
$$b_t = \beta (a_t - a_{t-1}) + (1-\beta) b_{t-1}$$
$$s_t = \gamma\frac{y_t}{a_t} + (1-\gamma)s_{t-m}$$
Разные модели с трендом и сезонностью
В примерах выше мы видели, что при изменении локального тренда ряда экспоненциальное сглаживание запазывает за значениями ряда при использовании сильного сглаживания. Если же использовать слабое сглаживание, то существенного запаздывания не проиходит, но ряд остается шумным. Если для ряда предполагаются значительные структурные изменения, можно использовать модель адаптивного экспоненциального сглаживания, в которой параметр сглаживания может меняться для разных отрезков временного ряда.
Пусть $\widehat{y}_t$ – прогноз значения $y_t$ в момент времени $t-1$ обычным экспоненциальным сглаживанием, а $\widehat{\varepsilon}_t = y_t - \widehat{y}_t$ – ошибка прогноза, сделанного на шаге $t-1$.
Определим следующие значения
$$E_t = \gamma \widehat{\varepsilon}_t + (1-\gamma)E_{t-1}.$$
$$A_t = \gamma \left|\widehat{\varepsilon}_t\right| + (1-\gamma)A_{t-1}.$$
Обычно берут значения $\gamma \in (0.05, 0.1)$. Чтобы экспоненциальное сглаживание быстро приспосабливалось к резким структурным изменениям берут $\alpha_t = \min \left(|K_t|, 1\right)$.
Однако, у данного подхода есть и недостатки, например,