Часть 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