【はじめてのSARIMAX】SARIMAXを使ってこのブログの記事クリック数を予測してみた

plt.figure(figsize = (12,6))
plt.plot(df["クリック数"])
plt.show()

from statsmodels.graphics.tsaplots import plot_acf, plot_pacf

# 自己相関と偏自己相関のプロット
plt.figure(figsize=(12, 6))
plot_acf(df["クリック数"], lags=30, ax=plt.subplot(121))  # 自己相関
plot_pacf(df["クリック数"], lags=30, ax=plt.subplot(122))  # 偏自己相関
plt.tight_layout()
plt.show()

# 必要なライブラリをインポート
import matplotlib.pyplot as plt

# 曜日の情報を追加する
df["曜日"] = df.index.dayofweek

# 曜日ごとのデータを集計
weekly_data = df.groupby("曜日")["クリック数"].apply(list)

# 箱ひげ図をプロット
plt.figure(figsize=(10, 6))
plt.boxplot(weekly_data, labels=["Mon", "Tue", "Wed", "Thu", "Fri", "Sat", "Sun"])
plt.title("Distribution of data by day of the week")
plt.xlabel("Week")
plt.ylabel("the number of clicks")
plt.grid(axis="y")
plt.show()

import pandas as pd
import matplotlib.pyplot as plt

# 大学休み期間の設定
df["大学休み"] = 0
df.loc["2024-02-01":"2024-03-31", "大学休み"] = 1  # 春休み
df.loc["2024-08-01":"2024-09-30", "大学休み"] = 1  # 夏休み
df.loc["2023-12-25":"2024-01-05", "大学休み"] = 1  # 冬休み
df.loc["2024-12-25":"2024-12-29", "大学休み"] = 1  # 冬休み

# 比較用に月別データを抽出
df["month"] = df.index.month
comparison_data = df[(df.index >= "2023-12-20") & (df.index <= "2024-01-10")]  # 冬休みの期間とその前後を比べる
comparison_data["period"] = "Non Vacation"
comparison_data.loc["2023-12-25":"2024-01-05", "period"] = "Vacation"

# 箱ひげ図の作成
plt.figure(figsize=(8, 6))
comparison_data.boxplot(column="クリック数", by="period", grid=False)
plt.title("Comparison of Clicks: Vacation vs Non-Vacation (Winter Vacation)")
plt.suptitle("")  # 不要なタイトルを削除
plt.xlabel("Period")
plt.ylabel("Clicks")
plt.show()

# 曜日の情報を追加する
df["曜日"] = df.index.dayofweek
# 金曜日または土曜日を判定するダミー変数を作成
df["金土"] = df["曜日"].isin([4, 5]).astype(int)

# 大学の休み期間の設定
df["大学休み"] = 0
df.loc["2024-02-01":"2024-03-31", "大学休み"] = 1
df.loc["2024-08-01":"2024-09-30", "大学休み"] = 1
df.loc["2023-12-25":"2024-01-05", "大学休み"] = 1
df.loc["2024-12-25":"2024-12-29", "大学休み"] = 1

# Train, Test の分割
train_data = df.loc[:"2024-10-31"].copy()
test_data_11 = df.loc["2024-11-01":"2024-11-30"].copy()
test_data_12 = df.loc["2024-12-01":"2024-12-29"].copy()

# 考えられるパラメータの組み合わせを全て作成
max_p = 2
max_q = 2
max_d = 2
max_sp = 1
max_sd = 1
max_sq = 1
params = []

for p in range(0, max_p + 1):
    for q in range(0, max_q + 1):
        for d in range(0, max_d + 1):
            for sp in range(0, max_sp + 1):
                for sd in range(0, max_sd + 1): 
                    for sq in range(0, max_sq + 1): 
                        params.append([p,d,q,sp,sd,sq])

import warnings
from statsmodels.tsa.statespace.sarimax import SARIMAX

# ベストパラメータの初期化
best_param = None
best_param_seasonal = None
best_aic = float("inf")

# 警告をキャッチするためのリスト
excluded_params = []

# パラメータループ
for param in params:
    try:
        with warnings.catch_warnings(record=True) as w:
            # 警告を一時的にキャッチ
            warnings.simplefilter("always")
            
            # SARIMAモデルを構築
            model = SARIMAX(
                train_data["クリック数"],
                order=(param[0], param[1], param[2]),
                seasonal_order=(param[3], param[4], param[5], 7),
                exog=train_data[["大学休み", "金土"]]
            )
            results = model.fit(disp=False)

            # 警告が発生したか確認
            if any("Non-stationary" in str(warn.message) or "Non-invertible" in str(warn.message) or "failed to converge" in str(warn.message) for warn in w):
                print(f"Excluded ARIMA{(param[0], param[1], param[2])}x{(param[3], param[4], param[5], 7)} due to warnings.")
                excluded_params.append(param)
                continue  # このパラメータをスキップ

            # AICを比較してベストパラメータを更新
            print("SARIMAX{}x{}7 - AIC:{}".format((param[0], param[1], param[2]), (param[3], param[4], param[5], 7), results.aic))
            if best_aic > results.aic:
                best_param = (param[0], param[1], param[2])
                best_param_seasonal = (param[3], param[4], param[5], 7)
                best_aic = results.aic

    except Exception as e:
        print(f"Failed for SARIMAX{(param[0], param[1], param[2])}x{(param[3], param[4], param[5], 7)}: {e}")
        continue  # エラーが発生した場合も次のパラメータに進む

# ベストパラメータの結果を表示
print(f"BEST SARIMAX{best_param}x{best_param_seasonal}7 - AIC:{best_aic}")
if excluded_params:
    print(f"Excluded parameters: {excluded_params}")

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from statsmodels.tsa.statespace.sarimax import SARIMAX
from sklearn.metrics import mean_absolute_percentage_error
# 11月
test_data = test_data_11

# 予測結果と実際の値を格納するリスト
predictions = []
actuals = []

# 学習データをコピーしてローリングに使用
history = train_data.copy()

# テストデータを1日ずつ予測
for date in test_data.index:
    # SARIMAモデルの再構築・再学習
    model = SARIMAX(
        history["クリック数"],
        order=(1, 1, 1),
        seasonal_order=(0, 1, 1, 7),
        exog=history[["大学休み", "金土"]]
    )
    model_fit = model.fit(disp=False)

    # テストデータの1日分を予測
    exog_test = test_data.loc[[date], ["大学休み", "金土"]]
    forecast = model_fit.get_forecast(steps=1, exog=exog_test)
    prediction = forecast.predicted_mean.iloc[0]

    # 予測値と実際値を記録
    predictions.append(prediction)
    actuals.append(test_data.loc[date, "クリック数"])

    # 実測値を履歴（学習データ）に追加
    new_row = test_data.loc[[date], ["クリック数", "大学休み", "金土"]]
    history = pd.concat([history, new_row])

# 結果をデータフレームにまとめる
results = pd.DataFrame({
    "Actual": actuals,
    "Prediction": predictions
}, index=test_data.index)

# 結果のプロット
plt.figure(figsize=(20, 6))
plt.plot(results.index, results["Actual"], label="Actual", color="blue")
plt.plot(results.index, results["Prediction"], label="Prediction", color="red")
plt.title("Rolling Forecast with SARIMAX")
plt.legend()
plt.xlabel("Date")
plt.ylabel("Click Counts")
plt.grid()
plt.show()

# MAPEの計算
mape_first = mean_absolute_percentage_error(results["2024-11-01":"2024-11-15"]["Actual"], results["2024-11-01":"2024-11-15"]["Prediction"])
mape_last = mean_absolute_percentage_error(results["2024-11-16":"2024-11-30"]["Actual"], results["2024-11-16":"2024-11-30"]["Prediction"])
mape_all = mean_absolute_percentage_error(results["Actual"], results["Prediction"])
print(f"MAPE(11月1日～11月15日): {mape_first:.2%}")
print(f"MAPE(11月16日～11月30日): {mape_last:.2%}")
print(f"MAPE(11月全体): {mape_all:.2%}")