RStan:ローカルレベルモデル+カルマンフィルタ-パラメーター未知-季節性データ

# シミュレーションの条件設定

seed <- 20250628
set.seed(seed)

T <- 120 # データ期間 (10年 * 12ヶ月)
S <- 12 # 季節周期 (12ヶ月)

# 真のパラメータ
mu <- numeric(T)
season <- numeric(T)
y <- numeric(T)

# ノイズの標準偏差
sigma_w_true <- 15 # システムノイズ(レベル): 競合やメディア露出などによる変動
sigma_s_true <- 5 # システムノイズ(季節): 季節パターンの微小な変動
sigma_v_true <- 50 # 観測ノイズ: 天候などによるランダムな変動

# 初期値
mu[1] <- 300 # 初期の売上レベル
# 初期の季節効果（合計が0になるようにsin波で作成）
initial_season <- 150 * sin(2 * pi * (1:S) / S)
initial_season <- initial_season - mean(initial_season)
season[1:S] <- initial_season

# 時系列データの生成
for (t in 2:T) {
  # レベルの更新
  mu[t] <- mu[t - 1] + rnorm(1, mean = 0, sd = sigma_w_true)

  # 季節効果の更新 (t > S の場合)
  if (t > S) {
    # 過去(S-1)個の季節効果の合計にマイナスをつけたものが期待値
    expected_season <- -sum(season[(t - S + 1):(t - 1)])
    season[t] <- expected_season + rnorm(1, mean = 0, sd = sigma_s_true)
  }
}

# 観測データの生成
y <- mu + season + rnorm(T, mean = 0, sd = sigma_v_true)

# 負の売上は不自然なので0にクリップ
y[y < 0] <- 0

# 生成したデータの可視化
sim_data <- data.frame(
  time = 1:T,
  y = y,
  true_level = mu,
  true_season = season
)

library(ggplot2)
ggplot(sim_data, aes(x = time)) +
  geom_line(aes(y = true_level, color = "真のレベル (μ)")) +
  geom_line(aes(y = y, color = "観測売上 (y)")) +
  labs(
    title = "シミュレーションデータ：ジェラート店の月間売上",
    x = "月",
    y = "売上 (万円)",
    color = "データ"
  ) +
  theme_bw()

# stanfit オブジェクトの読み込み
setwd(stan_output)
fit <- readRDS("stan_fit.rds")

# --- 結果の確認 ---
# 収束診断
print(fit, pars = c("s_w", "s_v", "s_s", "lp__"), digits_summary = 3)

# トレースプロットで収束を確認
rstan::traceplot(fit, pars = c("s_w", "s_v", "s_s"))

# --- 結果の可視化 ---
# MCMCサンプルの抽出
samples <- rstan::extract(fit)

# プロット用のデータフレームを作成
# シミュレーションデータ(sim_data)と推定結果を結合する
library(dplyr)
result_df <- sim_data %>%
  # 推定されたレベルの中央値と信用区間を追加
  mutate(
    mu_median = apply(samples$mu, 2, median),
    mu_lower = apply(samples$mu, 2, quantile, probs = 0.025),
    mu_upper = apply(samples$mu, 2, quantile, probs = 0.975),
    # 推定された季節性の中央値と信用区間を追加
    season_median = apply(samples$season, 2, median),
    season_lower = apply(samples$season, 2, quantile, probs = 0.025),
    season_upper = apply(samples$season, 2, quantile, probs = 0.975)
  )

# レベル成分をプロット
ggplot(result_df, aes(x = time)) +
  # 推定されたレベルの95%信用区間をリボンで表示
  geom_ribbon(aes(ymin = mu_lower, ymax = mu_upper, fill = "推定値の95%信用区間"), alpha = 0.3) +
  # 推定されたレベルの中央値を線で表示
  geom_line(aes(y = mu_median, color = "推定されたレベル (中央値)")) +
  # 真のレベルを破線で表示
  geom_line(aes(y = true_level, color = "真のレベル")) +

  # ラベルとタイトルの設定
  labs(
    title = "推定されたレベル成分と真の値の比較",
    x = "月 (Time)",
    y = "レベル (Level)",
    color = "系列",
    fill = ""
  ) +

  # 色と凡例の設定
  scale_color_manual(
    name = "系列",
    values = c("推定されたレベル (中央値)" = "darkorange", "真のレベル" = "black")
  ) +
  scale_fill_manual(
    name = "",
    values = c("推定値の95%信用区間" = "darkorange")
  ) +

  # テーマの設定
  theme_bw() +
  theme(legend.position = "bottom")


# 季節成分をプロット
ggplot(result_df, aes(x = time)) +
  # 推定された季節成分の95%信用区間をリボンで表示
  geom_ribbon(aes(ymin = season_lower, ymax = season_upper, fill = "推定値の95%信用区間"), alpha = 0.3) +
  # 推定された季節成分の中央値を線で表示
  geom_line(aes(y = season_median, color = "推定された季節成分 (中央値)")) +
  # 真の季節成分を破線で表示
  geom_line(aes(y = true_season, color = "真の季節成分")) +

  # ラベルとタイトルの設定
  labs(
    title = "推定された季節成分と真の値の比較",
    x = "月 (Time)",
    y = "季節効果 (Seasonal Effect)",
    color = "系列",
    fill = ""
  ) +

  # 色と凡例の設定
  scale_color_manual(
    name = "系列",
    values = c("推定された季節成分 (中央値)" = "steelblue", "真の季節成分" = "black")
  ) +
  scale_fill_manual(
    name = "",
    values = c("推定値の95%信用区間" = "steelblue")
  ) +

  # テーマの設定
  theme_bw() +
  theme(legend.position = "bottom")

# 観測値と推定されたレベル＋季節性のプロット
ggplot(result_df, aes(x = time)) +
  geom_line(aes(y = y, color = "観測売上")) +
  geom_line(aes(y = mu_median + season_median, color = "モデル予測 (μ+season)")) +
  labs(
    title = "観測値とモデルによる予測値の比較",
    x = "月", y = "売上 (万円)", color = ""
  ) +
  theme_bw() +
  theme(legend.position = "bottom")