業務でAIばかり考えることになっているのですが、データからモデルを学習し、そのモデルを使って未知のデータに対して推論を行う技術群をAIと呼ぶとすると、いわゆる統計モデルでの機械学習も当然その中には含まれます。従って、統計モデルを用いた異常検知も当然AIという広大なスペクトラム上にある。

異常検知の第一歩は「正常データを予測できるモデル」を作ることです。そこで本稿ではまず、自己回帰和分移動平均モデル（ARIMA）およびその季節拡張版（SARIMA）を使い、時系列データの予測を行ってみました。

題材

今回のデータセットはKaggleのとしました。 1981年から1990年の日次気温データであり、確かメルボルンのデータだったはずです。

概形

時系列データを予測に使う上では、まずはそのデータについてよく知らねばなりません。データは「日時（分単位）」「温度」の2カラムのみ構成されていますが、温度値に先頭?がつく欠損表現が混在していました。まずはこの前処理を行い、時系列プロットでデータの概観を確認します。

import pandas as pd

df = pd.read_csv(
    "data/daily-minimum-temperatures-in-me.csv",
    parse_dates=["Date"],
    index_col="Date",
).rename(columns={
    'Daily minimum temperatures in Melbourne, Australia, 1981-1990': 'Temperature'
}).asfreq('D')

# 「?」を空文字に置き換え → 数値変換
df['Temperature'] = (
    df['Temperature']
      .str.replace(r'\?', '', regex=True)      # 「?」を削除
      .str.strip()                             # 前後空白を削除
      .replace('', pd.NA)                      # 空文字は欠損に
      .astype(float)                           # float にキャスト
)

df.plot(title='Daily Minimum Temperatures', xlabel='Date', ylabel='Temperature')

結果がこちら。年次の周期性があり、夏冬で上下するという季節性も持った面白いデータです。

Daily Minimum Temperatures

データは3652個から構成されており、これをそのまま解析にかけるのは私のローカルPCだと難しそうです。難しいというか、実際にSARIMAモデルの学習を進めてみたんですが、メモリが60GBくらい持って行かれて SIGKILL されました。

print(df.shape) 
# (3652, 1)

ダウンサンプリング

原データは3,652件あります。モデル学習時のメモリ負荷を軽減するため、週次にダウンサンプリングして取り扱う件数を減らします。ここでは週の平均気温を当該週の気温としました。

weekly = df['Temperature'].resample('W').mean()
weekly.plot(title='Weekly Average Temperatures', xlabel='Date', ylabel='Temperature')

print(weekly.info())
# DatetimeIndex: 523 entries, 1981-01-04 to 1991-01-06
# Freq: W-SUN
# Series name: Temperature
# Non-Null Count  Dtype  
# --------------  -----  
# 523 non-null    float64
# dtypes: float64(1)
# memory usage: 24.3 KB

STL分解

時系列予測では、データを「長期トレンド」「季節成分」「残差」に分解するとモデル構築が容易になります。STL（Seasonal–Trend decomposition using Loess）は、これらを局所回帰（Loess）で推定し、各成分を明示的に抽出する手法です。データの長期傾向や周期性を明示的に抽出し、残差部分を異常検知や予測モデルの誤差解析に利用できるようになります。

周期成分が１年（52週）であることは容易に想像できるので、それを前提にしてSTL分解してみます。

from statsmodels.tsa.seasonal import STL
stl = STL(weekly, period=52, robust=True)
stlres = stl.fit()

_ = stlres.plot()

STL分解

グラフからの定性分析ですが、およそ1年(52週)ごとに「山（夏）→ 谷（冬）」のサイクルが繰り返されているわけですが、1981年初めから1984年頃までは緩やかに低下し、その後1985年前後で最も低くなり、1986年以降は再び上昇基調にあることがわかります。一方でSeason（季節成分）に目を向けると、当たり前ですが毎年同じ時期に高温・低温が現れる形がほぼ一定幅で繰り返されており、年ごとの山の高さ・谷の深さの差は比較的小さく、大きな年変動は見られません。残差（Resid）の部分は大半は−2 ～ +2 の範囲に収まっていますが、外れ値が散発的に発生していることがわかります。

ADF検定

時系列データ解析において、モデル構築の前提条件として「定常性」の確認は欠かせません。定常性とは、データの平均・分散・自己共分散が時間経過に依存せず一定である性質を指します。多くの統計モデル（ARMAやSARIMAなど）は、定常過程を前提としており、非定常なデータにそのまま適用すると誤った推定や予測精度の低下を招きます。

ADF検定（Augmented Dickey–Fuller test）は、時系列に「単位根」があるかを調べる検定です。単位根があると非定常過程と判断され、差分を取らずにモデルを適用すると誤った推定を招くため、検定で定常性を確認します。このため、ADF検定を実行し、以下の帰無仮説と対立仮説を検証します。

帰無仮説 $H_0$：時系列に単位根がある（非定常過程）
対立仮説 $H_1$：時系列に単位根がない（定常過程）

from statsmodels.tsa.stattools import adfuller

adf_result = adfuller(weekly)
print("ADF Statistic: %f" % adf_result[0])
print("p-value: %f" % adf_result[1])
for key, critical_value in adf_result[4].items():
    print(f"Critical Value ({key}): {critical_value:.3f}")
print("#Lags Used:", adf_result[2])

ADF検定の結果は以下のとおりになりました。

項目	値
ADF Statistic	-9.273936
p-value	0.000000
Critical Value (1%)	-3.443
Critical Value (5%)	-2.867
Critical Value (10%)	-2.570
使用ラグ数 (#Lags)	19

ADF検定における検定統計量 (–9.27) はすべての臨界値（–3.44, –2.87, –2.57）よりも小さく、P値もほぼ0であるため、帰無仮説（非定常過程である）は強く棄却されます。つまり、ダウンサンプリングした週次時系列は統計的に定常と判断できます。

ACF/PACF

時系列データの場合、現在の値がどのくらい過去の値と相関があるのかは重要な構造情報です。ここで活用するのがACF（Autocorrelation Function：自己相関関数）とPACF（Partial Autocorrelation Function：偏自己相関関数）のグラフです。

百聞は一見に如かずなので、まずはグラフを見てみましょう。

import matplotlib.pyplot as plt
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
fig, axes = plt.subplots(2, 1, figsize=(10, 8))

# 自己相関
plot_acf(weekly, lags=60, ax=axes[0])
axes[0].set_title("Weekly Series ACF")

# 偏自己相関
plot_pacf(weekly, lags=60, ax=axes[1], method='ywm')
axes[1].set_title("Weekly Series PACF")

plt.tight_layout()
plt.show()

ACF/PACF

ACF（自己相関関数）は、時系列の値とラグ$k$だけ離れた値との相関を示します。PACF（偏自己相関関数）は、ラグ間の中間相関を除いた「純粋な」相関を示し、AR項の次数決定に役立ちます。

ACFからわかることは以下のようなことです。

ラグ1〜数週にかけて高い正の相関があり、直近の週の値が強く次週に影響すること
約25〜27週で相関が負のピーク、その後 50〜52 週で再び正のピークになるため、約52週の年周期を強く示唆すること
青色の背景色で表される95%信頼区間を大きく逸脱しており、ホワイトノイズではない構造的な相関が存在すること

また、PACFにはラグ1に0.8程度の大きなピークがあり、ラグ2以降は相関が急速に減衰するため、前週「のみ」が効く構造が強い（AR(1)）ことが読み取れます。

SARIMAモデル

SARIMAモデル（Seasonal ARIMA）は、非定常過程を差分で定常化しつつ、季節成分まで含めた複雑な時系列パターンを統計的にモデル化します。異常検知では「予測値と実測値のズレ」を残差として扱うことで、大きく外れた観測を検出できます。

SARIMAを理解する上では、そのSARIMAの構成要素となる以下のモデルを順に押さえる必要があります。

AR（自己回帰）モデル

過去の自分自身の値を使って現在の値を予測するモデル。

例：$y_t = \phi_1 y_{t-1} + \phi_2 y_{t-2} + \cdots + \epsilon_t$

MA（移動平均）モデル

過去の誤差（ノイズ）を使って現在の値を予測するモデル。

例： $y_t = \theta_1 \epsilon_{t-1} + \theta_2 \epsilon_{t-2} + \cdots + \epsilon_t$

ARMAモデル

ARとMAを組み合わせたモデル。ただし「定常性がある」ことが前提。

ARIMAモデル（自己回帰和移動平均モデル）

非定常なデータに対して、差分を取ることで定常化してからARMAを適用するモデル。

例： $\Delta y_t = \phi_1 \Delta y_{t-1} + \cdots + \theta_1 \epsilon_{t-1} + \cdots + \epsilon_t$

差分の回数（階数）を「d」と呼び、$\text{ARIMA}(p,d,q)$ と表記されます。

SARIMA

SARIMAはARIMA を拡張し、季節性（たとえば週次・年次などの周期性）を加味できるようにしたモデルで、次のように表記されます。

$\text{SARIMA}(p,d,q) \times (P,D,Q,m)$

$p,d,q$：通常の ARIMA の次数（非季節成分）
$P,D,Q$：季節性ARIMAの次数（季節成分）
$m$：季節の周期（例：日次で年次季節性なら 365）

ハイパーパラメータの選定

SARIMAのパラメータである$p,d,q,P,D,Q,s$を決めねばモデルは構築できません。どのパラメータが適切かは（PACFから少しわかるにしても）決めるのがなかなか難しいため、ここではAIC（Akaike Information Criterion）と呼ばれる基準を使い、総当たりで探索します。

AICは時系列モデルや回帰モデルの選定において「どのモデルが最もよいか」を判断するための指標で、モデルのあてはまりの良さ（対数尤度）とモデルの複雑さ（パラメータ数）のバランスを考慮し、過学習を防ぎつつ予測力の高いモデルを選ぶものです。基本的には値が小さいほど良いモデルです。

import itertools
import numpy as np
import statsmodels.api as sm

# 探索範囲設定
ps = range(0, 2)
qs = range(0, 2)
Ps = range(0, 2)
Qs = range(0, 2)
d, D, m = 1, 1, 52

best_aic = np.inf
best_order = None
best_seasonal_order = None

for p, q, P, Q in itertools.product(ps, qs, Ps, Qs):
    try:
        print("Trying parameters:", p, q, P, Q)
        mod = sm.tsa.SARIMAX(
            weekly,
            order=(p, d, q),
            seasonal_order=(P, D, Q, m),
            enforce_stationarity=False,
            enforce_invertibility=False,
            low_memory=True
        ).fit()
        aic = mod.aic
        print("aic=", aic)
        if aic < best_aic:
            best_aic = aic
            best_order = (p, d, q)
            best_seasonal_order = (P, D, Q, m)
    except Exception as e:
        print(f"An error occurred: {e!r}")

print("Best AIC:", best_aic)
print("Order:", best_order)
print("Seasonal order:", best_seasonal_order)

このスクリプトの結果は次のようになり、$\text{SARIMA}(1,1,1) \times (0,1,1,52)$のモデルが適当であろうということになりました。

Best AIC: 1655.0679118864932
Order: (1, 1, 1)
Seasonal order: (0, 1, 1, 52)

予測モデルの構築

それではSARIMAのモデルを構築してみます。

import statsmodels.api as sm

model = sm.tsa.SARIMAX(
    endog=weekly,
    order=(1, 1, 1),
    seasonal_order=(0, 1, 1, 52),
    enforce_stationarity=False,
    enforce_invertibility=False
)

results = model.fit(disp=False)

これを元に、1990年〜1992年までの２年間のデータを予測してみましょう。

import pandas as pd
import matplotlib.pyplot as plt

# 予測開始インデックスを取得
base = pd.to_datetime("1990-01-07")
start_idx = weekly.index[weekly.index >= base][0]

# 予測ステップ数（2年＝104週）
steps = 104

# 予測
end_idx = start_idx + pd.Timedelta(weeks=steps-1)
pred = results.get_prediction(start=start_idx, end=end_idx, dynamic=False)

pred_mean = pred.predicted_mean
pred_ci   = pred.conf_int()

fig, ax = plt.subplots(figsize=(12, 6))
obs_plot = weekly.loc[start_idx:end_idx]
obs_plot.plot(ax=ax, label='observed')

pred_mean.plot(ax=ax, label='forecast (2 years)', color='orange')
ax.fill_between(
    pred_ci.index,
    pred_ci.iloc[:, 0],
    pred_ci.iloc[:, 1],
    color='lightgrey',
    alpha=0.5,
    label='95% CI'
)

ax.set_xlim([start_idx, start_idx + pd.Timedelta(weeks=steps-1)])
ax.set_xlabel('Date')
ax.set_ylabel('Weekly Temperature')
ax.legend()
plt.show()