【ベンチマーク】複数の周期が重なるデータ

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

rng = np.random.default_rng(42)
n = 730  # 2年分
dates = pd.date_range("2022-01-01", periods=n, freq="D")
t = np.arange(n)

# 週次周期 + 年次周期 + 緩やかなトレンド
weekly = 10 * np.sin(2 * np.pi * t / 7)
yearly = 25 * np.sin(2 * np.pi * t / 365.25)
trend = 0.02 * t
noise = rng.normal(0, 5, n)
series = pd.Series(200 + trend + weekly + yearly + noise, index=dates)

h = 28
train_raw = series.iloc[:-h]
test_raw = series.iloc[-h:]

fig, axes = plt.subplots(2, 1, figsize=(12, 6), gridspec_kw={"height_ratios": [2, 1]})
axes[0].plot(series.index, series.values, color="#94a3b8", linewidth=0.5)
axes[0].axvline(series.index[-h], color="#ef4444", linestyle="--")
axes[0].set_title("週次 + 年次の複合季節データ（赤線より右がテスト）")
axes[0].grid(alpha=0.2)

# 直近2か月を拡大
recent = series.iloc[-60:]
axes[1].plot(recent.index, recent.values, color="#2563eb", linewidth=1)
axes[1].set_title("直近2か月の拡大（週次の波が見える）")
axes[1].grid(alpha=0.2)

fig.tight_layout()
plt.show()

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80

import warnings
from statsmodels.tsa.holtwinters import ExponentialSmoothing, SimpleExpSmoothing
from statsmodels.tsa.statespace.sarimax import SARIMAX
from statsmodels.tsa.arima.model import ARIMA

warnings.filterwarnings("ignore")

def compute_mase(actual, predicted, train, m=7):
    e = np.abs(actual - predicted)
    scale = np.abs(train.values[m:] - train.values[:-m]).mean()
    return e.mean() / scale if scale > 0 else np.nan

def make_fourier(index, period, order):
    """sin/cos のフーリエ特徴量を生成。"""
    t = np.arange(len(index))
    cols = {}
    for k in range(1, order + 1):
        cols[f"sin_{period}_{k}"] = np.sin(2 * np.pi * k * t / period)
        cols[f"cos_{period}_{k}"] = np.cos(2 * np.pi * k * t / period)
    return pd.DataFrame(cols, index=index)

# フーリエ特徴量（年次周期をsin/cosで表現）
fourier_all = make_fourier(series.index, period=365.25, order=3)
fourier_train = fourier_all.iloc[:-h]
fourier_test = fourier_all.iloc[-h:]

results = {}

# 1. 季節ナイーブ
last7 = train_raw.values[-7:]
results["季節ナイーブ(7)"] = np.tile(last7, 4)[:h]

# 2. SES
results["SES"] = SimpleExpSmoothing(
    train_raw, initialization_method="estimated"
).fit().forecast(h).values

# 3. HW加法 (sp=7)
results["HW加法(7)"] = ExponentialSmoothing(
    train_raw, trend="add", seasonal="add", seasonal_periods=7,
    initialization_method="estimated",
).fit().forecast(h).values

# 4. ARIMA(2,1,1)
results["ARIMA"] = ARIMA(train_raw, order=(2, 1, 1)).fit().forecast(steps=h).values

# 5. SARIMAX(1,0,1)(1,1,1,7)
results["SARIMAX(7)"] = SARIMAX(
    train_raw, order=(1, 0, 1), seasonal_order=(1, 1, 1, 7),
    enforce_stationarity=False, enforce_invertibility=False,
).fit(disp=False).forecast(steps=h).values

# 6. SARIMAX + フーリエ（年次）
model_f = SARIMAX(
    train_raw, exog=fourier_train, order=(1, 0, 1), seasonal_order=(1, 1, 1, 7),
    enforce_stationarity=False, enforce_invertibility=False,
).fit(disp=False)
results["SARIMAX+フーリエ"] = model_f.forecast(steps=h, exog=fourier_test).values

# MASE 計算
mase_scores = {name: compute_mase(test_raw.values, pred, train_raw, m=7)
               for name, pred in results.items()}

# 棒グラフ
names = list(mase_scores.keys())
scores = [mase_scores[n] for n in names]
colors = ["#10b981" if s < 1 else "#ef4444" for s in scores]

fig, ax = plt.subplots(figsize=(10, 5))
bars = ax.barh(names, scores, color=colors, edgecolor="white")
ax.axvline(1.0, color="#1e40af", linestyle="--", linewidth=1, label="MASE = 1（季節ナイーブ水準）")
ax.set_xlabel("MASE（低いほど良い）")
ax.set_title("複数季節データ: モデル別 MASE 比較")
ax.legend()
for bar, s in zip(bars, scores):
    ax.text(bar.get_width() + 0.05, bar.get_y() + bar.get_height() / 2,
            f"{s:.2f}", va="center", fontsize=9)
ax.grid(alpha=0.2, axis="x")
fig.tight_layout()
plt.show()

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30

from statsmodels.graphics.tsaplots import plot_acf

targets = [
    ("SARIMAX+フーリエ", "#2563eb"),
    ("SARIMAX(7)", "#f97316"),
    ("SES", "#ef4444"),
]

fig, axes = plt.subplots(len(targets), 2, figsize=(13, 3.5 * len(targets)))
test_idx = test_raw.index

for i, (name, color) in enumerate(targets):
    pred = results[name]
    residuals = test_raw.values - pred
    mase = mase_scores[name]

    ax_l = axes[i, 0]
    ax_l.plot(test_idx, test_raw.values, "o-", color="#94a3b8", markersize=3, label="実測")
    ax_l.plot(test_idx, pred, "s--", color=color, markersize=3, label="予測")
    ax_l.set_title(f"{name}  (MASE={mase:.2f})")
    ax_l.legend(fontsize=8)
    ax_l.grid(alpha=0.2)

    ax_r = axes[i, 1]
    plot_acf(residuals, ax=ax_r, lags=14, color=color, title="")
    ax_r.set_title("残差ACF")

fig.suptitle("複数季節データにおける誤差パターン比較", fontsize=13, y=1.02)
fig.tight_layout()
plt.show()