Временные ряды и ARIMA: Как предсказывать будущее без хрустального шара

от автора

Часть 1

Что такое временной ряд, модель ARIMA и как к ней подбирать параметры.

Временной ряд — собранный в разные моменты времени статистический материал о значении каких-либо параметров (в простейшем случае одного) исследуемого процесса. (с) Википедия

Простым словами, временной ряд — это просто последовательность событий, которая как-то зависит от времени. Для начала будем считать, что ряд самый простецкий и нас просто есть скачущие туда-сюда точки, которые распределены по временной шкале.

Википедия предлагает такую страшную картинку

Википедия предлагает такую страшную картинку

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

Стационарный временной ряд на примере синуса

Стационарный временной ряд на примере синуса

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

Коротко о фичах временных рядов

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

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

1. Проверка на стационарность

from statsmodels.tsa.stattools import adfuller  def test_stationarity(series, title=''):     print(f"Results of ADF Test on {title}:")     result = adfuller(series, autolag='AIC')     print(f"ADF Statistic: {result[0]}")     print(f"p-value: {result[1]}")     if result[1] > 0.05:         print("Non-stationary")     else:         print("stationary")     for key, value in result[4].items():         print(f"Critical Value ({key}): {value}")     print("\n")  test_stationarity(df['count_1'], 'Counts')

Я юзала обычный тест Дики-Фуллера, который за меня красиво оформил chatGPT и библиотечная функция. Тут, как и везде в статистике, важно выбрать заранее p-value.
Про p-value и сам тест можно почитать дополнительно, а я хочу сосредоточиться на трёх моделях, которые использовала в работе.

2. ARIMA model

В любом руководстве по моделям ARIMA/SARIMAX вам первым делом скажут, что надо смотреть на автокорреляцию и частичную автокорреляцию, чтобы хоть как-то адекватно задать начальные параметры p и q. p и q смотрятся по выбросам: выбираются последние значения на горизонтальной оси, у которых точка выходит за границу синей области. Значение q смотрятся по графику ACF, a p — по графику PACF.
Смотреть их нужно на уже стационарном ряде. Но есть нюанс.

Идеальные ситуации из интернета

Идеальные ситуации из интернета
#Отрисовка ACF и PACF import statsmodels.api as sm  diff1 = dataset.diff() diff1.dropna(inplace=True) #delete all rows with NaN value  fig = plt.figure(figsize=(12,8)) ax1 = fig.add_subplot(211) half_size = len(diff1) // 2 fig = plot_acf(diff1, lags=half_size, ax=ax1) ax2 = fig.add_subplot(212) fig = sm.graphics.tsa.plot_pacf(diff1, lags=half_size, ax=ax2)

В моём случае это почти не сработало. Конечно, какие-то отправные p и q мне программа выдала, но для гиперпараметров модели они совсем не годились.

# лучше юзать библиотечные функции, но если даже делаете руками  # делайте длину через переменную, а не цифры length = int(len(diff1) * 0.8) train_data=diff1[0:length] test_data=diff1[length:]   model = sm.tsa.arima.ARIMA(diff1, order=(p, d, q)) model_fit = model.fit()  plt.figure(figsize=(10,6)) plt.plot(diff1, label='Исходные данные')  # Plot predicted data predicted_data = model_fit.predict() plt.plot(predicted_data, label='Предсказанные данные')  plt.legend() plt.show()  #выводит красивые картинки и прочую статистику по обучению print(model_fit.summary())

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

p = q = range(1, 40) d = range(1, 5) pdq = list(itertools.product(p, d, q))   def objective_arima(trial):     order=trial.suggest_categorical('order',pdq)     model=ARIMA(train_data,order=order)     mdl = model.fit() #disp=0     predictions = mdl.forecast(len(test_data))     predictions = pd.Series(predictions.values, index=test_data.index)     # predictions     residuals = test_data['count'] - predictions     mse=np.sqrt(np.mean(residuals**2))     accuracy=mse     return accuracy study=optuna.create_study(direction="minimize") study.optimize(objective_arima,n_trials=20)

Про что здесь нельзя забывать?

1. pd.Series(predictions.values, index=test_data.index) точно требует predictions.values, я с этим намучалась, потому что постоянны выскакивала ошибка из-за неверного типа данных
2. Обратное интегрирование: у меня долго предсказания висели внизу, а я ещё удивлялась, почему они так сильно не совпадают с тестовыми данными. Как вы уже догадались, дело было в константе интегрирования. Решается просто добавлением константы из хвоста тренировочных данных.
3. Надо аккуратнее оперировать индексами в DataFrame Pandas, большую часть error-ов я могла избежать, если бы вдумчивее подошла к вопросу.

Пример интегрирования и финального кода с метриками ниже:

inverted_forecast = pd.Series(predictions).cumsum() inverted_forecast += dataset[length - 1] inverted_true = dataset_sends[length:]  plt.figure(figsize=(10,6)) plt.plot(inverted_true, color = "green") plt.plot(inverted_forecast, color = "red") plt.legend(["Actual","prediction"])  plt.title("Predicted vs True Value") plt.xlabel("date") plt.ylabel("count_sends") plt.show()
from sklearn.metrics import mean_absolute_percentage_error  print("mean absolute error: ", round(mean_absolute_error(inverted_true['count'], inverted_forecast), 5)) print("mean squared error: ", round(mean_squared_error(inverted_true['count'], inverted_forecast), 5)) print("Root mean squared error: ", round(mean_squared_error(inverted_true['count'], inverted_forecast, squared=False), 5)) print("mean absolute percentage error", round(mean_absolute_percentage_error(inverted_true['count'], inverted_forecast)), "%")
Итог моих предсказаний

Итог моих предсказаний

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

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

Ермак Марина

Аналитик, SENSE IT


ссылка на оригинал статьи https://habr.com/ru/articles/821231/