ローカルレベルモデルは、状態空間モデルの中で最も簡単な設定といえるモデルです。
状態空間モデルの中でも、状態が全てガウス分布で表現でき、さらにシステム行列が線形である状態空間モデルを動的線形モデルや線形・ガウシアンモデルなどと表現しますが、ローカルレベルモデルは、これよりもさらに簡単な設定です。
今回は、ローカルレベルモデルを、Pythonのstatsmodelsライブラリを利用してサクッと実装してみます。
本記事の内容
- ローカルレベルモデルの数学的定義
- statsmodelsを利用した実装
- 合成データと実データに対するフィッティングと可視化
ローカルレベルモデルの数学的定義
ローカルレベルモデルは、以下の2つの方程式で定義されます。
観測方程式:
$$ y_t = \mu_t + \epsilon_t, \quad \epsilon_t \sim \mathcal{N}(0, \sigma_\epsilon^2) $$
状態遷移方程式:
$$ \mu_{t+1} = \mu_t + \eta_t, \quad \eta_t \sim \mathcal{N}(0, \sigma_\eta^2) $$
ここで、$y_t$ は時刻 $t$ における観測値、$\mu_t$ は時刻 $t$ における状態（レベル）、$\epsilon_t$ は観測ノイズ、$\eta_t$ はシステムノイズです。
これを一般的な状態空間モデルの形式で書くと、
$$ \begin{align} y_t &= \bm{Z} \alpha_t + \epsilon_t \\ \alpha_{t+1} &= \bm{T} \alpha_t + \bm{R} \eta_t \end{align} $$
ローカルレベルモデルでは、$\bm{Z} = 1$, $\bm{T} = 1$, $\bm{R} = 1$ となり、最もシンプルな設定です。
このモデルはカルマンフィルタで状態推定を行うことができます。推定すべきパラメータは観測ノイズの分散 $\sigma_\epsilon^2$ とシステムノイズの分散 $\sigma_\eta^2$ の2つだけです。
合成データでの実装
まずは合成データを生成して、ローカルレベルモデルの挙動を確認しましょう。
import numpy as np
import matplotlib.pyplot as plt
import statsmodels.api as sm
np.random.seed(42)
# パラメータ設定
n = 200
sigma_eta = 0.5 # システムノイズの標準偏差
sigma_eps = 1.0 # 観測ノイズの標準偏差
# 状態と観測の生成
mu = np.zeros(n)
y = np.zeros(n)
mu[0] = 0.0
y[0] = mu[0] + np.random.normal(0, sigma_eps)
for t in range(1, n):
mu[t] = mu[t-1] + np.random.normal(0, sigma_eta)
y[t] = mu[t] + np.random.normal(0, sigma_eps)
# 可視化
plt.figure(figsize=(12, 4))
plt.plot(y, 'o', markersize=2, alpha=0.5, label='Observations (y_t)')
plt.plot(mu, 'r-', linewidth=2, label='True state (mu_t)')
plt.title("Local Level Model: Synthetic Data")
plt.xlabel("Time")
plt.ylabel("Value")
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
赤線が真の状態 $\mu_t$ で、青い点が観測値 $y_t$ です。観測値は状態にノイズが加わったものなので、状態の周りにばらついていることが分かります。
statsmodelsでモデルをフィッティング
statsmodelsの UnobservedComponents クラスを使って、ローカルレベルモデルをフィッティングします。
import numpy as np
import matplotlib.pyplot as plt
import statsmodels.api as sm
np.random.seed(42)
# 合成データの生成（前セクションと同じ）
n = 200
sigma_eta = 0.5
sigma_eps = 1.0
mu_true = np.zeros(n)
y = np.zeros(n)
mu_true[0] = 0.0
y[0] = mu_true[0] + np.random.normal(0, sigma_eps)
for t in range(1, n):
mu_true[t] = mu_true[t-1] + np.random.normal(0, sigma_eta)
y[t] = mu_true[t] + np.random.normal(0, sigma_eps)
# ローカルレベルモデルの構築とフィッティング
model = sm.tsa.UnobservedComponents(y, level='local level')
result = model.fit(disp=False)
# 推定結果の表示
print(result.summary())
# 推定された状態の取得
estimated_level = result.level.smoothed
# 可視化
fig, axes = plt.subplots(2, 1, figsize=(12, 8))
# 上段: 観測値・真の状態・推定された状態
axes[0].plot(y, 'o', markersize=2, alpha=0.4, label='Observations')
axes[0].plot(mu_true, 'g-', linewidth=2, alpha=0.7, label='True state')
axes[0].plot(estimated_level, 'r-', linewidth=2, label='Estimated state (smoothed)')
axes[0].set_title("Local Level Model: State Estimation")
axes[0].set_xlabel("Time")
axes[0].set_ylabel("Value")
axes[0].legend()
axes[0].grid(True, alpha=0.3)
# 下段: フィルタリングとスムージングの比較
filtered_level = result.level.filtered
axes[1].plot(y, 'o', markersize=2, alpha=0.3, label='Observations')
axes[1].plot(filtered_level, 'b-', linewidth=1.5, label='Filtered state')
axes[1].plot(estimated_level, 'r-', linewidth=1.5, label='Smoothed state')
axes[1].set_title("Filtered vs Smoothed State")
axes[1].set_xlabel("Time")
axes[1].set_ylabel("Value")
axes[1].legend()
axes[1].grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
print(f"\n推定されたパラメータ:")
print(f" 観測ノイズ分散 (sigma_eps^2): {result.params[0]:.4f} (真値: {sigma_eps**2:.4f})")
print(f" システムノイズ分散 (sigma_eta^2): {result.params[1]:.4f} (真値: {sigma_eta**2:.4f})")
上段のグラフでは、推定された状態（赤）が真の状態（緑）に近い値を追跡していることが確認できます。下段では、フィルタリング（過去のデータのみ使用）とスムージング（全データ使用）の違いが見て取れます。スムージングの方がより滑らかで真の状態に近い推定を与えます。
実データでの適用: 航空旅客数データ
実データとして、有名な航空旅客数データセット（AirPassengers）にローカルレベルモデルを適用してみましょう。
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import statsmodels.api as sm
# AirPassengersデータの読み込み
data = sm.datasets.get_rdataset("AirPassengers")
passengers = data.data
passengers.columns = ['time', 'value']
y = passengers['value'].values
# 対数変換（分散安定化のため）
y_log = np.log(y)
# ローカルレベルモデルのフィッティング
model = sm.tsa.UnobservedComponents(y_log, level='local level')
result = model.fit(disp=False)
# 状態推定
level_smoothed = result.level.smoothed
# 可視化
fig, axes = plt.subplots(2, 1, figsize=(12, 8))
# 上段: 対数スケール
axes[0].plot(y_log, 'b-', linewidth=1, alpha=0.7, label='log(AirPassengers)')
axes[0].plot(level_smoothed, 'r-', linewidth=2, label='Estimated level (smoothed)')
axes[0].set_title("Local Level Model on AirPassengers (log scale)")
axes[0].set_xlabel("Month index")
axes[0].set_ylabel("log(Passengers)")
axes[0].legend()
axes[0].grid(True, alpha=0.3)
# 下段: 元のスケール
axes[1].plot(y, 'b-', linewidth=1, alpha=0.7, label='AirPassengers')
axes[1].plot(np.exp(level_smoothed), 'r-', linewidth=2, label='Estimated level')
axes[1].set_title("Local Level Model on AirPassengers (original scale)")
axes[1].set_xlabel("Month index")
axes[1].set_ylabel("Passengers")
axes[1].legend()
axes[1].grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
航空旅客数データはトレンドと季節性を持つため、ローカルレベルモデルだけでは季節変動を捉えることができません。推定されたレベルはトレンド成分のみを追跡しています。季節性も扱いたい場合は、ローカルリニアトレンドモデルや構造時系列モデル（季節成分を含むモデル）を利用する必要があります。
ローカルレベルモデルの限界
ローカルレベルモデルの特徴と限界を整理しておきましょう。
ローカルレベルモデルは、状態が「ランダムウォーク + 観測ノイズ」という最もシンプルな構造を仮定しています。そのため、以下のケースでは、より複雑なモデルが必要になります。
- トレンド（傾き）がある場合: ローカルリニアトレンドモデルを使う
- 季節性がある場合: 季節成分を追加した構造時系列モデルを使う
- 非線形なダイナミクス: 拡張カルマンフィルタやパーティクルフィルタを使う
まとめ
本記事では、状態空間モデルの中で最もシンプルなローカルレベルモデルをPythonのstatsmodelsで実装しました。
- ローカルレベルモデルはランダムウォーク状態 + 観測ノイズという最小構成の状態空間モデル
- statsmodelsの
UnobservedComponentsで簡単に実装できる - フィルタリング（過去データのみ使用）とスムージング（全データ使用）の2種類の状態推定が得られる
- 季節性やトレンドの傾きを扱うには、より複雑なモデルへの拡張が必要