Пошаговое руководство по статистике и Python для прогнозирования временных рядов
Вы когда-нибудь представляли себе предсказание будущего? Что ж, мы еще не достигли этого, но модели прогнозирования (с уровнем неопределенности) дают нам отличную ориентацию для более уверенного планирования нашего бизнеса, когда мы смотрим в будущее. В этом посте мы продемонстрируем подход к прогнозированию временных рядов продаж в автомобильной промышленности с использованием модели SARIMA.
Объяснение модели
SARIMA используется для нестационарных рядов, то есть, когда данные не колеблются вокруг одного и того же среднего значения, дисперсии и ковариации. Эта модель может определять тренд и сезонность, что делает ее очень важной. SARIMA состоит из других моделей прогнозирования:
AR: авторегрессивная модель (может быть простой, множественной или нелинейной регрессией).
MA: модель скользящих средних. В моделях скользящего среднего могут использоваться весовые коэффициенты, где наблюдения взвешиваются с помощью коэффициента обрезки (для самых старых данных в серии) и с более высоким весом для самых последних наблюдений.
Композиция AR и MA вместе несут модель ARMA, но эта модель используется только для стационарных рядов (среднее значение, постоянная дисперсия во времени).
Если ряд имеет тенденцию, необходимо использовать модель ARIMA.
ARIMA используется для нестационарных рядов. В этой модели шаг дифференцирования I (d) используется для устранения нестационарности.
Интегрированный элемент «I» для дифференцирования позволяет методу поддерживать временные ряды с трендом. Но все же эта модель не определяет сезонность.
Наконец, мы приходим к модели SARIMA, которая имеет сезонную корреляцию и может определять сезонность временного ряда. Теперь мы можем перейти к кодам!
Обработка данных
Мы собираемся использовать набор данных по продажам автомобилей, который можно скачать здесь.
import warnings
import itertools
import numpy as np
import matplotlib.pyplot as plt
warnings.filterwarnings("ignore")
plt.style.use('fivethirtyeight')
import pandas as pd
import statsmodels.api as sm
import matplotlib
matplotlib.rcParams['axes.labelsize'] = 14
matplotlib.rcParams['xtick.labelsize'] = 12
matplotlib.rcParams['ytick.labelsize'] = 12
matplotlib.rcParams['text.color'] = 'G'
df = pd.read_excel("Retail2.xlsx")
Этот шаг предназначен только для импорта библиотек, таких как numpy, pandas, matplotlib и statsmodels, то есть библиотеки, содержащей модель SARIMA и другие функции статистики. Эта часть кода используется для настройки диаграмм matplotlib.
Исходный набор данных и код немного сложны, но для упрощения загрузки файл, доступный для загрузки, имеет только столбцы даты и продажи, чтобы избежать предварительной обработки данных.
y = df.set_index(['Date']) y.head(5)
Команда set_index устанавливает дату столбца как индекс, а заголовок печатает первые 5 строк набора данных.
y.plot(figsize=(19, 4)) plt.show()
Анализируя график, мы можем заметить, что временной ряд имеет сезонный характер. Октябрь - пик продаж, по крайней мере, за последние 3 года. С годами также наблюдается тенденция к росту.
from pylab import rcParams rcParams['figure.figsize'] = 18, 8 decomposition = sm.tsa.seasonal_decompose(y, model='additive') fig = decomposition.plot() plt.show()
Используя команду «sm.tsa.seasonal_decompose» из библиотеки pylab, мы можем разложить временной ряд на три отдельных компонента: тренд, сезонность и шум.
SARIMA для прогнозирования временных рядов
Давайте использовать SARIMA. Обозначение модели SARIMA(p, d, q).(P,D,Q)m. Эти три параметра учитывают сезонность, тренд и шум в данных.
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 for SARIMA...')
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]))
Примеры параметра для SARIMA…
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)
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
ARIMA (0, 0, 0) x (0, 0, 1, 12) 12 - AIC: 410.521537786262
ARIMA (0, 0, 0) x (1, 0, 0, 12) 12 - AIC: 363.03322198787765
ARIMA (0, 0, 0) x (1, 0, 1, 12) 12 - AIC: 348.13731052896617
ARIMA (0, 0, 0) x (1, 1, 0, 12) 12 - AIC: 237.2069523573001
ARIMA (0, 0, 1) x (0, 1, 1, 12) 12 - AIC: 1005.9793011993983
ARIMA (0, 0, 1) x (1, 0, 0 , 12) 12 - AIC: 364.8331172721984
ARIMA (0, 0, 1) x (1, 0, 1, 12) 12 - AIC: 341.5095766512678
ARIMA (0, 0, 1) x (1 , 1, 0, 12) 12 - AIC: 239.13525566698246
ARIMA (0, 0, 1) x (1, 1, 1, 12) 12 - AIC: 223.43134436185036
Согласно Петерсону Т. (2014), AIC (информационный критерий Акаике) является оценкой относительного качества статистических моделей для заданного набора данных. Учитывая набор моделей для данных, AIC оценивает качество каждой модели относительно каждой из других моделей. Чем меньше значение AIC, тем лучше. Наши результаты показывают, что SARIMAX(0, 0, 1)x(1, 1, 1, 12) со значением AIC 223,43 - лучшая комбинация, поэтому мы должны считать это оптимальным вариантом.
mod = sm.tsa.statespace.SARIMAX(y,
order=(0, 0, 1),
seasonal_order=(1, 1, 1, 12),
enforce_stationarity=False,
enforce_invertibility=False)
results = mod.fit()
print(results.summary().tables[1])
в команде «mod = sm.tsa.statespace.SARIMAX» нам нужно настроить выбранную комбинацию.
results.plot_diagnostics(figsize=(18, 8)) plt.show()
С помощью приведенной выше диагностики мы можем визуализировать важную информацию в виде распределения и функции автокорреляции ACF (коррелограмма). Значения выше «0» имеют некоторую корреляцию с данными временного ряда. Значения, близкие к «1», демонстрируют наиболее сильную корреляцию.
pred = results.get_prediction(start=pd.to_datetime('2018-06-01'), dynamic=False)
pred_ci = pred.conf_int()
ax = y['2015':].plot(label='observed')
pred.predicted_mean.plot(ax=ax, label='One-step ahead Forecast', alpha=.7, figsize=(14, 4))
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('Retail_sold')
plt.legend()
plt.show()
Этот шаг заключается в сравнении истинных значений с предсказаниями прогноза. Наши прогнозы очень хорошо соответствуют истинным значениям. Команда «pred = results.get_prediction (start = pd.to_datetime (‘ 2018–06–01 ’)» определяет период, который вы бы прогнозировали при сравнении с истинными данными.
y_forecasted = pred.predicted_mean
y_truth = y['2018-06-01':]
mse = ((y_forecasted - y_truth) ** 2).mean()
print('The Mean Squared Error is {}'.format(round(mse, 2)))
print('The Root Mean Squared Error is {}'.format(round(np.sqrt(mse), 2)))
Среднеквадратичная ошибка 595,97
Среднеквадратичная ошибка 24,41
Замечания: как в MSE, так и в RMSE значения, близкие к нулю, лучше. Они являются мерой точности.
pred_uc = results.get_forecast(steps=12)
pred_ci = pred_uc.conf_int()
ax = y.plot(label='observed', figsize=(14, 4))
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('Sales')
plt.legend()
plt.show()
Здесь мы прогнозируем продажи на ближайшие 12 месяцев. Этот параметр можно изменить в строке кода «pred_uc = results.get_forecast (steps = 12)».
y_forecasted = pred.predicted_mean y_forecasted.head(12)
Этот шаг показывает прогнозируемые значения теста, который мы запускали ранее.
y_truth.head(12)
На этом шаге указываются истинные значения набора данных. Мы можем сравнить две приведенные выше серии, чтобы измерить точность модели.
pred_ci.head(24)
В таблице выше мы можем визуализировать нижнее и верхнее значения, которые модель указывает в качестве границ для прогнозирования.
forecast = pred_uc.predicted_mean forecast.head(12)
Наконец, у нас есть прогноз продаж на ближайшие 12 месяцев!
Мы можем заметить, что продажи со временем увеличиваются. Январь и февраль - худшие месяцы, как и в прошлые годы.
Это просто эксперимент со статистической моделью. Я действительно рекомендую попробовать повторяющиеся нейронные сети для прогнозирования, мой совет - просто используйте гораздо больше данных. Но это, безусловно, тема для другой статьи.
Ссылка:
Руководство по прогнозированию временных рядов с ARIMA в Python 3
