statsmodels надає багато інструментів для статистичного моделювання, включаючи роботу з часовими рядами. Я писав про інші бібліотеки, тому подумав, що варто також згадати про найпопулярнішу бібліотеку, яка є на ринку. Якщо ви використовуєте ARIMA або SARIMA, це буде логічним початком.
statsmodels охоплює моделювання як одновимірних, так і багатовимірних часових рядів. Вона містить безліч статистичних тестів для оцінки припущень моделей та їхніх результатів. Виведення результатів нагадує роботу з "MS Excel", де ви бачите ключові показники. І вона чудово працює з pandas та numpy.
Давайте розглянемо приклад. Ми проаналізуємо згенерований набір даних часових рядів, щоб продемонструвати деякі ключові можливості statsmodels.
Імпортуємо необхідні бібліотеки:
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.tsa.api import SimpleExpSmoothing, Holt, ExponentialSmoothing
from statsmodels.tsa.stattools import adfuller, acf, pacf
from statsmodels.tsa.arima.model import ARIMA
from statsmodels.tsa.seasonal import seasonal_decompose
Генерація або завантаження даних часового ряду
Згенеруємо часову серію з трендом і сезонністю:
"""
Генерація або завантаження даних часового ряду
Згенеруємо часову серію з трендом і сезонністю
"""
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
# Генерація згенерованих даних часового ряду
np.random.seed(42)
n = 200
time = pd.date_range(start="2023-01-01", periods=n, freq="D")
trend = np.linspace(10, 50, n)
seasonality = 10 * np.sin(np.linspace(0, 2 * np.pi, n))
noise = np.random.normal(0, 2, n)
data = trend + seasonality + noise
# Створення DataFrame; розподіл даних
df = pd.DataFrame({"date": time, "value": data})
df.set_index("date", inplace=True)
hold_out_days = 30
train = df.iloc[:-hold_out_days]
hold_out = df.iloc[-hold_out_days:]
# Побудова графіку
plt.figure(figsize=(10, 6))
plt.plot(df.index, df["value"], label="Повний набір даних", color="Blue")
plt.plot(hold_out.index, hold_out["value"], label="Тестові значення (Hold-Out)", color="Green")
plt.title("Згенерований часовой ряд з трендом і тестовими наборами")
plt.xlabel("Дата")
plt.ylabel("Значення")
plt.legend()
plt.grid()
plt.savefig("simulated_time_series.png")
plt.show()
Крок 3: Декомпозиція часового ряду
Використаємо seasonal_decompose, щоб розділити серію на тренд, сезонні та залишкові компоненти.
"""
Декомпозиція часового ряду
Використовуємо seasonal_decompose для розділення серії на тренд, сезонні та залишкові компоненти.
"""
from statsmodels.tsa.seasonal import seasonal_decompose
# Декомпозиція часового ряду
decomposition = seasonal_decompose(df["value"], model="additive", period=30)
# Побудова компонент
fig = decomposition.plot()
fig.set_size_inches(10, 8) # Налаштування розміру графіка
plt.suptitle("Декомпозиція часового ряду", fontsize=16, y=0.95) # Налаштування позиції заголовка
plt.tight_layout(rect=[0, 0, 1, 0.96]) # Забезпечуємо, щоб заголовок не накладався на підплоти
plt.savefig("time_series_decomposition.png")
plt.show()
Крок 4: Перевірка на стаціонарність
Використаємо тест на розширення Дікі-Фуллера (ADF), щоб оцінити стаціонарність.
"""
Перевірка на стаціонарність
Використовуємо тест на розширення Дікі-Фуллера (ADF) для оцінки стаціонарності.
"""
from statsmodels.tsa.stattools import adfuller
# Виконуємо тест ADF
result = adfuller(df["value"])
print(f"Статистика ADF: {result[0]:.4f}")
print(f"P-значення: {result[1]:.4f}")
if result[1] > 0.05:
print("Часовий ряд є нестатіонарним.")
else:
print("Часовий ряд є статіонарним.")
Статистика ADF: -0.5022
P-значення: 0.8916
Часовий ряд є нестатіонарним.
Автокореляція та часткова автокореляція
Візуалізуємо ACF та PACF, щоб визначити залежності затримок.
"""
Автокореляція та часткова автокореляція
Візуалізуємо ACF та PACF для визначення залежностей затримок.
"""
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
# Побудова ACF та PACF
fig, axes = plt.subplots(1, 2, figsize=(12, 6))
plot_acf(df["value"], lags=30, ax=axes[0])
plot_pacf(df["value"], lags=30, ax=axes[1])
plt.suptitle("Графіки ACF та PACF", fontsize=16)
plt.savefig("acf_pacf_plots.png")
plt.show()
Підганяємо модель ARIMA
Підганяємо модель ARIMA до даних для прогнозування.
"""
Підганяємо модель ARIMA
Підганяємо модель ARIMA до даних для прогнозування.
"""
from statsmodels.tsa.arima.model import ARIMA
# Підганяємо модель ARIMA(2,1,2)
model = ARIMA(df["value"], order=(2, 1, 2))
arima_result = model.fit()
print(arima_result.summary())
# Побудова діагностики залишків
arima_result.plot_diagnostics(figsize=(10, 6))
plt.savefig("arima_residuals_diagnostics.png")
plt.show()
Прогнозування майбутніх значень
Використовуйте підганяну модель для прогнозування залишкових значень.
"""
ARIMA
Прогнозуємо на 30 днів, які були відкладені
"""
# Підганяємо модель ARIMA на навчальних даних
model = ARIMA(train["value"], order=(2, 1, 2), freq="D") # Явно задаємо freq="D"
arima_result = model.fit()
# Прогнозуємо майбутні значення для періоду залишків
forecast = arima_result.get_forecast(steps=hold_out_days)
forecast_index = hold_out.index # Використовуємо той самий індекс, що й у періоді залишків
forecast_mean = forecast.predicted_mean
forecast_ci = forecast.conf_int()
# Розрахунок MAPE на періоді залишків
mape = mean_absolute_percentage_error(hold_out["value"], forecast_mean)
print(f"Середня абсолютна помилка у процентах (MAPE): {mape:.3%}")
# Побудова результатів
plt.figure(figsize=(10, 6))
plt.plot(train.index, train["value"], label="Навчальні дані", color="Blue")
plt.plot(hold_out.index, hold_out["value"], label="Відкладені (реальні значення)", color="Green")
plt.plot(forecast_index, forecast_mean, label="Прогноз", color="Red")
plt.fill_between(forecast_index, forecast_ci.iloc[:, 0], forecast_ci.iloc[:, 1], color="Red", alpha=0.2, label="Довірчий інтервал")
plt.title(f"Прогноз ARIMA (MAPE: {mape:.3%})")
plt.xlabel("Дата")
plt.ylabel("Значення")
plt.legend()
plt.grid()
plt.savefig("arima_forecast_holdout.png")
plt.show()
Наша модель прогнозує набагато повільніше зростання, ніж насправді має серія.
Експоненціальне згладжування Холта-Вінтерса
"""
Експоненціальне згладжування Холта-Вінтерса
"""
# Застосовуємо експоненціальне згладжування Холта-Вінтерса
hw_model = ExponentialSmoothing(
train["value"],
seasonal="add",
seasonal_periods=30
).fit()
hw_forecast = hw_model.forecast(steps=hold_out_days)
# Розрахунок MAPE на періоді залишків
mape_hw = mean_absolute_percentage_error(hold_out["value"], hw_forecast)
print(f"Holt-Winters MAPE: {mape_hw:.3%}")
# Побудова результатів
plt.figure(figsize=(10, 6))
plt.plot(train.index, train["value"], label="Навчальні дані", color="Blue")
plt.plot(hold_out.index, hold_out["value"], label="Відкладені (реальні значення)", color="Green")
plt.plot(hold_out.index, hw_forecast, label="Прогноз Холта-Вінтерса", color="Red")
plt.title(f"Прогноз Холта-Вінтерса \n MAPE: {mape_hw:.3%}")
plt.xlabel("Дата")
plt.ylabel("Значення")
plt.legend()
plt.grid()
plt.savefig("holt_winters_forecast.png")
plt.show()
Прогноз Холта-Вінтерса також не передбачає зміни дуже точно.
Тому я вирішив перевірити це за допомогою Tensorflow (Keras) з використанням моделі LSTM.
import tensorflow as tf
from tensorflow.keras.layers import LSTM
from sklearn.metrics import mean_absolute_percentage_error
from sklearn.preprocessing import MinMaxScaler
# Підготовка даних для LSTM
scaler = MinMaxScaler()
df["value"] = scaler.fit_transform(df["value"].values.reshape(-1, 1))
def create_lagged_features(data, lag):
X, y = [], []
for i in range(len(data) - lag):
X.append(data[i:i+lag])
y.append(data[i+lag])
return np.array(X), np.array(y)
lag = 10 # Кількість минулих спостережень для прогнозування
X, y = create_lagged_features(df["value"].values, lag)
X = X.reshape(X.shape[0], X.shape[1], 1)
# Розділяємо на навчальні та тестові набори
train_size = int(0.85 * len(X))
X_train, X_test = X[:train_size], X[train_size:]
y_train, y_test = y[:train_size], y[train_size:]
# Створюємо, підганяємо, прогнозуємо та оцінюємо модель LSTM
model = tf.keras.Sequential([
LSTM(50, activation='relu', input_shape=(lag, 1)),
tf.keras.layers.Dense(1)
])
model.compile(optimizer='adam', loss='mse')
model.summary()
model.fit(X_train, y_train, epochs=50, batch_size=8, verbose=1, validation_split=0.1)
y_pred_lstm = model.predict(X_test)
y_pred_lstm_inverse = scaler.inverse_transform(y_pred_lstm) # Обернене масштабування для порівняння
y_test_inverse = scaler.inverse_transform(y_test.reshape(-1, 1))
# Відновлення прогнозів для навчання для побудови графіку
train_predictions = model.predict(X_train)
train_predictions_inverse = scaler.inverse_transform(train_predictions)
# Розрахунок MAPE для тестового набору
mape = mean_absolute_percentage_error(y_test_inverse, y_pred_lstm_inverse)
print(f"LSTM MAPE: {mape:.3%}")
# Побудова результатів
plt.figure(figsize=(12, 8))
plt.plot(df.index, scaler.inverse_transform(df["value"].values.reshape(-1, 1)), label="Реальні дані", color="Blue")
train_index = df.index[lag:train_size + lag]
plt.plot(train_index, train_predictions_inverse, label="Прогнози для навчання", color="Orange")
test_index = df.index[train_size + lag:]
plt.plot(test_index, y_test_inverse, label="Відкладені (реальні значення)", color="Green")
plt.plot(test_index, y_pred_lstm_inverse, label="Прогнози для тестування", color="Red")
plt.title(f'Прогноз LSTM')
MAPE: {mape:.3%}')
plt.xlabel("Дата")
plt.ylabel("Значення")
plt.legend()
plt.grid()
plt.savefig("LSTM_forecast_with_holdout.png")
plt.show()
Ця модель дала набагато кращі результати, ніж моделі ARIMA та Holt-Winters. Хоча вона все ще недооцінює значення, тренд набагато ближчий до реальності, ніж у інших моделей.
Мене це не сильно турбує. У ARIMA є свої обмеження, але я все одно люблю statsmodels як бібліотеку — кращі прогнози з Tensorflow є результатом іншої моделі, а не кращого API.
Основні висновки
statsmodels — це орієнтир, який я використовую для порівняння з іншими бібліотеками для роботи з часовими рядами. Вона не є ідеальною, але досить гарна, і якщо б я міг використовувати лише одну бібліотеку на все життя, я б вибрав statsmodels.
Перекладено з: Time Series Analysis with statsmodels in Python