Временные ряды позволяют прогнозировать значения. На основании предыдущих значений временных рядов можно спрогнозировать тенденции в экономике и погоде или спланировать пропускную способность. Ввиду особенных свойств данных временных рядов для работы с ними применяются специализированные статистические методы и подходы.
Данное руководство поможет создать прогнозы временных рядов, ознакомит вас с понятиями автокорреляции, стационарности и сезонности, а также научит пользоваться инструментом для прогнозирования временных рядов ARIMA.
Один из методов моделирования и прогнозирования будущих точек временного ряда в Python называется SARIMAX. Данное руководство сосредоточено на модели ARIMA, которая используется для обработки данных временных рядов, что позволяет лучше понять и спрогнозировать будущие точки временного ряда.
Требования
- Данное руководство можно выполнить как на локальной, так и на удалённой машине.
- Для работы с объёмными данными потребуется много памяти (не меньше 2 Гб).
- Базовые навыки работы с временными рядами и статистикой.
- Предварительно установленный инструмент Jupyter Notebook (инструкции по установке можно найти ).
1: Установка пакетов
Разверните среду программирования:
cd environments
. my_env/bin/activate
Создайте в ней каталог для нового проекта. В руководстве он называется ARIMA. Если вы выбрали другое имя каталога, откорректируйте команды. Откройте каталог.
mkdir ARIMA
cd ARIMA
В дальнейшей работе вам понадобятся библиотеки warnings, itertools, pandas, numpy, matplotlib и statsmodels. Warnings и itertools включены в Python по умолчанию, остальные библиотеки (а также их зависимости) нужно установить.
pip install pandas numpy statsmodels matplotlib
2: Импорт пакетов и загрузка данных
Чтобы начать работу с данными, запустите Jupyter Notebook:
jupyter notebook
Чтобы создать новый документ, откройте New → 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).
Для работы мы используем набор данных «Atmospheric CO2 from Continuous Air Samples at Mauna Loa Observatory, Hawaii, U.S.A.», в котором находятся данные о CO2 с марта 1958 года по декабрь 2001 года. Загрузите эти данные в документ:
data = sm.datasets.co2.load_pandas()
y = data.data
Прежде чем продолжить, предварительно обработайте данные. С данными за неделю сложнее работать, так как это более короткий промежуток времени; поэтому лучше использовать средние данные за месяц. Преобразуйте данные с помощью функции resample. Также можно использовать функцию fillna (), чтобы заполнить пропущенные значения во временных рядах.
# 'MS' группирует месячные данные
y = y['co2'].resample('MS').mean()
# bfill значит, что нужно использовать значение до заполнения пропущенных значений
y = y.fillna(y.bfill())
print(y)
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()
При построении графика на основе этих данных можно заметить некоторые явные шаблоны. Временные ряды имеют очевидную сезонность и общий тренд на увеличение.
3: Модель ARIMA
Модель ARIMA (AutoregRessive Integrated Moving Average) – один из наиболее распространённых методов анализа и прогнозирования временных рядов. Эта модель позволяет обработать данные временного ряда, чтобы лучше понять этот ряд или предсказать его развитие.
ARIMA использует три основных параметра (p, d, q), которые выражаются целыми числами. Потому модель также записывается как ARIMA(p, d, q). Вместе эти три параметра учитывают сезонность, тенденцию и шум в наборах данных:
- p – порядок авторегрессии (AR), который позволяет добавить предыдущие значения временного ряда. Этот параметр можно проиллюстрировать утверждением «завтра, вероятно, будет тепло, если в последние три дня было тепло».
- d – порядок интегрирования (I; т. е. порядок разностей исходного временного ряда). Он добавляет в модель понятия разности временных рядов (определяет количество прошлых временных точек, которые нужно вычесть из текущего значения). Этот параметр иллюстрирует такое утверждение: «завтра, вероятно, будет такая же температура, если разница в температуре за последние три дня была очень мала».
- q – порядок скользящего среднего (MA), который позволяет установить погрешность модели как линейную комбинацию наблюдавшихся ранее значений ошибок.
Для отслеживания сезонности используется сезонная модель 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).
Чтобы итерировать различные комбинации параметров, используйте сеточный поиск.
Для каждой комбинации параметров функция SARIMAX () из модуля statsmodels может подобрать новую сезонную модель ARIMA и оценить ее общее качество. Оптимальным набором параметров будет тот, в котором нужные критерии наиболее производительны. Для начала сгенерируйте различные комбинации параметров:
# Определите p, d и q в диапазоне 0-2
p = d = q = range(0, 2)
# Сгенерируйте различные комбинации p, q и q
pdq = list(itertools.product(p, d, q))
# Сгенерируйте комбинации сезонных параметров p, q и 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]))
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)
Теперь можно использовать определенные выше триплеты параметров для автоматизации процесса оценки моделей ARIMA по различным комбинациям. В статистике и машинном обучении этот процесс известен как поиск по сетке параметров (сетчатый поиск, или оптимизация гиперпараметров).
При оценке и сравнении статистических моделей, соответствующих различным параметрам, учитывается, насколько та или иная модель соответствует данным и насколько точно она способна прогнозировать будущие точки данных. Используйте значение AIC (Akaike Information Criterion), которое подходит для работы с моделями ARIMA на основе statsmodels. AIC оценивает, насколько хорошо модель соответствует данным, принимая во внимание общую сложность модели. Чем меньше функций использует модель, чтобы достичь соответствия данным, тем выше её показатель AIC. Поэтому нужно найти модель с наименьшим значением AIC.
Предложенный ниже код итерирует комбинации параметров и использует функцию SARIMAX из statsmodels, чтобы найти соответствия с сезонной моделью ARIMA. Аргумент order определяет параметры (p, d, q), а аргумент seasonal_order указывает сезонный компонент модели ARIMA (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_seasonal,
enforce_stationarity=False,
enforce_invertibility=False)
results = mod.fit()
print('ARIMA{}x{}12 - AIC:{}'.format(param, param_seasonal, results.aic))
except:
continue
Поскольку некоторые комбинации параметров могут выдавать неточный результат, нужно явно отключить предупреждающие сообщения, чтобы избежать перегрузки предупреждений. Эти неточности также могут приводить к ошибкам, поэтому комбинации параметров, которые вызывают подобные проблемы, нужно игнорировать.
Вышеприведённый код выведет такой результат:
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),
enforce_stationarity=False,
enforce_invertibility=False)
results = mod.fit()
print(results.summary().tables[1])
==============================================================================
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 возвращает много информации, но мы сосредоточим наше внимание на таблице коэффициентов. Столбец 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) . Это хороший признак того, что остатки нормально распределены.
- График q-q в левом нижнем углу показывает, что упорядоченное распределение остатков (синие точки) следует линейному тренду выборок, взятых из стандартного распределения 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 (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)))
The Mean Squared Error of our forecasts is 0.07
MSE прогнозов на один шаг вперед дает значение 0,07 (это очень низкое значение, так как оно близко к 0). 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()
Чтобы узнать точность прогноза, вычислите 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)))
The Mean Squared Error of our forecasts is 1.01
Прогнозируемые данные вернули MSE 1.01. Это немного больше, чем в предыдущем разделе (чего следовало ожидать, учитывая, что здесь используются менее точные данные временных рядов).
Как пошаговый прогноз, так и динамические прогнозы подтверждают, что эта модель временного ряда работает. Теперь можно попробовать спрогнозировать будущие значения ряда.
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()
И сгенерированные прогнозы, и связанный с ними интервал теперь можно использовать для дальнейшего анализа и прогнозирования временных рядов. Полученные данные показывают, что временные ряды будут продолжать стабильный рост.
Конечно, чем дальше строится прогноз, тем менее точны его значения. Это отражается на интервалах, генерируемых моделью (чем дальше прогноз, тем больше интервал).
Заключение
Теперь вы умеете работать с данными временных рядов, анализировать их и строить прогнозы с помощью сезонной модели ARIMA в Python.
Чтобы потренироваться, попробуйте загрузить другой временной ряд и выполнить руководство с новыми данными.
Специально для сайта ITWORLD.UZ. Новость взята с сайта www.8host.com