【ベンチマーク】構造変化のあるデータ

 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

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

rng = np.random.default_rng(42)
n = 48
dates = pd.date_range("2020-01-01", periods=n, freq="MS")
t = np.arange(n)

# 前半: 緩やかなトレンド, 後半: レベルシフト + 急トレンド
trend = np.where(t < 24, 0.5 * t, 0.5 * 24 + 30 + 1.0 * (t - 24))
noise = rng.normal(0, 4, n)
series = pd.Series(80 + trend + noise, index=dates)

h = 12
train_raw = series.iloc[:-h]
test_raw = series.iloc[-h:]
break_point = dates[24]

fig, ax = plt.subplots(figsize=(10, 5))
ax.plot(series.index, series.values, "o-", color="#94a3b8", markersize=4)
ax.axvline(break_point, color="#8b5cf6", linestyle="-", linewidth=1.5, label="構造変化点")
ax.axvline(series.index[-h], color="#ef4444", linestyle="--", linewidth=1, label="テスト開始")
ax.fill_betweenx([series.min() - 5, series.max() + 5],
                  series.index[0], break_point, alpha=0.08, color="#8b5cf6")
ax.set_title("構造変化を含む月次データ（紫=変化前, 赤破線=テスト開始）")
ax.set_ylabel("値")
ax.legend()
ax.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

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

warnings.filterwarnings("ignore")

def compute_mase(actual, predicted, train, m=1):
    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 forecast_model(train, h, name):
    if name == "季節ナイーブ":
        return np.full(h, train.values[-1])
    if name == "SES":
        return SimpleExpSmoothing(train, initialization_method="estimated").fit().forecast(h).values
    if name == "HW加法":
        return ExponentialSmoothing(
            train, trend="add", seasonal=None,
            initialization_method="estimated",
        ).fit().forecast(h).values
    if name == "ARIMA(1,1,1)":
        return ARIMA(train, order=(1, 1, 1)).fit().forecast(steps=h).values
    if name == "ARIMA(2,1,1)":
        return ARIMA(train, order=(2, 1, 1)).fit().forecast(steps=h).values
    raise ValueError(name)

models = ["季節ナイーブ", "SES", "HW加法", "ARIMA(1,1,1)", "ARIMA(2,1,1)"]

# 全期間で学習
train_all = series.iloc[:-h]
# 変化後のみで学習 (月25以降)
train_post = series.iloc[24:-h]

results_all = {}
results_post = {}
for m in models:
    try:
        results_all[m] = forecast_model(train_all, h, m)
    except Exception:
        results_all[m] = np.full(h, np.nan)
    try:
        results_post[m] = forecast_model(train_post, h, m)
    except Exception:
        results_post[m] = np.full(h, np.nan)

mase_all = {m: compute_mase(test_raw.values, results_all[m], train_all) for m in models}
mase_post = {m: compute_mase(test_raw.values, results_post[m], train_post) for m in models}

# グループ棒グラフ
x = np.arange(len(models))
width = 0.35

fig, ax = plt.subplots(figsize=(10, 5))
bars1 = ax.bar(x - width / 2, [mase_all[m] for m in models], width,
               label="全期間で学習", color="#ef4444", alpha=0.8)
bars2 = ax.bar(x + width / 2, [mase_post[m] for m in models], width,
               label="変化後のみで学習", color="#2563eb", alpha=0.8)
ax.axhline(1.0, color="#94a3b8", linestyle="--", linewidth=0.8)
ax.set_xticks(x)
ax.set_xticklabels(models, fontsize=9)
ax.set_ylabel("MASE")
ax.set_title("学習期間の選択がMASEに与える影響")
ax.legend()

for bar in bars1:
    h_val = bar.get_height()
    if not np.isnan(h_val):
        ax.text(bar.get_x() + bar.get_width() / 2, h_val + 0.05,
                f"{h_val:.2f}", ha="center", va="bottom", fontsize=8)
for bar in bars2:
    h_val = bar.get_height()
    if not np.isnan(h_val):
        ax.text(bar.get_x() + bar.get_width() / 2, h_val + 0.05,
                f"{h_val:.2f}", ha="center", va="bottom", fontsize=8)

ax.grid(alpha=0.2, axis="y")
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

fig, axes = plt.subplots(2, 2, figsize=(13, 8))
test_idx = test_raw.index

# 全期間 ARIMA vs 変化後 ARIMA
for i, (label, pred, color) in enumerate([
    ("全期間 ARIMA(1,1,1)", results_all["ARIMA(1,1,1)"], "#ef4444"),
    ("変化後 ARIMA(1,1,1)", results_post["ARIMA(1,1,1)"], "#2563eb"),
]):
    ax = axes[0, i]
    ax.plot(test_idx, test_raw.values, "o-", color="#94a3b8", markersize=4, label="実測")
    ax.plot(test_idx, pred, "s--", color=color, markersize=4, label="予測")
    m_val = compute_mase(test_raw.values, pred, train_all if "全期間" in label else train_post)
    ax.set_title(f"{label}  (MASE={m_val:.2f})")
    ax.legend(fontsize=8)
    ax.grid(alpha=0.2)

# 残差の分布
for i, (label, pred, train_used, color) in enumerate([
    ("全期間", results_all["ARIMA(1,1,1)"], train_all, "#ef4444"),
    ("変化後", results_post["ARIMA(1,1,1)"], train_post, "#2563eb"),
]):
    ax = axes[1, i]
    resid = test_raw.values - pred
    ax.bar(range(len(resid)), resid, color=color, alpha=0.7)
    ax.axhline(0, color="black", linewidth=0.5)
    ax.axhline(resid.mean(), color=color, linestyle="--", label=f"MBE={resid.mean():.1f}")
    ax.set_title(f"{label}学習の残差")
    ax.set_xlabel("テスト月")
    ax.legend(fontsize=8)
    ax.grid(alpha=0.2)

fig.suptitle("ARIMA(1,1,1) の学習期間別 予測 vs 残差", fontsize=13)
fig.tight_layout()
plt.show()