Прогнозируем временные ряды с ARIMA в Python 3

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

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

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

Что нужно

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

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

Для работы с данными я как всегда выбираю  Jupyter Notebook . Если у вас его нет или что-то не настроено, то идите по ссылке и начинайте от сюда.

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

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

[email protected]:~$ cd environments

 

[email protected]:~/environments$ . my_env/bin/activate

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

(my_env) [email protected]:~/environments$ mkdir ARIMA

 

(my_env) [email protected]:~/environments$ cd ARIMA

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

Ну а что насчет других библиотек, то тут вам в помощь pip. Мы должны установить pandas, statsmodels и пакет данных matplotlib:

(my_env) [email protected]:~/environments$ pip install pandas numpy statsmodels matplotlib

Ну теперь все готово. Продолжаем.

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

Чтобы начать работу просто запустите Jupyter Notebook:

(my_env) [email protected]:~/environments/ARIMA$ jupyter notebook

 Создавайте новый python file. New > Python 3   в выпадающем меню в правом верхнем углу:
Create a new Python 3 notebook

Сначала нужно импортировать все нужные библиотеки и сделаем мы это в самом начале документа:

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.

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

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

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

y = y[‘co2’].resample(‘MS’).mean()
y = y.fillna(y.bfill())
print(y)

Output
co2
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()

Figure 1: CO2 Levels Time Series

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

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

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

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

В общем-то, для прогнозирования временных времен 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 ) несезонные параметры и они описаны выше.

Термин s является периодичностью временного ряда (4 для квартальных периодов, 12 для годовых периодов и т. Д.).

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

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

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

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

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

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

Сначала создадим комбинации параметров, которые мы хотим оценить:

p = d = q = range(0, 2)
pdq = list(itertools.product(p, d, q))
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]))

Output
Examples 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)

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

Здесь аргумент order устанавливает параметры  (p, d, q), в то время как seasonal_order определяет параметры (P, D, Q, S).

После обучения каждой модели SARIMAX () выдает соответствующую оценку AIC:

warnings.filterwarnings(“ignore“)
for param in pdq:
for param_seasonal in seasonal_pdq:
try:
mod = sm.tsa.statespace.SARIMAX(y,
order=param,
seasonal_order=param_seaso
nal,
enforce_stationarity=False,
enforce_invertibility=False)
results = mod.fit()
print(‘ARIMA{}x{}12 – AIC:{}’.format(param, param_seasonal, results.aic))
except:
continue

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

Результат должен быть примерно таким:

Output
SARIMAX(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)
nal,
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, которая вызывается из results дает нам значительный объем информации,  о мы сосредоточим наше внимание на таблице коэффициентов.

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

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

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

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

Figure 2: Model Diagnostics

На этих моделях диагностики мы можем видеть следующее:

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

Шаг 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()

Figure 3: CO2 Levels Static Forecast

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

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)))

Output
The Mean Squared Error of our forecasts is 0.07

Результат – 0, 07. Он очень мал. Почти близок нулю.

Если у вас mse выдает 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()

Figure 4: CO2 Levels Dynamic Forecast

Ну и вычисляем MSE:

# извлекаем прогнозированные и верные значения нашего временного ряда
y_forecasted = pred_dynamic.predicted_mean
y_truth = y[‘1998-01-01‘:]
# вычисляем среднеквадратичную ошибку
mse = ((y_forecasted – y_truth) ** 2).mean()
print(‘The Mean Squared Error of our forecasts is {}‘.format(round(mse, 2)))

 

Output
The Mean Squared Error of our forecasts is 1.01

В это результат MSE у нас получился намного выше, в первом случае. Ну этого и следовало ожидать, так как здесь мы почти не используем исторические данные.

Шаг 7 – Подготовка и визуализация прогнозов

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

# получаем прогноз на 500 шагов
pred_uc = results.get_forecast(steps=500)
# получаем доверительные интервалы
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()

Figure 5: Time Series and Forecast of Future Values

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