Руководство по прогнозированию временных рядов с помощью ARIMA в Python 3

Введение

Временные ряды дают возможность прогнозировать будущие значения. Основываясь на предыдущих значениях, временные ряды можно использовать для прогнозирования тенденций в экономике, погоде и планировании пропускной способности, и это лишь некоторые из них. Специфические свойства данных временных рядов означают, что обычно требуются специализированные статистические методы.

В этом уроке мы будем стремиться производить надежные прогнозы временных рядов. Мы начнем с введения и обсуждения понятий автокорреляции, стационарности и сезонности, а затем перейдем к применению одного из наиболее часто используемых методов прогнозирования временных рядов, известного как ARIMA.

Один из методов, доступных в Python для моделирования и прогнозирования будущих точек временного ряда, известен как SARIMAX, что означает сезонные авторегрессивные интегрированные скользящие средние с экзогенными регрессорами. Здесь мы в первую очередь сосредоточимся на компоненте ARIMA, который используется для сопоставления данных временных рядов, чтобы лучше понять и спрогнозировать будущие точки временного ряда.

Предпосылки

В этом руководстве рассказывается, как выполнять анализ временных рядов на локальном рабочем столе или на удаленном сервере. Работа с большими наборами данных может потребовать большого объема памяти, поэтому в любом случае компьютеру потребуется не менее 2 ГБ памяти для выполнения некоторых расчетов в этом руководстве.

Чтобы максимально использовать этот учебник, может быть полезно некоторое знакомство с временными рядами и статистикой.

В этом уроке мы будем использовать Jupyter Notebook для работы с данными. Если у вас его еще нет, вы должны следовать нашему руководству, чтобы установить и настроить Jupyter Notebook для Python 3.

Шаг 1 — Установка пакетов

Чтобы настроить нашу среду для прогнозирования временных рядов, давайте сначала перейдем к нашей серверной среде программирования:

cd environments

. my_env/bin/activate

Отсюда давайте создадим новый каталог для нашего проекта. Мы назовем его ARIMA, а затем перейдем в каталог. Если вы называете проект другим именем, обязательно замените свое имя на ARIMA во всем руководстве.

mkdir ARIMA
cd ARIMA

Для этого руководства потребуются warnings, itertools, pandas, numpy, matplotlib и библиотеки statsmodels. Библиотеки warnings и itertools входят в стандартный набор библиотек Python, поэтому вам не нужно их устанавливать.

Как и в случае с другими пакетами Python, мы можем установить эти требования с помощью pip.

pip install pandas numpy statsmodels matplotlib

На данный момент мы готовы начать работу с установленными пакетами.

Шаг 2 — Импорт пакетов и загрузка данных

Чтобы начать работу с нашими данными, мы запустим Jupyter Notebook:

jupyter notebook

Чтобы создать новый файл записной книжки, выберите «Создать» > «Python 3» в правом верхнем раскрывающемся меню:

Это откроет блокнот.

Рекомендуется начать с импорта необходимых библиотек в верхнюю часть блокнота:

import warnings
import itertools
import pandas as pd
import numpy as np
import statsmodels.api as sm
import matplotlib.pyplot as plt
plt.style.use('fivethirtyeight')

Мы также определили для наших графиков стиль matplotlib, равный FivethirtyEight.

Мы будем работать с набором данных под названием «Атмосферный CO2 из непрерывных проб воздуха в обсерватории Мауна-Лоа, Гавайи, США», в котором были собраны пробы CO2 с марта 1958 г. по декабрь 2001 г. Мы можем ввести эти данные следующим образом:

data = sm.datasets.co2.load_pandas()
y = data.data

Давайте немного предварительно обработаем наши данные, прежде чем двигаться дальше. Еженедельные данные могут быть сложными для работы, поскольку это более короткий промежуток времени, поэтому давайте вместо этого будем использовать среднемесячные значения. Мы сделаем преобразование с помощью функции resample. Для простоты мы также можем использовать функцию fillna(), чтобы убедиться, что в нашем временном ряду нет пропущенных значений.

# The 'MS' string groups the data in buckets by start of the month
y = y['co2'].resample('MS').mean()

# The term bfill means that we use the value before filling in missing values
y = y.fillna(y.bfill())

print(y)

Outputco2
1958-03-01  316.100000
1958-04-01  317.200000
1958-05-01  317.433333
...
2001-11-01  369.375000
2001-12-01  371.020000

Давайте рассмотрим этот временной ряд как визуализацию данных:

y.plot(figsize=(15, 6))
plt.show()

Когда мы наносим данные на график, появляются некоторые различимые закономерности. Временной ряд имеет очевидную сезонность, а также общую тенденцию к увеличению.

Чтобы узнать больше о предварительной обработке временных рядов, обратитесь к «Руководству по визуализации временных рядов с помощью Python 3», где приведенные выше шаги описаны более подробно.

Теперь, когда мы преобразовали и изучили наши данные, давайте перейдем к прогнозированию временных рядов с помощью ARIMA.

Шаг 3 — Модель временных рядов ARIMA

Один из наиболее распространенных методов, используемых в прогнозировании временных рядов, известен как модель ARIMA, что означает AutoregRessive Integrated Moving Average. ARIMA — это модель, которую можно адаптировать к данным временных рядов, чтобы лучше понять или предсказать будущие точки ряда.

Существует три различных целых числа (p, d, q), которые используются для параметризации моделей ARIMA. Из-за этого модели ARIMA обозначаются нотацией ARIMA(p, d, q). Вместе эти три параметра учитывают сезонность, тенденцию и шум в наборах данных:

p — это авторегрессивная часть модели. Это позволяет нам включить эффект прошлых значений в нашу модель. Интуитивно это было бы похоже на заявление о том, что завтра будет тепло, если последние 3 дня было тепло.
d — это интегрированная часть модели. Это включает в себя термины в модели, которые включают количество различий (т. е. количество прошлых моментов времени, которые нужно вычесть из текущего значения) для применения к временному ряду. Интуитивно это было бы похоже на заявление о том, что завтра будет такая же температура, если разница температур за последние три дня была очень небольшой.
q — это скользящая средняя часть модели. Это позволяет нам установить ошибку нашей модели как линейную комбинацию значений ошибок, наблюдаемых в предыдущие моменты времени в прошлом.

При работе с сезонными эффектами мы используем сезонный ARIMA, который обозначается как ARIMA(p,d,q)(P,D,Q)s. Здесь (p, d, q) — несезонные параметры, описанные выше, а (P, D, Q) следуют тому же определению, но применяются к сезонным параметрам. компонент временного ряда. Термин s — это периодичность временного ряда (4 для квартальных периодов, 12 для годовых периодов и т. д.).

Сезонный метод ARIMA может показаться сложным из-за задействования множества параметров настройки. В следующем разделе мы опишем, как автоматизировать процесс определения оптимального набора параметров для сезонной модели временных рядов ARIMA.

Шаг 4 — Выбор параметров для модели временных рядов ARIMA

При поиске соответствия данных временных рядов сезонной модели ARIMA наша первая цель — найти значения ARIMA(p,d,q)(P,D,Q)s, которые оптимизируют показатель интерес. Существует множество руководств и лучших практик для достижения этой цели, однако правильная параметризация моделей ARIMA может быть кропотливым ручным процессом, требующим знаний в предметной области и времени. Другие языки статистического программирования, такие как R, предоставляют автоматизированные способы решения этой проблемы, но их еще предстоит перенести на Python. В этом разделе мы решим эту проблему, написав код Python для программного выбора оптимальных значений параметров для нашей модели временных рядов ARIMA(p,d,q)(P,D,Q)s.

Мы будем использовать «поиск по сетке» для многократного изучения различных комбинаций параметров. Для каждой комбинации параметров мы подбираем новую сезонную модель ARIMA с помощью функции SARIMAX() из statsmodels и оцените его общее качество. После того, как мы изучим весь спектр параметров, нашим оптимальным набором параметров будет тот, который обеспечивает наилучшую производительность по интересующим нас критериям. Давайте начнем с создания различных комбинаций параметров, которые мы хотим оценить:

# Define the p, d and q parameters to take any value between 0 and 2
p = d = q = range(0, 2)

# Generate all different combinations of p, q and q triplets
pdq = list(itertools.product(p, d, q))

# Generate all different combinations of seasonal p, q and q triplets
seasonal_pdq = [(x[0], x[1], x[2], 12) for x in list(itertools.product(p, d, q))]

print('Examples of parameter combinations for Seasonal ARIMA...')
print('SARIMAX: {} x {}'.format(pdq[1], seasonal_pdq[1]))
print('SARIMAX: {} x {}'.format(pdq[1], seasonal_pdq[2]))
print('SARIMAX: {} x {}'.format(pdq[2], seasonal_pdq[3]))
print('SARIMAX: {} x {}'.format(pdq[2], seasonal_pdq[4]))

OutputExamples of parameter combinations for Seasonal ARIMA...
SARIMAX: (0, 0, 1) x (0, 0, 1, 12)
SARIMAX: (0, 0, 1) x (0, 1, 0, 12)
SARIMAX: (0, 1, 0) x (0, 1, 1, 12)
SARIMAX: (0, 1, 0) x (1, 0, 0, 12)

Теперь мы можем использовать триплеты параметров, определенные выше, для автоматизации процесса обучения и оценки моделей ARIMA для различных комбинаций. В статистике и машинном обучении этот процесс известен как поиск по сетке (или оптимизация гиперпараметров) для выбора модели.

При оценке и сравнении статистических моделей, оснащенных различными параметрами, каждую из них можно ранжировать относительно друг друга в зависимости от того, насколько хорошо она соответствует данным или ее способности точно предсказывать будущие точки данных. Мы будем использовать значение AIC (информационный критерий Akaike), которое удобно возвращать с моделями ARIMA, подобранными с помощью statsmodels. AIC измеряет, насколько хорошо модель соответствует данным, принимая во внимание общую сложность модели. Модель, которая очень хорошо соответствует данным при использовании большого количества функций, получит более высокий балл AIC, чем модель, которая использует меньше функций для достижения такого же качества соответствия. Поэтому нас интересует поиск модели, дающей наименьшее значение AIC.

В приведенном ниже фрагменте кода перебираются комбинации параметров и используется функция SARIMAX из statsmodels для соответствия соответствующей сезонной модели ARIMA. Здесь аргумент order указывает параметры (p, d, q), а аргумент seasonal_order указывает (P, D , Q, S) сезонный компонент модели Seasonal ARIMA. После подгонки каждой модели SARIMAX() код выводит соответствующую оценку AIC.

warnings.filterwarnings("ignore") # specify to ignore warning messages

for param in pdq:
    for param_seasonal in seasonal_pdq:
        try:
            mod = sm.tsa.statespace.SARIMAX(y,
                                            order=param,
                                            seasonal_order=param_seasonal,
                                            enforce_stationarity=False,
                                            enforce_invertibility=False)

            results = mod.fit()

            print('ARIMA{}x{}12 - AIC:{}'.format(param, param_seasonal, results.aic))
        except:
            continue

Поскольку некоторые комбинации параметров могут привести к неправильным числовым спецификациям, мы явно отключили предупреждающие сообщения, чтобы избежать перегрузки предупреждающих сообщений. Эти неправильные спецификации также могут привести к ошибкам и вызвать исключение, поэтому мы обязательно перехватываем эти исключения и игнорируем комбинации параметров, которые вызывают эти проблемы.

Приведенный выше код должен дать следующие результаты, это может занять некоторое время:

OutputSARIMAX(0, 0, 0)x(0, 0, 1, 12) - AIC:6787.3436240402125
SARIMAX(0, 0, 0)x(0, 1, 1, 12) - AIC:1596.711172764114
SARIMAX(0, 0, 0)x(1, 0, 0, 12) - AIC:1058.9388921320026
SARIMAX(0, 0, 0)x(1, 0, 1, 12) - AIC:1056.2878315690562
SARIMAX(0, 0, 0)x(1, 1, 0, 12) - AIC:1361.6578978064144
SARIMAX(0, 0, 0)x(1, 1, 1, 12) - AIC:1044.7647912940095
...
...
...
SARIMAX(1, 1, 1)x(1, 0, 0, 12) - AIC:576.8647112294245
SARIMAX(1, 1, 1)x(1, 0, 1, 12) - AIC:327.9049123596742
SARIMAX(1, 1, 1)x(1, 1, 0, 12) - AIC:444.12436865161305
SARIMAX(1, 1, 1)x(1, 1, 1, 12) - AIC:277.7801413828764

Вывод нашего кода показывает, что SARIMAX(1, 1, 1)x(1, 1, 1, 12) дает наименьшее значение AIC, равное 277,78. Поэтому мы должны считать этот вариант оптимальным из всех рассмотренных нами моделей.

Шаг 5 — Подгонка модели временных рядов ARIMA

Используя поиск по сетке, мы определили набор параметров, который обеспечивает наилучшее соответствие модели нашим данным временного ряда. Мы можем перейти к более глубокому анализу этой конкретной модели.

Мы начнем с добавления оптимальных значений параметров в новую модель SARIMAX:

mod = sm.tsa.statespace.SARIMAX(y,
                                order=(1, 1, 1),
                                seasonal_order=(1, 1, 1, 12),
                                enforce_stationarity=False,
                                enforce_invertibility=False)

results = mod.fit()

print(results.summary().tables[1])

Output==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
ar.L1          0.3182      0.092      3.443      0.001       0.137       0.499
ma.L1         -0.6255      0.077     -8.165      0.000      -0.776      -0.475
ar.S.L12       0.0010      0.001      1.732      0.083      -0.000       0.002
ma.S.L12      -0.8769      0.026    -33.811      0.000      -0.928      -0.826
sigma2         0.0972      0.004     22.634      0.000       0.089       0.106
==============================================================================

Атрибут summary, являющийся результатом вывода SARIMAX, возвращает значительный объем информации, но мы сосредоточим наше внимание на таблице коэффициентов. Столбец coef показывает вес (то есть важность) каждой функции и то, как каждая из них влияет на временной ряд. Столбец P>|z| информирует нас о значимости каждого веса функции. Здесь каждый вес имеет p-значение меньше или близко к 0,05, поэтому разумно сохранить их все в нашей модели.

При подборе сезонных моделей ARIMA (и любых других моделей в этом отношении) важно запустить диагностику модели, чтобы убедиться, что ни одно из допущений, сделанных моделью, не было нарушено. Объект plot_diagnostics позволяет нам быстро генерировать диагностику модели и исследовать любое необычное поведение.

results.plot_diagnostics(figsize=(15, 12))
plt.show()

Нашей главной задачей является обеспечение того, чтобы остатки нашей модели были некоррелированными и нормально распределенными с нулевым средним. Если сезонная модель ARIMA не удовлетворяет этим свойствам, это хороший признак того, что ее можно улучшить.

В этом случае наша диагностика модели предполагает, что остатки модели нормально распределяются на основе следующего:

На верхнем правом графике мы видим, что красная строка KDE следует за линией N(0,1) (где N(0,1)) — это стандартное обозначение нормального распределения со средним значением 0 и стандартным отклонением 1). Это хороший показатель нормального распределения остатков.
График qq в левом нижнем углу показывает, что упорядоченное распределение остатков (синие точки) соответствует линейному тренду выборок, взятых из стандартного нормального распределения с N(0, 1). Опять же, это убедительный признак того, что остатки распределены нормально.
Остатки во времени (верхний левый график) не имеют явной сезонности и кажутся белым шумом. Это подтверждается графиком автокорреляции (то есть коррелограммой) в правом нижнем углу, который показывает, что остатки временного ряда имеют низкую корреляцию с запаздывающими версиями самого себя.

Эти наблюдения приводят нас к выводу, что наша модель дает удовлетворительное соответствие, которое может помочь нам понять наши данные временного ряда и прогнозировать будущие значения.

Хотя у нас есть удовлетворительное соответствие, некоторые параметры нашей сезонной модели ARIMA можно изменить, чтобы улучшить соответствие нашей модели. Например, наш поиск по сетке рассматривал только ограниченный набор комбинаций параметров, поэтому мы можем найти лучшие модели, если расширим поиск по сетке.

Шаг 6 — Проверка прогнозов

Мы получили модель для нашего временного ряда, которую теперь можно использовать для получения прогнозов. Мы начинаем со сравнения прогнозируемых значений с реальными значениями временного ряда, что поможет нам понять точность наших прогнозов. Атрибуты get_prediction() и conf_int() позволяют нам получать значения и соответствующие доверительные интервалы для прогнозов временных рядов.

pred = results.get_prediction(start=pd.to_datetime('1998-01-01'), dynamic=False)
pred_ci = pred.conf_int()

Приведенный выше код требует, чтобы прогнозы начинались с января 1998 года.

Аргумент dynamic=False гарантирует, что мы создаем прогнозы на один шаг вперед, что означает, что прогнозы в каждой точке генерируются с использованием полной истории до этой точки.

Мы можем построить реальные и прогнозируемые значения временного ряда CO2, чтобы оценить, насколько хорошо мы справились. Обратите внимание, как мы увеличили конец временного ряда, разрезав индекс даты.

ax = y['1990':].plot(label='observed')
pred.predicted_mean.plot(ax=ax, label='One-step ahead Forecast', alpha=.7)

ax.fill_between(pred_ci.index,
                pred_ci.iloc[:, 0],
                pred_ci.iloc[:, 1], color='k', alpha=.2)

ax.set_xlabel('Date')
ax.set_ylabel('CO2 Levels')
plt.legend()

plt.show()

В целом, наши прогнозы очень хорошо совпадают с истинными значениями, демонстрируя общую тенденцию к увеличению.

Также полезно количественно оценить точность наших прогнозов. Мы будем использовать MSE (среднеквадратическую ошибку), которая суммирует среднюю ошибку наших прогнозов. Для каждого предсказанного значения мы вычисляем его расстояние до истинного значения и возводим результат в квадрат. Результаты должны быть возведены в квадрат, чтобы положительные/отрицательные различия не компенсировали друг друга, когда мы вычисляем общее среднее значение.

y_forecasted = pred.predicted_mean
y_truth = y['1998-01-01':]

# Compute the mean square error
mse = ((y_forecasted - y_truth) ** 2).mean()
print('The Mean Squared Error of our forecasts is {}'.format(round(mse, 2)))

OutputThe Mean Squared Error of our forecasts is 0.07

СКО наших прогнозов на один шаг вперед дает значение 0,07, что очень мало, поскольку оно близко к 0. СКО, равное 0, означает, что оценщик предсказывает наблюдения параметра с идеальной точностью. , что было бы идеальным сценарием, но обычно это невозможно.

Однако лучшее представление о нашей истинной предсказательной способности можно получить, используя динамические прогнозы. В этом случае мы используем информацию из временного ряда только до определенного момента, а после этого прогнозы генерируются с использованием значений из предыдущих прогнозируемых моментов времени.

В приведенном ниже фрагменте кода мы указываем, чтобы начать вычисление динамических прогнозов и доверительных интервалов с января 1998 года.

pred_dynamic = results.get_prediction(start=pd.to_datetime('1998-01-01'), dynamic=True, full_results=True)
pred_dynamic_ci = pred_dynamic.conf_int()

Нанося на график наблюдаемые и прогнозируемые значения временного ряда, мы видим, что общие прогнозы точны даже при использовании динамических прогнозов. Все прогнозируемые значения (красная линия) довольно близко соответствуют действительности (синяя линия) и находятся в пределах доверительных интервалов нашего прогноза.

ax = y['1990':].plot(label='observed', figsize=(20, 15))
pred_dynamic.predicted_mean.plot(label='Dynamic Forecast', ax=ax)

ax.fill_between(pred_dynamic_ci.index,
                pred_dynamic_ci.iloc[:, 0],
                pred_dynamic_ci.iloc[:, 1], color='k', alpha=.25)

ax.fill_betweenx(ax.get_ylim(), pd.to_datetime('1998-01-01'), y.index[-1],
                 alpha=.1, zorder=-1)

ax.set_xlabel('Date')
ax.set_ylabel('CO2 Levels')

plt.legend()
plt.show()

Еще раз, мы количественно оцениваем прогностическую эффективность наших прогнозов, вычисляя MSE:

# Extract the predicted and true values of our time series
y_forecasted = pred_dynamic.predicted_mean
y_truth = y['1998-01-01':]

# Compute the mean square error
mse = ((y_forecasted - y_truth) ** 2).mean()
print('The Mean Squared Error of our forecasts is {}'.format(round(mse, 2)))

OutputThe Mean Squared Error of our forecasts is 1.01

Прогнозные значения, полученные из динамических прогнозов, дают MSE 1,01. Это немного выше, чем на один шаг вперед, чего и следовало ожидать, учитывая, что мы полагаемся на менее исторические данные из временных рядов.

Как прогнозы на один шаг вперед, так и динамические прогнозы подтверждают правильность этой модели временных рядов. Тем не менее, большая часть интереса к прогнозированию временных рядов заключается в возможности прогнозировать будущие значения заблаговременно.

Шаг 7 — Создание и визуализация прогнозов

На последнем этапе этого руководства мы опишем, как использовать нашу сезонную модель временных рядов ARIMA для прогнозирования будущих значений. Атрибут get_forecast() нашего объекта временных рядов может вычислять прогнозируемые значения на заданное количество шагов вперед.

# Get forecast 500 steps ahead in future
pred_uc = results.get_forecast(steps=500)

# Get confidence intervals of forecasts
pred_ci = pred_uc.conf_int()

Мы можем использовать выходные данные этого кода для построения временных рядов и прогнозов их будущих значений.

ax = y.plot(label='observed', figsize=(20, 15))
pred_uc.predicted_mean.plot(ax=ax, label='Forecast')
ax.fill_between(pred_ci.index,
                pred_ci.iloc[:, 0],
                pred_ci.iloc[:, 1], color='k', alpha=.25)
ax.set_xlabel('Date')
ax.set_ylabel('CO2 Levels')

plt.legend()
plt.show()

И прогнозы, и связанный с ними доверительный интервал, которые мы сгенерировали, теперь можно использовать для дальнейшего понимания временных рядов и прогнозирования того, чего ожидать. Наши прогнозы показывают, что временной ряд, как ожидается, будет продолжать расти устойчивыми темпами.

По мере того, как мы прогнозируем будущее, для нас естественно становиться менее уверенными в наших ценностях. Это отражено в доверительных интервалах, созданных нашей моделью, которые увеличиваются по мере того, как мы продвигаемся дальше в будущее.

Заключение

В этом руководстве мы описали, как реализовать сезонную модель ARIMA в Python. Мы широко использовали библиотеки pandas и statsmodels и показали, как запускать диагностику модели, а также как создавать прогнозы временных рядов CO2.

Вот еще несколько вещей, которые вы можете попробовать:

Измените дату начала ваших динамических прогнозов, чтобы увидеть, как это повлияет на общее качество ваших прогнозов.
Попробуйте больше комбинаций параметров, чтобы увидеть, сможете ли вы улучшить качество подгонки вашей модели.
Выберите другой показатель, чтобы выбрать лучшую модель. Например, мы использовали показатель AIC, чтобы найти наилучшую модель, но вместо этого вы могли попытаться оптимизировать среднеквадратичную ошибку вне выборки.

Для большей практики вы также можете попробовать загрузить другой набор данных временных рядов, чтобы получить свои собственные прогнозы.