Rの関数:dlmModARMA {dlm}

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

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

# シミュレーションの設定
time_len <- 1000
ar_true <- c(0.7, -0.2) # AR(2)の真の係数
# ma_true <- NULL # MA項はなし
sigma2_true <- 1.5 # ノイズの真の分散

# arima.sim を使ってARMA過程に従う時系列データを生成します
y_t <- arima.sim(
  n = time_len,
  model = list(ar = ar_true),
  sd = sqrt(sigma2_true)
)

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

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

# データをプロットして確認します
p_data <- ggplot(sim_data, aes(x = time, y = y)) +
  geom_line(color = "royalblue") +
  labs(
    title = "生成したAR(2)のサンプルデータ",
    x = "時間",
    y = "値"
  ) +
  theme_bw()
print(p_data)

# ヘルパー関数：AR係数を部分自己相関係数(PACF)に変換
# Durbin-Levinsonの再帰アルゴリズムを実装します。
ar_to_pacf <- function(ar) {
  p <- length(ar)
  if (p == 0) {
    return(numeric(0))
  }

  # AR係数を行列として初期化
  phi <- matrix(0, nrow = p, ncol = p)
  phi[1, 1] <- ar[1]

  if (p > 1) {
    for (k in 2:p) {
      # より次数の低い係数を計算
      for (j in 1:(k - 1)) {
        phi[k - 1, j] <- phi[k, j] + phi[k, k] * phi[k, k - j]
      }
    }
  }

  # 対角要素がPACFに対応します
  return(diag(phi))
}


# 1. モデルの構造を定義する関数 (build function) を作成します
build_ar2 <- function(par) {
  ar_coeffs <- ARtransPars(par[1:2])
  sigma2_val <- exp(par[3])
  dlmModARMA(ar = ar_coeffs, sigma2 = sigma2_val)
}

# 2. 最尤推定 (MLE) でパラメータを推定します

# === 初期値の設定 ===
ar_fit_init <- ar.mle(y_t, order.max = 2, demean = FALSE)
ar_coeffs_init <- ar_fit_init$ar
sigma2_init <- ar_fit_init$var.pred

# 上記で定義したヘルパー関数を使い、AR係数をPACFに変換します
pacf_init <- ar_to_pacf(ar_coeffs_init)
# PACFをatanh()で制約のないパラメータに変換します
par_ar_init <- atanh(pacf_init)

par_sigma2_init <- log(sigma2_init)
par_init <- c(par_ar_init, par_sigma2_init)

cat("計算された初期値:", par_init, "\n")

mle_fit <- dlmMLE(
  y = y_t,
  parm = par_init, # 計算した初期値を設定
  build = build_ar2,
  method = "BFGS"
)

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

# 推定されたパラメータを元のスケールに戻します
ar_est <- ARtransPars(mle_fit$par[1:2])
sigma2_est <- exp(mle_fit$par[3])

cat("推定されたAR(1)係数:", ar_est[1], "(真の値:", ar_true[1], ")\n")
cat("推定されたAR(2)係数:", ar_est[2], "(真の値:", ar_true[2], ")\n")
cat("推定されたノイズの分散 (sigma2):", sigma2_est, "(真の値:", sigma2_true, ")\n\n")

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

# 4. カルマンフィルタを実行します
filtered_res <- dlmFilter(y = y_t, mod = final_model)

# 5. 結果を整形し、可視化します
one_step_ahead_forecasts <- filtered_res$f
results_df <- data.frame(
  time = sim_data$time,
  observed = sim_data$y,
  forecast = as.numeric(one_step_ahead_forecasts)
)

p_forecast <- ggplot(results_df, aes(x = time)) +
  geom_line(aes(y = observed, color = "観測値"), alpha = 0.8) +
  geom_line(aes(y = forecast, color = "1期先予測値"), linetype = "dotted", linewidth = 1) +
  scale_color_manual(values = c("観測値" = "black", "1期先予測値" = "red")) +
  labs(
    title = "観測値とフィルタリングによる1期先予測値の比較",
    x = "時間",
    y = "値",
    color = "系列"
  ) +
  theme_bw() +
  theme(legend.position = "bottom")

print(p_forecast)