Rの関数:dlmModReg {dlm}

# 必要なパッケージをロードします
library(ggplot2)
library(dplyr)
library(tidyr)

# 再現性のために乱数シードを設定します
seed <- 20251001
set.seed(seed)

# シミュレーションの設定
time_len <- 1000
v_var_true <- 0.8 # 観測ノイズの真の分散
dw0_true <- 0.1 # 切片係数の変動の大きさ（システムノイズの分散）
dw1_true <- 0.2 # xの係数の変動の大きさ（システムノイズの分散）

# 説明変数 X を作成します (sinカーブ)
x_t <- sin(2 * pi * (1:time_len) / 50)

# 真の時変係数を作成します
# 切片項 (beta_0,t) はランダムウォークに従います
b0_true <- cumsum(rnorm(time_len, mean = 0, sd = sqrt(dw0_true)))
# xの係数 (beta_1,t) もランダムウォークに従います
b1_true <- cumsum(rnorm(time_len, mean = 0, sd = sqrt(dw1_true)))
# b1_true の平均を調整して、変動を分かりやすくします
b1_true <- b1_true - mean(b1_true) + 2

# 観測ノイズ v_t を作成します
v_t <- rnorm(time_len, mean = 0, sd = sqrt(v_var_true))

# 観測値 y_t を生成します
y_t <- b0_true + b1_true * x_t + v_t

# 作成したデータをデータフレームにまとめます
sim_data <- data.frame(
  time = 1:time_len,
  y = y_t,
  x = x_t,
  b0_true = b0_true,
  b1_true = b1_true
)

# 生成した観測値と真の係数を確認します
cat("--- サンプルデータの一部を確認 ---\n")
print(head(sim_data))

# データをプロットして確認します
p_data <- ggplot(sim_data, aes(x = time)) +
  geom_line(aes(y = y, color = "観測値 y")) +
  geom_line(aes(y = b0_true, color = "真の切片係数")) +
  geom_line(aes(y = b1_true, color = "真のxの係数")) +
  labs(
    title = "生成したサンプルデータ",
    x = "時間",
    y = "値",
    color = "系列"
  ) +
  theme_bw()
print(p_data)

# 1. モデルの構造を定義する関数 (build function) を作成します
# この関数は、推定したいパラメータ(par)を引数に取り、dlmオブジェクトを返します
# 分散は正の値である必要があるため、対数を取った値を推定し、exp()で元に戻します
build_model <- function(par) {
  dlmModReg(
    X = x_t,
    addInt = TRUE, # 切片項を追加します
    dV = exp(par[1]), # 観測ノイズの分散
    dW = exp(par[2:3]) # 状態ノイズの分散（切片項、xの係数）
  )
}

# 2. 最尤推定 (MLE) でパラメータを推定します
# parmにはパラメータの初期値を設定します
mle_fit <- dlmMLE(y = y_t, parm = rep(0, 3), build = build_model)

# 推定結果を確認します
cat("推定の収束状況 (0は正常終了):", mle_fit$conv, "\n")
cat("推定された対数尤度:", -mle_fit$value, "\n\n")

# 推定されたパラメータ
dV_est <- exp(mle_fit$par[1])
dW_est <- exp(mle_fit$par[2:3])
cat("推定された観測ノイズの分散 (dV):", dV_est, "(真の値:", v_var_true, ")\n")
cat("推定されたシステムノイズの分散 (dW for b0):", dW_est[1], "(真の値:", dw0_true, ")\n")
cat("推定されたシステムノイズの分散 (dW for b1):", dW_est[2], "(真の値:", dw1_true, ")\n\n")

# 3. 推定されたパラメータを使って最終的なモデルを構築します
final_model <- build_model(mle_fit$par)

# 4. スムージングを実行して、時変係数を推定します
# dlmSmoothは、全期間のデータを使って各時点の状態を推定します
smoothed_res <- dlmSmooth(y = y_t, mod = final_model)

# 5. 結果を整形します
# スムージングされた状態（係数）の平均値を取得します
# dropFirst=TRUEは、t=0の初期状態を除外するためのオプションです
smooth_means <- dropFirst(smoothed_res$s)
colnames(smooth_means) <- c("b0_est", "b1_est")

# スムージングされた状態の分散共分散行列から、各係数の標準偏差を計算します
smooth_vars <- dlmSvd2var(smoothed_res$U.S, smoothed_res$D.S)

# sapplyの結果(p x T+1)を転置(T+1 x p)し、t=0の初期状態を除外します
# これにより、seは smooth_means と同じ T x p (200 x 2) の行列になります
se <- dropFirst(t(sqrt(sapply(smooth_vars, diag))))
# ======================

# 95%信頼区間を計算します
ci_upper <- smooth_means + 1.96 * se
ci_lower <- smooth_means - 1.96 * se
colnames(ci_upper) <- c("b0_upper", "b1_upper")
colnames(ci_lower) <- c("b0_lower", "b1_lower")

# 全ての推定結果と真の値を一つのデータフレームにまとめます
results_df <- data.frame(
  time = 1:time_len,
  b0_true = b0_true,
  b1_true = b1_true,
  b0_est = smooth_means[, "b0_est"],
  b1_est = smooth_means[, "b1_est"],
  b0_lower = ci_lower[, "b0_lower"],
  b0_upper = ci_upper[, "b0_upper"],
  b1_lower = ci_lower[, "b1_lower"],
  b1_upper = ci_upper[, "b1_upper"]
)

# 6. 結果を可視化します
# 切片係数の比較
p_b0 <- ggplot(results_df, aes(x = time)) +
  geom_ribbon(aes(ymin = b0_lower, ymax = b0_upper), fill = "skyblue", alpha = 0.5) +
  geom_line(aes(y = b0_true, color = "真の値"), size = 1) +
  geom_line(aes(y = b0_est, color = "推定値"), size = 1) +
  scale_color_manual(values = c("真の値" = "red", "推定値" = "blue")) +
  labs(
    title = "時変係数の推定結果: 切片 (β_0,t)",
    x = "時間",
    y = "係数の値",
    color = "系列"
  ) +
  theme_bw() +
  theme(legend.position = "bottom")

# xの係数の比較
p_b1 <- ggplot(results_df, aes(x = time)) +
  geom_ribbon(aes(ymin = b1_lower, ymax = b1_upper), fill = "lightgreen", alpha = 0.5) +
  geom_line(aes(y = b1_true, color = "真の値"), size = 1) +
  geom_line(aes(y = b1_est, color = "推定値"), size = 1) +
  scale_color_manual(values = c("真の値" = "red", "推定値" = "darkgreen")) +
  labs(
    title = "時変係数の推定結果: xの係数 (β_1,t)",
    x = "時間",
    y = "係数の値",
    color = "系列"
  ) +
  theme_bw() +
  theme(legend.position = "bottom")

# プロットを表示します
print(p_b0)
print(p_b1)