В прошлом посте я рассказывала про свои мучения с моделькой 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