Анализ временных рядов для машинного обучения
Тренд, выбросы, стационарность, сезонность
Резюме
В описательной статистике временной ряд определяется как набор случайных величин, упорядоченных по времени. Временные ряды изучаются как для интерпретации явления, определения компонентов тренда, цикличности, сезонности, так и для прогнозирования его будущих значений. Я думаю, что они являются лучшим примером соединения между областью экономики и науки о данных (цены на акции, экономические циклы, бюджеты и денежные потоки,…).
В этой статье я шаг за шагом объясню стандартный подход к анализу временных рядов и представлю некоторые полезные инструменты (код Python), которые можно легко использовать в других подобных случаях (просто скопируйте, вставьте и запустите). Я рассмотрю каждую строку кода с комментариями, чтобы вы могли легко воспроизвести этот пример (ссылка на полный код ниже).
Мы будем использовать набор данных из конкурса Kaggle «Прогноз будущих продаж» (ссылка ниже), в котором вам будут предоставлены ежедневные исторические данные о продажах, а задача состоит в том, чтобы спрогнозировать общее количество проданных продуктов. Набор данных представляет собой интересный временной ряд, поскольку он очень похож на варианты использования, которые можно найти в реальном мире, поскольку мы знаем, что ежедневные продажи любого продукта никогда не бывают стационарными и всегда сильно зависят от сезонности.
При этом основная цель данного руководства - понять основные этапы анализа временных рядов перед тем, как приступить к разработке и тестированию моделей для прогнозирования (обратите внимание, что это статья предполагает базовые знания темы, поэтому я не буду вдаваться в определения, но вставлю все полезные гиперссылки). В частности, мы пройдем через:
- Анализ тренда, чтобы определить, является ли он линейным или нет, поскольку большинство моделей требуют эту информацию в качестве входных данных.
- Обнаружение выбросов, чтобы понять, как их обнаруживать и обрабатывать
- Тест стационарности, чтобы понять, можем ли мы предположить, что временной ряд является стационарным или нет
- Анализ сезонности для определения наилучшего сезонного параметра для использования при моделировании.
Полный код (Github):
Набор данных (Kaggle):
Прогнозирование будущих продаж
Финальный проект курса Coursera« Как победить в конкурсе по науке о данных www.kaggle.com»
Настройка
Прежде всего, мы импортируем следующие библиотеки
## For data import pandas as pd import numpy as np ## For plotting import matplotlib.pyplot as plt ## For outliers detection from sklearn import preprocessing, svm ## For stationarity test and decomposition import statsmodels.tsa.api as smt import statsmodels.api as sm
Затем мы прочитаем данные в фрейм данных pandas
dtf = pd.read_csv('data.csv')
dtf.head()
В исходном наборе данных есть разные столбцы, однако для целей этого руководства нам нужен только следующий столбец: дата и количество проданных продуктов (item_cnt_day). Другими словами, мы будем создавать серию панд (называемую продажи) с ежедневным индексом периодичности datetime, используя только дневную сумму продаж.
## format datetime column
dtf["date"] = pd.to_datetime(dtf['date'], format='%d.%m.%Y')
## create time series
ts = dtf.groupby("date")["item_cnt_day"].sum().rename("sales")
ts.head()
ts.tail()
Таким образом, временной ряд варьируется от 2013–01–01 до 2015–10–31, он содержит 1034 наблюдения, среднее значение 3528 и стандартное отклонение 1585. Это выглядит так:
ts.plot()
Давайте начнем, ладно?
Анализ тенденций
Цель этого анализа - понять, есть ли в данных тренд и является ли эта закономерность линейной или нет. Лучший инструмент для этой работы - визуализация.
Давайте напишем функцию, которая поможет нам понять тенденцию и движения временного ряда. Мы хотим видеть на графике некоторую скользящую статистику, такую как:
- M скользящее среднее: невзвешенное среднее предыдущих n данных (также называемое скользящим средним).
- Полосы Боллинджера: верхняя полоса в k раз на n -периодное стандартное отклонение выше скользящей средней и нижняя полоса в k раз на стандартное отклонение N ниже. скользящая средняя.
'''
Plot ts with rolling mean and 95% confidence interval with rolling std.
:parameter
:param ts: pandas Series
:param window: num - for rolling stats
:param plot_ma: bool - whether plot moving average
:param plot_intervals: bool - whether plot upper and lower bounds
'''
def plot_ts(ts, plot_ma=True, plot_intervals=True, window=30,
figsize=(15,5)):
rolling_mean = ts.rolling(window=window).mean()
rolling_std = ts.rolling(window=window).std()
plt.figure(figsize=figsize)
plt.title(ts.name)
plt.plot(ts[window:], label='Actual values', color="black")
if plot_ma:
plt.plot(rolling_mean, 'g', label='MA'+str(window),
color="red")
if plot_intervals:
lower_bound = rolling_mean - (1.96 * rolling_std)
upper_bound = rolling_mean + (1.96 * rolling_std)
plt.fill_between(x=ts.index, y1=lower_bound, y2=upper_bound,
color='lightskyblue', alpha=0.4)
plt.legend(loc='best')
plt.grid(True)
plt.show()
Когда набор данных имеет хотя бы полный год наблюдения, я всегда начинаю с скользящего окна 30 дней:
plot_ts(ts, window=30)
Глядя на красную линию на графике, вы можете легко заметить закономерность: временной ряд следует линейному нисходящему тренду с сильными сезонными пиками каждый январь. Тенденция становится очевидной, если использовать скользящее окно не менее 1 года.
plot_ts(ts, window=365)
Как видите, это четкий линейный нисходящий тренд. Это полезно при разработке моделей, поскольку большинство моделей требует, чтобы вы указали, существует ли компонент тренда и является ли он линейным (также называется «аддитивным») или нелинейным (также называется «мультипликативный»).
Обнаружение выбросов
Цель этого раздела - выявить выбросы и решить, как с ними бороться. На практике часто используются традиционные детерминированные методы, такие как построение графика распределения и обозначение выбросом всех наблюдений выше или ниже выбранного порогового значения. Например:
## Plot histogram ts.hist(color="black", bins=100)
## Boxplot ts.plot.box()
Этот метод работает, особенно когда вы хорошо знакомы со своими данными и уже знаете, какому процессу и распределению они следуют и, следовательно, какой порог работает лучше. Однако мне лично проще позволить алгоритму машинного обучения делать это за меня с любым набором данных временных рядов.
Давайте напишем функцию для автоматического обнаружения выбросов во временном ряду с использованием алгоритма кластеризации из библиотеки scikit-learn: машина векторов поддержки одного класса, она изучает границы распределения (называемые « support ») и, следовательно, может классифицировать любые точки, лежащие за пределами границы, как выбросы.
'''
Find outliers using sklearn unsupervised support vetcor machine.
:parameter
:param ts: pandas Series
:param perc: float - percentage of outliers to look for
:return
dtf with raw ts, outlier 1/0 (yes/no), numeric index
'''
def find_outliers(ts, perc=0.01, figsize=(15,5)):
## fit svm
scaler = preprocessing.StandardScaler()
ts_scaled = scaler.fit_transform(ts.values.reshape(-1,1))
model = svm.OneClassSVM(nu=perc, kernel="rbf", gamma=0.01)
model.fit(ts_scaled)
## dtf output
dtf_outliers = ts.to_frame(name="ts")
dtf_outliers["index"] = range(len(ts))
dtf_outliers["outlier"] = model.predict(ts_scaled)
dtf_outliers["outlier"] = dtf_outliers["outlier"].apply(lambda
x: 1 if x==-1 else 0)
## plot
fig, ax = plt.subplots(figsize=figsize)
ax.set(title="Outliers detection: found"
+str(sum(dtf_outliers["outlier"]==1)))
ax.plot(dtf_outliers["index"], dtf_outliers["ts"],
color="black")
ax.scatter(x=dtf_outliers[dtf_outliers["outlier"]==1]["index"],
y=dtf_outliers[dtf_outliers["outlier"]==1]['ts'],
color='red')
ax.grid(True)
plt.show()
return dtf_outliers
С помощью этой функции мы сможем выявлять выбросы, но что нам делать с ними после обнаружения? Здесь нет оптимальной стратегии: прогнозировать временные ряды проще без точек данных, которые значительно отличаются от других наблюдений, но удаление этих точек может глубоко изменить распределение данных. Если вы решили исключить выбросы, наиболее удобный способ их удаления - интерполяция.
Давайте напишем функцию для удаления выбросов после их обнаружения путем интерполяции значений до и после выброса.
'''
Interpolate outliers in a ts.
'''
def remove_outliers(ts, outliers_idx, figsize=(15,5)):
ts_clean = ts.copy()
ts_clean.loc[outliers_idx] = np.nan
ts_clean = ts_clean.interpolate(method="linear")
ax = ts.plot(figsize=figsize, color="red", alpha=0.5,
title="Remove outliers", label="original", legend=True)
ts_clean.plot(ax=ax, grid=True, color="black",
label="interpolated", legend=True)
plt.show()
return ts_clean
Теперь давайте воспользуемся этими функциями. Сначала мы обнаруживаем выбросы:
dtf_outliers = find_outliers(ts, perc=0.05)
Затем обработайте их:
## outliers index position outliers_index_pos = dtf_outliers[dtf_outliers["outlier"]==1].index ## exclude outliers ts_clean = remove_outliers(ts, outliers_idx=outliers_index_pos)
Для целей этого руководства я продолжу использовать необработанные временные ряды (включая выбросы), но удаление выбросов и построение модели на чистом временном ряду (без выбросов) также является хорошей стратегией.
Тест стационарности
Стационарный процесс - это случайный процесс, безусловное совместное распределение вероятностей которого не изменяется при сдвиге во времени. Следовательно, такие параметры, как среднее значение и дисперсия, также не меняются со временем , поэтому стационарные временные ряды легче прогнозировать.
Есть несколько способов установить, является ли временной ряд стационарным или нет, наиболее распространенными являются старая добрая визуализация, просмотр автокорреляции и выполнение статистических тестов.
Наиболее распространенным тестом является тест Дики-Фуллера (также называемый тест ADF), где нулевая гипотеза состоит в том, что временной ряд имеет единичный корень, другими словами, временной ряд не является стационарным . Мы проверим, можно ли отвергнуть нулевую гипотезу, сравнив значение p с выбранным порогом (α), чтобы, если значение p меньше, мы могли отклонить нулевую гипотезу и предположить, что временной ряд с уверенностью является стационарным. уровень 1-α (технически мы просто не можем сказать, что это не так):
Давайте напишем функцию, которая объединяет все эти методы и показывает фигуру, состоящую из следующего:
- Результат 95% (α = 0,05) теста ADF (который будет напечатан в заголовке выходного рисунка).
- Первый график будет отображать среднее значение и дисперсию первых x% данных, это графический тест: если свойства временного ряда постоянны, мы увидим, что 1-x% данных колеблется вокруг среднего. и в пределах диапазона дисперсии первых x% наблюдений
- Последние 2 графика отображают PACF и ACF.
'''
Test stationarity by:
- running Augmented Dickey-Fuller test wiht 95%
- plotting mean and variance of a sample from data
- plottig autocorrelation and partial autocorrelation
'''
def test_stationarity_acf_pacf(ts, sample=0.20, maxlag=30, figsize=
(15,10)):
with plt.style.context(style='bmh'):
## set figure
fig = plt.figure(figsize=figsize)
ts_ax = plt.subplot2grid(shape=(2,2), loc=(0,0), colspan=2)
pacf_ax = plt.subplot2grid(shape=(2,2), loc=(1,0))
acf_ax = plt.subplot2grid(shape=(2,2), loc=(1,1))
## plot ts with mean/std of a sample from the first x%
dtf_ts = ts.to_frame(name="ts")
sample_size = int(len(ts)*sample)
dtf_ts["mean"] = dtf_ts["ts"].head(sample_size).mean()
dtf_ts["lower"] = dtf_ts["ts"].head(sample_size).mean()
+ dtf_ts["ts"].head(sample_size).std()
dtf_ts["upper"] = dtf_ts["ts"].head(sample_size).mean()
- dtf_ts["ts"].head(sample_size).std()
dtf_ts["ts"].plot(ax=ts_ax, color="black", legend=False)
dtf_ts["mean"].plot(ax=ts_ax, legend=False, color="red",
linestyle="--", linewidth=0.7)
ts_ax.fill_between(x=dtf_ts.index, y1=dtf_ts['lower'],
y2=dtf_ts['upper'], color='lightskyblue', alpha=0.4)
dtf_ts["mean"].head(sample_size).plot(ax=ts_ax,
legend=False, color="red", linewidth=0.9)
ts_ax.fill_between(x=dtf_ts.head(sample_size).index,
y1=dtf_ts['lower'].head(sample_size),
y2=dtf_ts['upper'].head(sample_size),
color='lightskyblue')
## test stationarity (Augmented Dickey-Fuller)
adfuller_test = sm.tsa.stattools.adfuller(ts, maxlag=maxlag,
autolag="AIC")
adf, p, critical_value = adfuller_test[0], adfuller_test[1],
adfuller_test[4]["5%"]
p = round(p, 3)
conclusion = "Stationary" if p < 0.05 else "Non-Stationary"
ts_ax.set_title('Dickey-Fuller Test 95%: '+conclusion+
'(p value: '+str(p)+')')
## pacf (for AR) e acf (for MA)
smt.graphics.plot_pacf(ts, lags=maxlag, ax=pacf_ax,
title="Partial Autocorrelation (for AR component)")
smt.graphics.plot_acf(ts, lags=maxlag, ax=acf_ax,
title="Autocorrelation (for MA component)")
plt.tight_layout()
Давай запустим:
test_stationarity_acf_pacf(ts, sample=0.20, maxlag=30)
Результат теста Дики-Фуллера показывает, что временной ряд может быть стационарным, поскольку мы можем отклонить нулевую гипотезу нестационарности с уровнем достоверности 95% (p-значение 0,015 ‹α 0,05). Однако это не совпадает с «глазным тестом», поскольку мы видим, что временной ряд уходит от среднего значения после января 2015 года. Более того, мы не могли отклонить нулевую гипотезу нестационарности с уровнем достоверности 99%. (p-значение 0,015 ›α 0,01), и автокорреляция не может сходиться к нулю.
Мы проведем те же тесты после сравнения временных рядов. Различие может помочь стабилизировать среднее значение, удалив изменения в уровне наблюдений и, следовательно, исключив (или уменьшив) тенденцию и сезонность. В основном мы применим следующее преобразование:
diff [t] = y [t] - y [t-lag]
Теперь давайте попробуем разграничить временные ряды на 1 задержку и снова запустим предыдущую функцию.
diff_ts = ts - ts.shift(1) diff_ts = diff_ts[(pd.notnull(diff_ts))] test_stationarity_acf_pacf(diff_ts, sample=0.20, maxlag=30)
На этот раз мы можем отвергнуть нулевую гипотезу нестационарности с уровнем достоверности 95% и 99% (p-значение 0,000). Можно сделать вывод, что лучше предположить, что временной ряд не является стационарным.
Что касается графика автокорреляции, очевидно, что отрицательная сезонность наблюдается каждые 2 дня, что означает, что в начале недели меньше продаж, а положительная сезонность каждые 7 дней (больше продаж в выходные).
Анализ сезонности
Цель этого последнего раздела - понять, какая сезонность влияет на данные (еженедельная сезонность, если колебания происходят каждые 7 дней, сезонность месяца, если колебания происходят каждые 30 дней, и т. Д.).
Это очень важно для части проектирования модели, которая следует после этого сеанса анализа. В частности, при работе с моделями сезонной авторегрессии необходимо указать количество наблюдений за сезон: я говорю о параметре «s» в SARIMA (p, d, q) x ( P, D, Q, s).
В библиотеке statsmodel есть суперполезная функция, которая позволяет декомпозировать временные ряды. Эта функция разбивает данные на 3 компонента: тренд, сезонность и остатки.
Давайте построим разложение временного ряда с сезонностью 7 дней.
decomposition = smt.seasonal_decompose(ts, freq=7)
trend = decomposition.trend
seasonal = decomposition.seasonal
residual = decomposition.resid
fig, ax = plt.subplots(nrows=4, ncols=1, sharex=True, sharey=False)
ax[0].plot(ts)
ax[0].set_title('Original')
ax[0].grid(True)
ax[1].plot(trend)
ax[1].set_title('Trend')
ax[1].grid(True)
ax[2].plot(seasonal)
ax[2].set_title('Seasonality')
ax[2].grid(True)
ax[3].plot(residual)
ax[3].set_title('Residuals')
ax[3].grid(True)
Я обычно выбираю сезонный параметр, который приводит к меньшим остаткам. В этом случае при попытке использовать 2 дня, 7 дней и 30 дней результаты будут лучше при еженедельной сезонности (s = 7).
Заключение
Эта статья представляет собой руководство о том, как анализировать временные ряды реального мира с помощью статистики и машинного обучения, прежде чем приступить к построению модели прогнозирования. Результаты этого анализа полезны для разработки модели, которая может хорошо соответствовать временным рядам (что будет сделано в следующих руководствах, ссылки вверху). Особенно:
- мы можем включить компонент линейного тренда в наши модели прогнозирования
- мы можем обучить модель на необработанных данных, включая выбросы, и на обработанных данных без выбросов и проверить, какая из них работает лучше
- мы знаем, что временной ряд не является стационарным, и поэтому мы должны использовать модель AR-I-MA вместо ARMA.
- мы можем включить еженедельную сезонную составляющую в наши модели прогнозирования.
Надеюсь, вам понравилось! Не стесняйтесь обращаться ко мне с вопросами и отзывами или просто для того, чтобы поделиться своими интересными проектами.
Эта статья является частью серии Прогнозирование временных рядов с помощью Python, см. Также:
