ベイズ構造時系列｜時系列分析

🎓 レベル：標準　|　重要度：A（必須）

📎 前提：ローカルレベルとローカルトレンドモデル　|　数理：状態空間モデルの枠組み・接続：信用区間と事後予測分布（ベイズ・事後予測）・ベイズ線形回帰（ベイズ・係数に事前）

要点（BLUF）

ベイズ構造時系列＝系列を「トレンド＋季節＋（回帰）＋ノイズ」に状態空間で分解し、各成分の未知量（分散・係数）に事前分布を置いてベイズ推定するモデル。
得られるのは点推定でなく事後分布。トレンドや季節の各成分が「95%信用区間つき」で出て、予測も事後予測分布として不確実性込みで出ます。
ローカルレベルとローカルトレンドモデルが分散 $Q,H$ を最尤で点推定したのに対し、ここは事前を置いて事後分布で推定します。真のトレンド・季節を仕込んだ擬似系列で、成分の復元と予測区間を pymc で確かめます。

1. 構造時系列をベイズで分解する

構造時系列は、観測 $y_t$ を解釈できる成分の和として状態空間で書きます（状態空間モデルの枠組みの構造モデル）。最小構成はトレンド $\mu_t$ ＋季節 $\gamma_t$ ＋観測ノイズ $\varepsilon_t$ ：

y_t = \mu_t + \gamma_t + \varepsilon_t, \qquad \varepsilon_t \sim N(0,\, \sigma_{\text{obs}}^2)

トレンドはランダムウォーク（ローカルレベル、ローカルレベルとローカルトレンドモデル）：

\mu_t = \mu_{t-1} + \eta_t, \qquad \eta_t \sim N(0,\, \sigma_{\text{level}}^2)

季節は周期 $s$ のフーリエ級数で表します（少ない係数で滑らかな季節を作れる。Prophet と同じ発想、Prophet（分解＋ベイズ））：

\gamma_t = \sum_{k=1}^{K} \left[ a_k \sin\!\frac{2\pi k t}{s} + b_k \cos\!\frac{2\pi k t}{s} \right]

ここまではローカルレベルとローカルトレンドモデルと同じ「構造」です。ベイズが変えるのは推定の仕方——未知量すべてに事前分布を置きます：

\sigma_{\text{level}} \sim \text{HalfNormal}(1),\quad \sigma_{\text{obs}} \sim \text{HalfNormal}(2),\quad a_k, b_k \sim N(0,\, 5^2)

そして観測を条件にした事後分布 $p(\mu_{1:n}, \gamma_{1:n}, \sigma_{\text{level}}, \sigma_{\text{obs}}, a, b \mid y_{1:n})$ を MCMC（NUTS、ハミルトニアンモンテカルロとNUTS）でサンプリングします。最尤との違いは決定的です。

flowchart LR
    PR["各成分の事前：トレンド分散・観測分散・季節フーリエ係数"] --> M["状態空間モデル y_t = トレンド + 季節 + ノイズ"]
    D["観測データ y_t"] --> M
    M --> MC["MCMC（NUTS）で同時サンプリング"]
    MC --> PO["事後分布：トレンド・季節の各成分 + 分散パラメータ"]
    PO --> CI["成分ごとに 95%信用区間"]
    PO --> PP["事後予測分布で予測 + 予測区間"]

最尤（ローカルレベルとローカルトレンドモデル）：分散 $Q,H$ を「最も尤もらしい1点」に固定し、状態はカルマンで条件付き推定。パラメータの不確実性は予測区間に入りません。
ベイズ（本ノート）：分散も係数も分布のまま。パラメータの不確実性が成分・予測区間に自動で伝播します。さらに事前で「トレンドは滑らか」等の知識を入れられます。

2. 擬似系列で成分を復元し、予測する

コード①：ベイズ構造時系列でトレンド・季節を復元し事後予測

真のランダムウォークトレンド（ $\sigma_{\text{level}}=0.5$ ）＋フーリエ季節（ $5\sin(2\pi t/12)+2\cos(2\pi\cdot 2t/12)$ ）＋観測ノイズ（ $\sigma_{\text{obs}}=1.5$ ）を仕込んだ系列に、pymc で上のモデルを当てます。トレンドのランダムウォークは非中心化（funnel 回避）で書きます。トレンド・季節の事後平均と95%信用区間が真値を捉えるか、フーリエ係数を復元するか、事後予測がホールドアウトを当てるかを数値で確かめます（PyMC 6.0.1・ArviZ 1.2.0、Windows は cores=1）。

import warnings
warnings.simplefilter("ignore")
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import japanize_matplotlib  # 日本語ラベル用
import pymc as pm
import pytensor.tensor as pt
import arviz as az

# 真の構造：ランダムウォークのトレンド + フーリエ季節(周期12) + 観測ノイズ
np.random.seed(7)
period, n_total, h = 12, 132, 12
n = n_total - h                       # 訓練 120
t_all = np.arange(n_total)
sigma_level_true, sigma_obs_true = 0.5, 1.5
trend_true = np.zeros(n_total); trend_true[0] = 20.0
for t in range(1, n_total):
    trend_true[t] = trend_true[t-1] + np.random.normal(0, sigma_level_true)
seasonal_true = 5*np.sin(2*np.pi*t_all/period) + 2*np.cos(2*np.pi*2*t_all/period)
y_all = trend_true + seasonal_true + np.random.normal(0, sigma_obs_true, n_total)
y, t_train = y_all[:n], t_all[:n]

# フーリエ季節の計画行列（K=2 ハーモニクス＝係数4本）
K = 2
def fourier(tt, period, K):
    cols = []
    for k in range(1, K+1):
        cols += [np.sin(2*np.pi*k*tt/period), np.cos(2*np.pi*k*tt/period)]
    return np.column_stack(cols)
X = fourier(t_train, period, K)

# ベイズ構造時系列：トレンド(RW・非中心化) + 季節(フーリエ係数に事前) + ノイズ
with pm.Model() as model:
    sigma_level = pm.HalfNormal("sigma_level", sigma=1.0)   # トレンド変化の大きさ（事前）
    sigma_obs = pm.HalfNormal("sigma_obs", sigma=2.0)       # 観測ノイズ（事前）
    mu0 = pm.Normal("mu0", mu=20.0, sigma=10.0)
    z = pm.Normal("z", 0.0, 1.0, shape=n-1)                 # RW増分の標準化変数
    trend = pm.Deterministic("trend", mu0 + pt.concatenate([[0.0], pt.cumsum(sigma_level*z)]))
    beta = pm.Normal("beta", 0.0, 5.0, shape=2*K)           # 季節フーリエ係数（事前）
    seasonal = pm.Deterministic("seasonal", pt.dot(X, beta))
    pm.Normal("y_obs", mu=trend + seasonal, sigma=sigma_obs, observed=y)
    idata = pm.sample(draws=1000, tune=1000, chains=2, cores=1,
                      random_seed=11, target_accept=0.95, progressbar=False)

print(f"発散数 = {int(idata.sample_stats['diverging'].values.sum())}")
print(f"R-hat sigma_level = {float(az.rhat(idata)['sigma_level'].values):.3f} / "
      f"sigma_obs = {float(az.rhat(idata)['sigma_obs'].values):.3f}")

# 成分の事後（平均と95%信用区間）を真値と対比
trend_post = idata.posterior["trend"].values.reshape(-1, n)
tr_mean = trend_post.mean(0)
tr_lo, tr_hi = np.percentile(trend_post, [2.5, 97.5], axis=0)
cover_tr = np.mean((trend_true[:n] >= tr_lo) & (trend_true[:n] <= tr_hi))
seas_post = idata.posterior["seasonal"].values.reshape(-1, n)
se_mean = seas_post.mean(0)
beta_mean = idata.posterior["beta"].values.reshape(-1, 2*K).mean(0)
print(f"トレンド事後平均RMSE vs 真値 = {np.sqrt(np.mean((tr_mean-trend_true[:n])**2)):.3f}")
print(f"トレンド95%信用帯の真値被覆 = {cover_tr*100:.0f}%")
print(f"季節事後平均RMSE vs 真値   = {np.sqrt(np.mean((se_mean-seasonal_true[:n])**2)):.3f}")
print(f"季節フーリエ係数の事後平均 = {np.round(beta_mean,2)}  （真値 sin1=5, cos1=0, sin2=0, cos2=2）")

# 事後予測による予測（RWを前進 + 季節 + 観測ノイズ）。区間は事後予測サンプルの分位点
sl = idata.posterior["sigma_level"].values.ravel()
so = idata.posterior["sigma_obs"].values.ravel()
beta_s = idata.posterior["beta"].values.reshape(-1, 2*K)
last_trend = trend_post[:, -1]
S = sl.shape[0]
t_future = t_all[n:n_total]
Xf = fourier(t_future, period, K)
rng = np.random.default_rng(123)
ypred = np.empty((S, h))
for i in range(S):
    fut_tr = last_trend[i] + np.cumsum(rng.normal(0, sl[i], h))
    ypred[i] = fut_tr + Xf @ beta_s[i] + rng.normal(0, so[i], h)
pred_mean = ypred.mean(0)
pred_lo, pred_hi = np.percentile(ypred, [2.5, 97.5], axis=0)
test = y_all[n:n_total]
print(f"予測RMSE = {np.sqrt(np.mean((pred_mean-test)**2)):.3f}  "
      f"95%事後予測区間カバレッジ = {int(np.sum((test>=pred_lo)&(test<=pred_hi)))}/{h}")
print(f"事後予測区間幅 1期先 = {pred_hi[0]-pred_lo[0]:.2f}  12期先 = {pred_hi[-1]-pred_lo[-1]:.2f}（先ほど広がる）")

# 図：トレンド/季節の事後と予測区間
fig, ax = plt.subplots(2, 1, figsize=(9, 7))
ax[0].plot(t_train, y, color="0.7", lw=0.8, marker="o", ms=2.5, label="観測 y_t")
ax[0].plot(t_train, trend_true[:n], color="k", ls="--", lw=1.4, label="真のトレンド")
ax[0].plot(t_train, tr_mean, color="C0", lw=1.8, label="トレンド事後平均")
ax[0].fill_between(t_train, tr_lo, tr_hi, color="C0", alpha=0.25, label="トレンド95%信用区間")
ax[0].plot(t_future, pred_mean, color="C1", lw=1.8, label="予測（事後予測平均）")
ax[0].fill_between(t_future, pred_lo, pred_hi, color="C1", alpha=0.25, label="95%事後予測区間")
ax[0].plot(t_future, test, color="r", lw=1.0, marker="x", ms=4, label="実測（検証）")
ax[0].axvline(n-0.5, ls=":", color="k")
ax[0].set_xlabel("時点 t"); ax[0].set_ylabel("値"); ax[0].legend(fontsize=8, ncol=2)
ax[0].set_title("ベイズ構造時系列：トレンドの事後＋予測区間")
ax[1].plot(t_train, seasonal_true[:n], color="k", ls="--", lw=1.4, label="真の季節")
ax[1].plot(t_train, se_mean, color="C2", lw=1.6, label="季節事後平均")
se_lo, se_hi = np.percentile(seas_post, [2.5, 97.5], axis=0)
ax[1].fill_between(t_train, se_lo, se_hi, color="C2", alpha=0.25, label="季節95%信用区間")
ax[1].set_xlabel("時点 t"); ax[1].set_ylabel("季節成分"); ax[1].legend(fontsize=8)
ax[1].set_title("季節成分の復元（フーリエ係数の事後）")
plt.tight_layout(); plt.show()

出力：

発散数 = 0
R-hat sigma_level = 1.001 / sigma_obs = 1.003
トレンド事後平均RMSE vs 真値 = 0.613
トレンド95%信用帯の真値被覆 = 94%
季節事後平均RMSE vs 真値   = 0.253
季節フーリエ係数の事後平均 = [ 5.11 -0.26 -0.18  1.89]  （真値 sin1=5, cos1=0, sin2=0, cos2=2）
予測RMSE = 1.877  95%事後予測区間カバレッジ = 12/12
事後予測区間幅 1期先 = 6.61  12期先 = 8.36（先ほど広がる）

出力の意味：まず収束——発散数 0、R-hat も $1.001/1.003$ で問題なし（収束診断）。本題は成分の復元です。トレンドは事後平均 RMSE $0.613$ 、しかも95%信用帯が真の隠れトレンドを 94% 被覆（名目 95% とよく整合）——観測ノイズに埋もれた潜在トレンドを、不確実性つきで掘り出せました。季節は事後平均 RMSE $0.253$ 、フーリエ係数の事後平均が $[5.11,\ -0.26,\ -0.18,\ 1.89]$ で、真値 $[\text{sin}_1{=}5,\ \text{cos}_1{=}0,\ \text{sin}_2{=}0,\ \text{cos}_2{=}2]$ をきれいに当てています。予測は事後予測分布から作り、RMSE $1.877$ 、95%事後予測区間が 12/12 を被覆。区間幅は $1$ 期先 $6.61$ → $12$ 期先 $8.36$ とホライズンとともに広がる——ランダムウォークトレンドの将来不確実性が累積するためで、ローカルレベルとローカルトレンドモデルの予測区間が広がるのと同じ理屈です。点だけでなく「成分も予測も、信用区間つき」で出るのがベイズ構造時系列の値打ちです。

3. 04-03 の最尤との違い：点推定 vs 事後分布

同じ「構造（トレンド＋季節）」でも、推定の哲学が違います。

	ローカルレベルとローカルトレンドモデル（最尤）	本ノート（ベイズ）
分散 $Q,H$	1点に固定（最尤値）	事後分布のまま
状態（トレンド・季節）	$Q,H$ 既知としてカルマンで推定	パラメータと同時にサンプル
パラメータの不確実性	予測区間に入らない	予測区間に自動で伝播
事前知識	入れにくい	事前で「滑らかさ」等を入れられる
計算	速い（カルマン＋最適化）	重い（MCMC）

最尤の予測区間は「 $Q,H$ がぴったり当たっている」前提で計算されるので、短い系列で $Q,H$ の推定が怪しいときは楽観的になりがちです。ベイズは「 $Q,H$ がこのくらいばらつく」ことごと事後予測に織り込むので、区間が正直に広がります。この差は次ノート状態空間のMCMC（カルマン×ベイズ）で、 $Q,H$ の事後分布を最尤の点と直接重ねて可視化します。

4. 数式の直観

「分解」と「不確実性」を同時に得る：状態空間は系列を解釈できる成分（水準・傾き・季節）に割り、ベイズは各成分を分布で返す。「トレンドはどの辺で、どれくらい確信があるか」が信用区間で読めます。
フーリエ季節＝少ない係数で滑らかな周期：周期 $s$ を $\sin/\cos$ の和で表すと、 $K$ ハーモニクスで $2K$ 個の係数だけ。係数に正規事前を置くのはベイズ線形回帰（ベイズ線形回帰）そのもの。 $K$ を上げるほど複雑な季節を表せますが、事前が過適合を抑えます。
事前は「正則化」： $\sigma_{\text{level}}$ に小さめの事前を置けば「トレンドは滑らか」を促し、係数の事前は季節の振幅を抑える。これは正則化との関係の通り罰則付き推定と表裏一体です。
事後予測＝パラメータを積分消去した予測： $p(y_{n+h}\mid y_{1:n})=\int p(y_{n+h}\mid\theta)\,p(\theta\mid y_{1:n})\,d\theta$ 。MCMC ではパラメータ事後からサンプルし、各サンプルで将来を生成して集めるだけ（コード①の for ループ）。これが信用区間と事後予測分布の事後予測分布です。

⚠️ よくある誤解・落とし穴

「事前は結果に効かない」ではない：データが少ない・成分が識別しにくいと、事前が事後を強く引っ張ります。 $\sigma_{\text{level}}$ の事前を変えてトレンドの滑らかさが動かないか感度分析を（事前の選び方と感度分析）。本コードの結果は弱情報事前＋十分なデータゆえ安定しています。
「信用区間＝信頼区間」ではない：信用区間は「パラメータがこの区間にある事後確率95%」というベイズの直接解釈。頻度論の信頼区間とは意味が違います（信用区間と信頼区間）。
「ランダムウォークトレンドは中心化のまま回せる」ではない： $\sigma_{\text{level}}$ と増分がfunnel を作り、発散・低 ESS を招きます。コードのように非中心化（ $\mu_t=\mu_0+\sigma_{\text{level}}\sum z_k,\ z_k\sim N(0,1)$ ）が定石（階層モデルの実例と再パラメータ化）。
「事後予測区間は外挿でも万能」ではない：区間は「モデル構造が正しい」前提。トレンドが将来屈折する・季節が変わるといった構造変化は折り込めません。遠未来ほど名目より楽観的になり得ます（Prophet（分解＋ベイズ）の changepoint はこの一部に対処）。
「最尤より常に良い」ではない：ベイズは MCMC ぶん重く、収束診断（R-hat・発散）を毎回確認する手間が要ります。素早い点予測ならローカルレベルとローカルトレンドモデルの UnobservedComponents で十分なことも多い。不確実性を正直に出したいときにベイズ、が使い分けです。