В этом видео мы разберём, как строить прогноз временного ряда с помощью стандартных и не очень стандартных функций, которые есть в Python'е. Задачи, которые мы будем решать, следующие: известны ежемесячные продажи австралийского вина в тысячах литров с января 1980 года по июль 1995. Мы хотим построить прогноз этого ряда на следующие три года. Давайте начнём. Импортируем все нужные модули и создадим ещё некоторые новые функции, например, функцию, которая делает обратное преобразование Бокса-Кокса. К сожалению, в Python'е нет функции, которая это делает в стандартных пакетах, поэтому мы напишем её сами, но в этом нет ничего сложного. Вот так выглядит ряд, который мы собираемся прогнозировать. Что можно заметить про этот ряд? У него достаточно сильно выраженная годовая сезонность, кроме того кажется, что здесь есть какой-то тренд, не линейный, но тем не менее, возможно, этот тренд можно как-то описать. Кроме того, у этого ряда, кажется, нестационарная дисперсия. Размах сезонных колебаний здесь в начале ряда немного меньше, чем вот здесь, ближе к концу ряда. Прежде чем начинать этот ряд моделировать, давайте ещё немного его предварительно поанализируем. Ну, во-первых, проверим критерием Дики-Фуллера формально стационарность. Мы видим, что достигаемый уровень значимости критерия Дики-Фуллера чуть-чуть больше 0,05, то есть гипотеза нестационарности для этого ряда не отвергается, однако значения — пограничны. Давайте сделаем ещё одно очень полезное действие. Сделаем STL-декомпозицию. STL-декомпозиция — это некоторый набор эвристик, который позволяет вам визуально посмотреть, из каких примерно компонент состоит ваш ряд. Вот результат применения STL-декомпозиции — перед вами, на графике. На верхнем из четырёх графиков — исходный ряд. На нижних трёх соответственно тренд, сезонность и остатки. То есть всё, что не определяется суммой вот этого тренда и сезонности, каким-то образом оцененных. Значит, что мы здесь видим? Сезонный профиль достаточно хорошо выражен. Тренд имеет достаточно сложную структуру. Сначала он повышается, затем доходит до некоторого пика, немного понижается и медленно возвращается на исходный уровень. Это ещё раз нам подтверждает, что ряд трудно считать стационарным. В качестве следующего шага давайте сделаем стабилизацию дисперсии методом Бокса-Кокса. Это делает функция boxcox, она возвращает преобразованный ряд и автоматически каким-то образом наилучшее подобранное значение параметра лямбда ( λ). В данном случае это значение равно примерно 0,24. А вот так выглядит преобразованный ряд. Критерий Дики-Фуллера на этом ряде даёт достигаемый уровень значимости примерно 0,03. То есть нулевая гипотеза нестационарности отвергается, и этот ряд критерий Дики-Фуллера считает стационарным. Тем не менее, мы очень хорошо видим, что стационарным этот ряд считать нельзя, потому что, во-первых, в нём очень сильно выражена сезонность, во-вторых, в нём есть вот тот самый тренд, который мы наблюдали на STL-декомпозиции. Критерий Дики-Фуллера — это всего лишь инструмент, и не очень совершенный. Поэтому, если вы видите, что ваш ряд явно не стационарен, не нужно доверять тому, что говорит вам критерий Дики-Фуллера. Вот и мы не будем. Давайте проведём дифференцирование ряда, для того чтобы добиться стационарности. Начнём с сезонного дифференцирования, сделаем это сезонное дифференцирование, снова применим к продифференцированному ряду критерий Дики-Фуллера и посмотрим на его STL-декомпозицию. На продифференцированном ряде гипотеза нестационарности не отвергается. Достигаемый уровень значимости — примерно 0,013. А вот так выглядит результат STL-декомпозиции. Посмотрите, как здесь изменился тренд. Он практически константный на первом участке ряда, но вот в месте, где происходил в исходном ряде перегиб тренда, мы видим существенный провал. То есть не удивительно, что этот ряд остался нестационарным. Давайте применим ещё одно дифференцирование. На сей раз — обычное. На продифференцированном ещё раз ряде критерий Дики-Фуллера гипотезу нестационарности уверенно отвергает с очень маленьким достигаемым уровнем значимости. А вот так выглядит STL-декомпозиция два раза продифференцированного ряда. На этот раз мы на трендовой компоненте не видим никакого систематического поведения. То есть она колеблется примерно вокруг константы. Это то, чего мы хотели добиться. Видимо, этот ряд уже можно считать стационарным. Перейдём теперь к подбору моделей ARIMA. Давайте посмотрим на автокорреляционную и частично автокорреляционную функцию вот этого два раза продифференцированного ряда. Они выглядят вот так. На верхнем графике — автокорреляционная функция, и по этому графику мы будем подбирать начальное приближение для параметров q и Q. Максимальный сезонный лаг, значимо отличающийся от ноля, то есть вылезающий из вот этого синего коридора, это лаг 12. Таким образом, начальным приближением для Q мы будем брать единицу. Что касается q — здесь значимо отличаются от нуля первые два несезонных лага. Таким образом, для q начальное приближение мы будем брать равным двойке. Перейдём к начальному приближению для P. На частично автокорреляционной функции единственный сезонный лаг, значимо отличающийся от нуля — это лаг 12. Значит P мы возьмём равным единице, для начала. Что касается p — в начале при маленьких лагах на частично автокорреляционной функции есть несколько значимо отличающихся от ноля. Возьмём, допустим, p равным четырём. Прекрасно. Давайте теперь с этим начальным приближением будем делать перебор. Переберём все значения параметров p, P и q, Q от нуля до вот этих приближений, которые мы подобрали по графикам автокорреляционных функций. Всего таких моделей будет 60. Давайте просто настроим все их и посчитаем для каждой значение информационного критерия Акаике, выберем ту модель, у которой это значение будет минимальным. Это делает следующий кусок кода. Сезонная модель ARIMA настраивается с помощью функции SARIMAX из пакета statsmodels. Если эта функция у вас не будет работать, вам нужно будет обновить пакет statsmodels, так чтобы у вас была его версия выше седьмой. Начиная с восьмой версии эта функция уже работает. Перед тем как делать перебор, давайте отключим Warning'и. Дело в том, что не все модели, которые мы будем перебирать, можно хорошо настроить. Часть из них будут расходящимися, и про все такие модели Python нам будет выдавать предупреждение о том, что сходимости не произошло. Отключим вывод этих предупреждений и будем просто печатать параметры каждой модели, которая не сошлась. Прекрасно. Этот код выполняется довольно долго. 60 моделей перебирается в течение примерно 20 секунд. Вот есть несколько комбинаций значений параметров, при которых сходимости не произошло. Давайте посмотрим на лучшие модели. Самый лучший результат, то есть результат с минимальным значением информационного критерия Акаике показала вот такая модель: С p = 2, q = 1, P = 0 и Q = 1. Если мы посмотрим на всякий случай на несколько следующих по AIC'у моделей, мы увидим, что эти модели хоть и близки по значению информационного критерия, но параметров в них больше. Таким образом, мы действительно можем убедиться в том, что наша лучшая модель с самым маленьким Акаике не станет немножко незначительно хуже, если мы её чуть-чуть упростим. То есть никакая простая модель, более простая модель здесь в топ не входит. Модель на пятом месте — достаточно простая. Параметров в ней на один меньше, но значение информационного критерия Акаике у неё уже больше на три единицы. Давайте посмотрим на нашу лучшую модель. Используем метод Summary и посмотрим, что наша функция из statsmodels нам настроила. Мы получили модель SARIMAX с параметрами (2, 1, 1) (0, 1, 1). Сезонный период мы фиксировали и задали равным 12. Что ещё здесь есть полезного в этом выводе? Здесь есть значение статистики критерия Льюнг — Бокса при каком-то автоматически выбранном Q и соответствующий достигаемый уровень значимости — мы видим, что достигаемый уровень значимости достаточно большой. То есть, по всей видимости, остатки модели не автокоррелированы. Вот всё, что в этом блоке выводится, вся информация, она касается именно остатков. Ну что здесь ещё можно посмотреть интересного? Пожалуй, всё остальное не представляет большого интереса. Давайте сделаем визуальный анализ остатков нашей модели. Вот так они выглядят. На первый взгляд кажется, что всё в порядке. В этих остатках не видно никакого тренда, никакой сезонности, они довольно-таки похожи на шум. Вот так выглядит автокорреляционная функция этих остатков, и мы видим, что здесь нет никакой существенной структуры. Да? Некоторые лаги действительно значимо отличаются от нуля. Но эти лаги находятся достаточно далеко, то есть эти лаги — большие. Ну, естественно, что какие-то лаги, когда вы строите автокорреляционную функцию, могут отличаться от нуля чисто случайно, просто в силу эффекта множественной проверки гипотез. Поэтому об этом можно не беспокоиться. Тем более вот здесь отличие от нуля возникает при лаге 22 — это лаг не сезонный, следовательно, можно закрыть на это глаза. Проверим свойства остатков с помощью формальных критериев. Критерий Стьюдента проверяет гипотезу несмещённости и её не отвергает с достигаемым уровнем значимости примерно 0,26. Критерий Дики-Фуллера уверенно отвергает гипотезу нестационарности. Таким образом, мы получаем, что остатки не смещены, стационарны и неавтокоррелированы. Следовательно, наша модель достаточно хороша. Давайте посмотрим, насколько хорошо эта модель описывает данные. Вот здесь синяя линия — это исходный ряд. А красное — это то, что предсказывает нашу модель в этом месте. Мы видим, что действительно модель и данные достаточно похожи друг на друга. Ну всё. Давайте теперь построим прогноз. Прогноз делается с помощью функции predict. Не забудем применить обратное преобразование к преобразованию Бокса-Кокса, потому что модель, которую мы подобрали, она не знает о том, что исходный ряд был преобразован. И всё в ней подбирается на ряд после преобразования Бокса-Кокса. Вот. Ну и вот результат нашего прогнозирования. Здесь синее — это информация, которая у нас была. Красное — это прогноз на три года вперёд. Этот прогноз выглядит достаточно адекватным. Он передаёт то, что мы знаем о сезонности, глядя на предыдущий кусок ряда, и в нём, в принципе, нет тренда. Ну, это тоже можно понять, потому что тренд имел достаточно сложную структуру, и похоже, вот в тот момент, в который мы делаем предсказание, очень сложно знать заранее, будет ли средний уровень ряда дальше повышаться или понижаться. Прогноз достаточно адекватный.