🕐 1. 時系列とは

時間順に並んだ観測値の列 $\{y_t\}_{t=1}^T$。 観測の順序と時間間隔に意味がある。 通常の機械学習の独立同分布仮定が成り立たない。

典型的な構成要素

• トレンド：長期的な傾き
• 季節性：周期的変動（年・週・日次）
• 循環：景気循環のような非周期変動
• 不規則変動：ノイズ

🧩 2. 時系列分解

加法モデル：$y_t = T_t + S_t + R_t$（トレンド + 季節 + 残差）

乗法モデル：$y_t = T_t \cdot S_t \cdot R_t$

1
2
3
4
5
6
7
8

import pandas as pd
from statsmodels.tsa.seasonal import seasonal_decompose

# 例：気象庁公開の月次気温データを想定（実データURL読み替え可）
ts = pd.read_csv('data/raw/temperature_monthly.csv',
                  index_col='date', parse_dates=True)['temp']
res = seasonal_decompose(ts, model='additive', period=12)
res.plot()

STL 分解（Seasonal-Trend Loess）

季節性が時間変動する場合に有効。 ロバストオプション付き。

⚖️ 3. 定常性

定常時系列：平均・分散・自己相関が時間によらず一定。 多くのモデルが定常性を要求する。

3.1 ADF（拡張ディッキー–フラー）検定

帰無仮説：単位根あり（非定常）。 p < 0.05 で棄却＝定常。

1
2
3

from statsmodels.tsa.stattools import adfuller
res = adfuller(ts)
print(f'ADF統計量: {res[0]:.3f}, p値: {res[1]:.4f}')

3.2 差分・対数差分

• 1階差分：$\Delta y_t = y_t - y_{t-1}$（トレンド除去）
• 季節差分：$y_t - y_{t-s}$（季節除去）
• 対数差分：$\Delta \log y_t$（対数収益率に対応）

1
2

ts_diff = ts.diff().dropna()
ts_log_diff = ts.apply(lambda x: x).pipe(lambda s: s.apply(__import__('numpy').log)).diff().dropna()

KPSS検定

ADF と逆で、 帰無仮説が定常。 両方を併用するのが推奨（Hyndman 流）。

📊 4. 自己相関関数 (ACF) と偏自己相関 (PACF)

$$\rho_k = \mathrm{Corr}(y_t, y_{t-k})$$

ACF は「ラグ k での相関」、 PACF は「他のラグの影響を取り除いた」相関。

ACF/PACF の読み方

1
2
3
4
5
6

from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
import matplotlib.pyplot as plt
fig, axes = plt.subplots(1, 2, figsize=(14, 4))
plot_acf(ts_diff, lags=24, ax=axes[0])
plot_pacf(ts_diff, lags=24, ax=axes[1])
plt.show()

🔬 5. ARIMA モデル

5.1 AR(p) モデル

$$y_t = c + \phi_1 y_{t-1} + \cdots + \phi_p y_{t-p} + \varepsilon_t$$

$\phi_i$ は「ファイ・サブ・アイ」自己回帰係数。 過去の自分を線形結合。

5.2 MA(q) モデル

$$y_t = \mu + \varepsilon_t + \theta_1 \varepsilon_{t-1} + \cdots + \theta_q \varepsilon_{t-q}$$

$\theta_j$ は「シータ・サブ・ジェイ」移動平均係数。 過去のショックを線形結合。

5.3 ARMA(p, q)

AR と MA の組合せ。 定常データ向け。

5.4 ARIMA(p, d, q)

d 階差分してから ARMA を適用。 非定常データに対応。

1
2
3
4
5
6

from statsmodels.tsa.arima.model import ARIMA
model = ARIMA(ts, order=(1, 1, 1))
fit = model.fit()
print(fit.summary())
forecast = fit.forecast(steps=12)
print(forecast)

5.5 SARIMA(p,d,q)(P,D,Q,s)

季節成分を別の (P,D,Q) で表現。 s は季節周期（月次なら 12、 週次なら 7）。

1
2
3
4

from statsmodels.tsa.statespace.sarimax import SARIMAX
model = SARIMAX(ts, order=(1, 1, 1), seasonal_order=(1, 1, 1, 12))
fit = model.fit(disp=False)
print(fit.aic)

5.6 自動次数選択（pmdarima）

1
2
3
4

import pmdarima as pm
auto = pm.auto_arima(ts, seasonal=True, m=12, stepwise=True, suppress_warnings=True)
print(auto.summary())
forecast = auto.predict(n_periods=12)

📈 6. 指数平滑法

古典的な予測手法。 過去ほど指数的に重みを下げる。

• 単純指数平滑：レベルのみ
• Holt法：レベル + トレンド
• Holt-Winters：レベル + トレンド + 季節（加法/乗法）
• ETS：状態空間表現で AICによる自動選択

1
2
3

from statsmodels.tsa.holtwinters import ExponentialSmoothing
hw = ExponentialSmoothing(ts, trend='add', seasonal='add', seasonal_periods=12).fit()
forecast = hw.forecast(12)

🌌 7. 状態空間モデル

観測 $y_t = Z \alpha_t + \varepsilon_t$、 状態 $\alpha_{t+1} = T\alpha_t + \eta_t$。 Kalman フィルタで隠れ状態を再帰的に推定。

• 欠損値を自動で扱える
• 構造的時系列モデル（ローカルレベル・ローカル線形トレンド）
• 動的線形モデル (DLM)
• UCM (Unobserved Components Model)

1
2
3
4

from statsmodels.tsa.statespace.structural import UnobservedComponents
ucm = UnobservedComponents(ts, level='local linear trend', seasonal=12)
fit = ucm.fit()
print(fit.summary())

🔮 8. Prophet（Meta）

「トレンド + 季節 + 休日効果」のシンプルな加法モデル。 欠損・外れ値に頑健。 ビジネス時系列で人気。

1
2
3
4
5
6
7

from prophet import Prophet
df_p = pd.DataFrame({'ds': ts.index, 'y': ts.values})
m = Prophet(yearly_seasonality=True, weekly_seasonality=False)
m.fit(df_p)
future = m.make_future_dataframe(periods=12, freq='M')
fcst = m.predict(future)
m.plot(fcst)

🔗 9. VAR（ベクトル自己回帰）

複数時系列を同時にモデル化：$\mathbf{y}_t = \mathbf{c} + \sum_i \mathbf{A}_i \mathbf{y}_{t-i} + \boldsymbol{\varepsilon}_t$。

応用：因果関係（Granger因果）、 インパルス応答、 ショック分解。

1
2
3
4

from statsmodels.tsa.api import VAR
multi_ts = pd.DataFrame({'temp': ts, 'humidity': ts2})
var = VAR(multi_ts).fit(maxlags=12, ic='aic')
print(var.summary())

📊 10. 予測評価指標

• MAE：平均絶対誤差
• RMSE：二乗平均平方根誤差
• MAPE：平均絶対パーセント誤差。 ただし y=0 で破綻
• sMAPE：対称版 MAPE
• MASE：単純予測（ナイーブ）との比較。 推奨 (Hyndman)

TimeSeriesSplit でのバックテスト

 1
 2
 3
 4
 5
 6
 7
 8
 9
10

from sklearn.model_selection import TimeSeriesSplit
import numpy as np
tscv = TimeSeriesSplit(n_splits=5, test_size=12)
errors = []
y = ts.values
for tr, te in tscv.split(y):
    model = ARIMA(y[tr], order=(1,1,1)).fit()
    pred = model.forecast(steps=len(te))
    errors.append(np.mean(np.abs(y[te] - pred)))
print(f'CV MAE: {np.mean(errors):.3f}')

📊 11. モデル比較

⚠️ 12. よくある落とし穴

🏋️ 13. 練習問題

注意：SSDSE-B は時系列ではないので、 別の時系列実データ（例：気象庁公開の日次気温、 e-Stat の月次小売販売、 政府統計局の月次失業率など）を使用してください。

1
2
3

from statsmodels.tsa.seasonal import STL
res = STL(ts, period=12, robust=True).fit()
res.plot()

1
2
3

from statsmodels.tsa.stattools import adfuller, kpss
print('ADF p:', adfuller(ts)[1])
print('KPSS p:', kpss(ts, nlags='auto')[1])

1
2
3
4
5
6
7
8

import pmdarima as pm
import numpy as np
train, test = ts[:-12], ts[-12:]
auto = pm.auto_arima(train, seasonal=True, m=12)
pred = auto.predict(n_periods=12)
naive = train.iloc[-12:].values  # 直近1年で繰り返し
mase = np.mean(np.abs(test - pred)) / np.mean(np.abs(np.diff(train)))
print(f'MASE = {mase:.3f}')

📝 14. 報告フォーマット

❌ NG例

「ARIMA を当てはめたら良い予測ができました。」

✅ OK例

「月次データ（2015-01 から 2025-12、 N=132）に対し、 1階差分で定常化（ADF p < .001）。 ACF/PACF と AIC 最小化により SARIMA(1,1,1)(1,1,1,12) を選択。 12-期 rolling-origin バックテストで MASE = 0.81（< 1 でナイーブより優れる）、 95% 予測区間のカバレッジ 92%（名目 95%）。 比較として Prophet（MASE=0.85）、 ETS（0.78）、 LightGBM with lag features（0.74）を実施し、 LightGBM を採用。」

🐍 15. ライブラリ早見表

📜 16. 時系列分析の歴史

• 1927：Yule の自己回帰モデル
• 1933：Khintchine, Wold の確率過程理論
• 1960：Kalman フィルタ
• 1970：Box & Jenkins の ARIMA 教科書 — 時系列モデリングを体系化
• 1980：Engle の ARCH（金融ボラティリティ）
• 1990：状態空間モデル普及（Harvey, Durbin & Koopman）
• 2017：Facebook の Prophet 公開
• 2020年代：Transformer ベース時系列モデル（TFT, Informer, PatchTST）

💼 17. 実務応用

• 需要予測：小売・製造の在庫管理
• 金融：株価・為替・ボラティリティ予測（GARCH 系）
• 気象・エネルギー：気温・電力需要・発電量予測
• マーケティング：トラフィック・コンバージョン予測
• IoT：センサー異常検知
• 疫学：感染症の流行予測（COVID-19）
• マクロ経済：GDP・物価・失業率の予測

✅ 18. 時系列分析チェックリスト

• □ プロットで全体の概要を把握したか？
• □ 欠損・外れ値を確認したか？
• □ 季節成分・トレンドを分解したか？
• □ ADF/KPSS で定常性を確認したか？
• □ 必要なら差分・対数変換したか？
• □ ACF/PACF を確認したか？
• □ 複数モデル（ARIMA / ETS / Prophet）を比較したか？
• □ TimeSeriesSplit でバックテストしたか？
• □ MASE 等の適切な評価指標を使ったか？
• □ 予測区間を出したか？
• □ 残差の自己相関 (Ljung–Box) を確認したか？

🩺 19. 残差診断

ARIMA 等のフィット後、 残差が「ホワイトノイズ」になっているかを確認。 残差にパターンが残っていれば、 モデルが情報を取り切れていない。

Ljung–Box 検定

残差の自己相関がすべて 0 か？ 帰無仮説：H₀ = 自己相関なし。 p > 0.05 で「ホワイトノイズ」と判断可。

1
2
3

from statsmodels.stats.diagnostic import acorr_ljungbox
lb = acorr_ljungbox(fit.resid, lags=[12], return_df=True)
print(lb)

残差プロット

• 残差の時系列プロット → 構造的パターンが残っていないか
• 残差の ACF → 自己相関の残存
• 残差のQQプロット → 正規性
• 残差ヒストグラム → 分布形

🔗 関連用語

• 時系列（個別）
• 自己相関
• 移動平均
• モデル選択 — TimeSeriesSplit
• 評価指標
• 回帰
• パネル・因果推論

🧮 SSDSE-B を使った時系列分析の実例

SSDSE-B のパネル構造（47都道府県 × 2003-2020年）から、 例えば東京都の「人口」や沖縄県の「合計特殊出生率」を取り出すと、 単一都道府県の年次時系列が得られます。 ここでは東京都の総人口（千人）を題材に、 ARIMA・指数平滑・トレンド分解の実値計算を見ます。

1
2
3
4
5
6
7
8
9

import pandas as pd
df = pd.read_csv('data/raw/SSDSE-B-2023.csv', encoding='shift_jis', header=[0,1])
df.columns = ['_'.join(c).strip() for c in df.columns]
tokyo = df[df['都道府県_Prefecture'] == '東京都'].set_index('年度_Year')
y = tokyo['総人口_Total population'] / 1000  # 千人単位
print(y.head())
# 2003 年: 12,420 千人
# 2010 年: 13,160 千人
# 2020 年: 14,048 千人 — 緩やかに増加

1
2
3
4
5

from statsmodels.tsa.seasonal import STL
# 年次データなのでトレンド + 残差のみ（周期なし）
stl = STL(y, period=None, robust=True).fit()
print(f'トレンド成分: 2003={stl.trend.iloc[0]:.0f}, 2020={stl.trend.iloc[-1]:.0f}')
print(f'残差 SD: {stl.resid.std():.1f}')

1
2
3
4
5
6
7
8
9

from statsmodels.tsa.stattools import adfuller
res = adfuller(y)
print(f'ADF stat = {res[0]:.3f}, p = {res[1]:.4f}')
# 例：ADF = -0.5, p = 0.88 → 単位根あり（非定常）

# 1階階差で再検定
res2 = adfuller(y.diff().dropna())
print(f'diff ADF stat = {res2[0]:.3f}, p = {res2[1]:.4f}')
# 例：ADF = -3.2, p = 0.02 → 階差で定常 → I(1) → ARIMA(p,1,q) 候補

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11

from statsmodels.tsa.arima.model import ARIMA
model = ARIMA(y, order=(1, 1, 0)).fit()
print(model.summary())
# AR(1) 係数 ≈ 0.4, σ² ≈ 1500 程度

forecast = model.get_forecast(steps=5)  # 2021-2025 を予測
ci = forecast.conf_int(alpha=0.05)
print(forecast.predicted_mean)
print(ci)
# 2021: 14,090 千人 [13,990, 14,190]（CI 幅 ±100）
# 2025: 14,180 千人 [13,890, 14,470]（先になるほど CI 拡大）

1
2
3
4
5

from statsmodels.tsa.holtwinters import ExponentialSmoothing
hw = ExponentialSmoothing(y, trend='add', seasonal=None).fit()
print(f'α (level) = {hw.params["smoothing_level"]:.3f}')
print(f'β (trend) = {hw.params["smoothing_trend"]:.3f}')
print(hw.forecast(5))  # ARIMA とほぼ同じ予測値が得られる

🐍 実装バリエーション — statsmodels / sktime / Prophet / Darts

時系列ライブラリは選択肢が豊富。 用途・速度・拡張性で使い分けます。

(A) statsmodels — 古典統計の標準

1
2
3
4

from statsmodels.tsa.arima.model import ARIMA
from statsmodels.tsa.statespace.sarimax import SARIMAX
from statsmodels.tsa.holtwinters import ExponentialSmoothing
# 古典的 ARIMA・SARIMAX・状態空間モデルが揃う

(B) pmdarima — auto_arima の決定版

1
2
3
4

from pmdarima import auto_arima
model = auto_arima(y, seasonal=False, stepwise=True, trace=True)
print(model.summary())
# (p,d,q) を AIC 最小で自動探索

(C) Prophet — トレンド + 休日 + 季節を柔軟に

1
2
3
4
5
6

from prophet import Prophet
df_p = pd.DataFrame({'ds': pd.to_datetime(y.index, format='%Y'), 'y': y.values})
m = Prophet(yearly_seasonality=False, daily_seasonality=False).fit(df_p)
future = m.make_future_dataframe(periods=5, freq='Y')
forecast = m.predict(future)
print(forecast[['ds', 'yhat', 'yhat_lower', 'yhat_upper']].tail(5))

(D) sktime — scikit-learn API で時系列

1
2
3
4
5
6
7

from sktime.forecasting.arima import ARIMA as SKARIMA
from sktime.forecasting.model_selection import temporal_train_test_split
y_train, y_test = temporal_train_test_split(y, test_size=3)
forecaster = SKARIMA(order=(1,1,0))
forecaster.fit(y_train)
y_pred = forecaster.predict(fh=[1, 2, 3])
print(y_pred)

(E) Darts — 機械学習・深層学習も統一API

1
2
3
4
5

from darts import TimeSeries
from darts.models import ExponentialSmoothing as DartsES, NBEATSModel
ts = TimeSeries.from_series(y)
es = DartsES(); es.fit(ts); print(es.predict(5))
# NBEATSModel など深層学習モデルも同じインターフェース

(F) scikit-learn TimeSeriesSplit — 時系列CV

1
2
3
4
5
6

from sklearn.model_selection import TimeSeriesSplit
from sklearn.metrics import mean_absolute_error
tscv = TimeSeriesSplit(n_splits=5)
for train_idx, test_idx in tscv.split(y):
    y_tr, y_te = y.iloc[train_idx], y.iloc[test_idx]
    # 各 fold でモデル訓練・評価。 通常CV と違い「過去→未来」順を維持

⚠️ 追加の落とし穴 — 時系列分析の実務

🔗 関連用語（拡張マップ）

📚 前提となる用語（4個以上）

🔀 並列・関連トピック（4個以上）

🚀 発展・応用先（4個以上）