SARIMAX vs Экспоненциальное сглаживание: Когда простота побеждает

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

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

Данные

У меня была статистика компании за 2.5 года по дням, для пары экспериментов я пробовала склеить её в недельные данные. К сожалению или к счастью никакой физики за процессом тут не лежит. Требовалась просто сносная модель, которая бы предсказывала данные на следующий месяц и при этом давала результат лучше, чем предсказания аналитиков через конверсии. И была возможность вытащить из данных ещё дополнительную колонку внешних параметров.

SARIMA

На модель SARIMA у меня были большие надежды, всё-таки в отличие от ARIMA она могла оперировать внешними параметрами и с сезонностью у неё дела обстояли сильно лучше. В итоге я остановилась на статистике по дням.

Чтобы работать со стационарным рядом, я дифференцировала изначальный:

 diff1 = df.copy(deep=True)  diff1['1'] = df['1'].diff() diff1['2'] = df['2'].diff() diff1.dropna(inplace=True)

Затем сделала пробный тычок пальцем в небо, чтобы посмотреть на работу модели:

from statsmodels.tsa.statespace.sarimax import SARIMAX # Ensure the SARIMAX model is initialized properly with your training data and exogenous variables model = SARIMAX(diff1['1'][:962], exog=diff1['2'][:962], order=(1, 1, 1), seasonal_order=(1, 0, 1, 12))  # Fit the model result = model.fit(disp=False) # Forecast using the correct steps and aligned exogenous data forecast = result.get_forecast(steps=len(df_test), exog=forecast_exog) predictions = result.forecast(len(df_test), exog=forecast_exog) # predictions predictions = pd.Series(predictions.values, index=df_test.index) # If the forecast returns NaNs, check alignment and data quality forecast_df = pd.DataFrame({     'predicted_mean': forecast.predicted_mean,     'lower_ci': forecast.conf_int().iloc[:, 0],     'upper_ci': forecast.conf_int().iloc[:, 1] }, index=df_test.index)

Оптимизация с помощью Optuna:

 #Generate all different combinations of p, d and q triplets #Generate all different combinations of p, d, q and s triplets import itertools import optuna  p = range(1, 3)  # Smaller range for non-seasonal AR component d = range(0, 2) q = range(0, 3) P = range(0, 2)  # Smaller range for seasonal AR component D = range(0, 2) Q = range(0, 2) s = 12  pdq = list(itertools.product(p, d, q)) pdqs = [(x[0], x[1], x[2], s) for x in list(itertools.product(P, D, Q))] #%% def objective_sarima(trial):     non_seasonal_order = trial.suggest_categorical('non_seasonal_order', pdq)     seasonal_order = trial.suggest_categorical('seasonal_order', pdqs)     trend = trial.suggest_categorical('trend', ['n', 'c', 't', 'ct', None])          # # Generate predictions     # predictions = mdl.forecast(len(df_test))     model = SARIMAX(diff1['1'][:962], exog=diff1['2'][:962], order=non_seasonal_order, seasonal_order=seasonal_order)     # Fit the model     result = model.fit(disp=True)     # Forecast using the correct steps and aligned exogenous data     forecast = result.get_forecast(steps=len(df_test), exog=forecast_exog)     predictions = result.forecast(len(df_test), exog=forecast_exog)     # predictions     predictions = pd.Series(predictions.values, index=df_test.index)     # If the forecast returns NaNs, check alignment and data quality     forecast_df = pd.DataFrame({         'predicted_mean': forecast.predicted_mean,         'lower_ci': forecast.conf_int().iloc[:, 0],         'upper_ci': forecast.conf_int().iloc[:, 1]     }, index=df_test.index)     # Calculate residuals and error metric     residuals = diff1['1'] - predictions     mse = np.sqrt(np.mean(residuals**2))     return mse  study=optuna.create_study(direction="minimize") study.optimize(objective_sarima,n_trials=10)

Затем я на всё это счастье визуализировала и смотрела. Получилось устрашающе:

plt.figure(figsize=(10, 6)) plt.plot(df_test['1'], label='Actual Future', marker='o', color='green') plt.plot(forecast_df['predicted_mean'], label='Forecasted', marker='o', color='red') plt.fill_between(forecast_df.index, forecast_df['lower_ci'], forecast_df['upper_ci'], color='red', alpha=0.2) plt.title('Forecasted vs Actual Counts') plt.xlabel('Date') plt.ylabel('Count') plt.legend() plt.grid(True) plt.xticks(rotation=45) plt.tight_layout() plt.show()

После я сравнила ещё и обратно проинтегрированный ряд с исходными данными, которые мне предоставили аналитики. Меня просили оформить всё с табличкой-сравнением, так что код получился монструозный. Тут важно опять-таки не забывать про константу интегрирования, которая в коде носит имя offset. Иначе это чеховское ружьё выстрелит в ногу.

#integral constant! offset = int(list(dataset_sends['count_offers'])[-len(df_test)]) print(offset) inverted_forecast = pd.Series(forecast_df['predicted_mean']).cumsum() inverted_forecast = pd.DataFrame(forecast_df['predicted_mean']) inverted_true = pd.DataFrame(list(dataset_sends['count_offers'][-len(df_test):]), index=range(0, len(df_test)))  inverted_forecast = inverted_forecast + offset inverted_forecast = pd.DataFrame(list(inverted_forecast['predicted_mean'][-len(df_test):]), index=range(0, len(df_test))) dates = pd.DataFrame(list(dataset_sends['date'][-len(df_test):]), index=range(0, len(df_test)))  df_merged = pd.concat([inverted_forecast, inverted_true, dates], axis=1)  df_merged.columns = ['predicted_mean', 'count_offers', 'date']  df_merged.dropna(inplace=True) df_merged['1'] = df_merged['1'].apply(int) df_merged['predicted_mean'] = df_merged['predicted_mean'].apply(int) # Absolute Error df_merged['error'] = abs(df_merged['count_offers'] - df_merged['predicted_mean']) df_merged['date'] = pd.to_datetime(df_merged['date']) # Set 'date' column as the index df_merged.set_index('date', inplace=True)  # Group by week and sum the other columns df_merged = df_merged.resample('W').sum()  # Reset index to make 'date' a column again df_merged.reset_index(inplace=True) # Relative Error df_merged['relative_error, %'] = abs(df_merged['count_offers'] - df_merged['predicted_mean']) / df_merged['count_offers'] * 100 df_merged.dropna(inplace=True)   display(df_merged[-17:-1]) plt.plot(df_merged.index, df_merged['count_offers'], label='count_offers') # Plot 'predicted_mean' plt.plot(df_merged.index, df_merged['predicted_mean'], label='predicted_mean')  # Set labels and title plt.xlabel('Time') plt.ylabel('Count') plt.title('Comparison of count_offers and predicted_mean') plt.legend()  # Show plot plt.show()

Увы, показать табличку я не могу из-за NDA, и так хожу по офигенно тонкому льду:)

Вывод: хотя я и питала большие надежды, результат SARIMA меня не удовлетворил. И я пошла гуглить дальше, что ещё умного делают люди для предсказания временных рядов.

Экспоненциальное сглаживание

Этот подход упоминают в Яндекс-учебнике по ML, но уделяют не так много внимания, как моделям ARIMA/SARIMA. А вот мне он дал неожиданно классные результаты.

Здесь я тоже юзала пакетную историю и работала с данным, лежащими в переменной aust. Также руками поставила период в 40 дней, который подбирала ручками, без всяких оптимизаторов. Не делайте так:)

import pandas as pd import numpy as np import matplotlib.pyplot as plt from datetime import datetime import statsmodels.api as sm from pandas.plotting import register_matplotlib_converters from statsmodels.tsa.holtwinters import SimpleExpSmoothing from statsmodels.tsa.holtwinters import ExponentialSmoothing from sklearn.metrics import mean_squared_error,mean_absolute_error import warnings from statsmodels.tsa.holtwinters import SimpleExpSmoothing, Holt warnings.filterwarnings('ignore')  period = 40 fit1 = ExponentialSmoothing(     aust,     seasonal_periods=period,     trend="add",     seasonal="add",     use_boxcox=True,     initialization_method="estimated", ).fit() fit2 = ExponentialSmoothing(     aust,     seasonal_periods=period,     trend="add",     seasonal="mul",     use_boxcox=True,     initialization_method="estimated", ).fit() fit3 = ExponentialSmoothing(     aust,     seasonal_periods=period,     trend="add",     seasonal="add",     damped_trend=True,     use_boxcox=True,     initialization_method="estimated", ).fit() fit4 = ExponentialSmoothing(     aust,     seasonal_periods=period,     trend="add",     seasonal="mul",     damped_trend=True,     use_boxcox=True,     initialization_method="estimated", ).fit() results = pd.DataFrame(     index=[r"$\alpha$", r"$\beta$", r"$\phi$", r"$\gamma$", r"$l_0$", "$b_0$", "SSE"] ) params = [     "smoothing_level",     "smoothing_trend",     "damping_trend",     "smoothing_seasonal",     "initial_level",     "initial_trend", ] results["Additive"] = [fit1.params[p] for p in params] + [fit1.sse] results["Multiplicative"] = [fit2.params[p] for p in params] + [fit2.sse] results["Additive Dam"] = [fit3.params[p] for p in params] + [fit3.sse] results["Multiplica Dam"] = [fit4.params[p] for p in params] + [fit4.sse]    ax = train_data.plot(marker="o",     color="black",     title="Forecasts from Holt-Winters' multiplicative method", )  ax.set_xlim(1697587200000000000, 1713139200000000000) ax.set_ylim(0, 20) ax.set_ylabel("counts") ax.set_xlabel("Year")  fit1.fittedvalues.plot(style="--",marker="o", color="red", ax=ax) fit2.fittedvalues.plot(style="--", marker="o",color="green", ax=ax) fit3.fittedvalues.plot(style="--", marker="o",color="blue", ax=ax) fit4.fittedvalues.plot(style="--", marker="o",color="green", ax=ax)   forecast1 = fit1.forecast(121).rename("Holt-Winters (add-add-seasonal)") forecast2 = fit2.forecast(121).rename("Holt-Winters (add-mul-seasonal)") forecast3 = fit3.forecast(121).rename("Holt-Winters (add-add-seasonal heuristic)") forecast4 = fit4.forecast(121).rename("Holt-Winters (add-mul-seasonal heuristic)")  predictions1 = pd.Series(forecast1.values, index=test_data.index) predictions2 = pd.Series(forecast2.values, index=test_data.index) predictions3 = pd.Series(forecast3.values, index=test_data.index) predictions4 = pd.Series(forecast4.values, index=test_data.index)  predictions1.dropna(inplace=True) predictions2.dropna(inplace=True) predictions3.dropna(inplace=True) predictions4.dropna(inplace=True)  print("1mean absolute error : ",round(mean_absolute_error(test_data[:len(predictions1)], predictions1), 5)) print("1mean squared error : ",round(mean_squared_error(test_data[:len(predictions1)], predictions1), 5)) print("1Root mean squared error : ",round(mean_squared_error(test_data[:len(predictions1)], predictions1,squared=False), 5))  print("2mean absolute error : ",round(mean_absolute_error(test_data[:len(predictions2)], predictions2),5)) print("2mean squared error : ",round(mean_squared_error(test_data[:len(predictions2)], predictions2),5)) print("2Root mean squared error : ",round(mean_squared_error(test_data[:len(predictions2)], predictions2,squared=False),5))  print("3mean absolute error : ",round(mean_absolute_error(test_data[:len(predictions3)], predictions3),5)) print("3mean squared error : ",round(mean_squared_error(test_data[:len(predictions3)], predictions3),5)) print("3Root mean squared error : ",round(mean_squared_error(test_data[:len(predictions3)], predictions3,squared=False),5))  # print("4mean absolute error : ",round(mean_absolute_error(test_data[:len(predictions4)], predictions4),5)) # print("4mean squared error : ",round(mean_squared_error(test_data[:len(predictions4)], predictions4),5)) # print("4Root mean squared error : ",round(mean_squared_error(test_data[:len(predictions4)], predictions4,squared=False),5))  forecast1.plot(ax=ax, style="--", marker="o", color="green", legend=True, figsize=(20, 10)) forecast2.plot(ax=ax, style="--", marker="o", color="red", legend=True) forecast3.plot(ax=ax, style="--", marker="o", color="blue", legend=True) forecast4.plot(ax=ax, style="--", marker="o", color="yellow", legend=True)  plt.plot(predictions1, marker="o", color = "green") plt.plot(predictions2, marker="o", color = "red") plt.plot(predictions3, marker="o", color = "blue") plt.plot(predictions4, marker="o", color = "yellow") plt.plot(test_data, marker="o", color = 'black') plt.plot(train_data, marker="o", color = 'black') # Show the plot plt.show()  print("Figure 7.6: Forecasting.") # results

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

К тому же я не смогла по-человечески вписать сюда историю с адекватным временем. Мне пришлось в самом начале перевести всё в цифирьки и довольствоваться малым:

aust['date'] = pd.to_numeric(pd.to_datetime(aust['date']))

Исходная картинка у меня пропала, но картинка была какая-то такая:

Затем я взяла среднее по всем эти предсказаниям и получила вполне сносную картинку

combined_predictions = pd.DataFrame({     'predictions1': predictions1,     'predictions2': predictions2,     'predictions3': predictions3,     'predictions4': predictions4 }) # Calculate the average prediction across all models average_prediction = combined_predictions.mean(axis=1)  from sklearn.metrics import mean_absolute_percentage_error plt.figure(figsize=(15, 5)) plt.plot(average_prediction, marker="o", color = "red") plt.plot(test_data, marker="o", color = 'black') plt.ylim(0, 20)  plt.show() print("mean absolute error : ",round(mean_absolute_error(test_data, average_prediction),5)) print("mean squared error : ",round(mean_squared_error(test_data, average_prediction),5)) print("Root mean squared error : ",round(mean_squared_error(test_data, average_prediction,squared=False),5)) print("mean relative error", round(mean_absolute_percentage_error(test_data, average_prediction)), "%")

Вывод:
1. Вероятно, такую хорошую картинку я получила именно из-за сглаживания. Про это писали среди комментов к прошлому посту.
2. Временную шкалу мне пришлось дополнительно обрабатывать, чтобы всё заработало. Можно ли было обойтись без этого — не знаю.

В заключительной части этой «временной» серии расскажу про мой опыт использования пакета FEDOT. И если успею потестить что-то из предложенного ранее в комментах — тоже напишу.

Ермак Марина

Аналитик, SENSE IT