Rでインパルス応答における変数順序の重要性確認

library(vars)
library(MASS)
library(tidyverse)
library(patchwork)

seed <- 20250625
set.seed(seed)

# データ生成のパラメータ設定
n_obs <- 200 # サンプルサイズ
# 誤差項の共分散行列（非対角成分が0でないことが重要）
sigma <- matrix(c(
  1.0, 0.6,
  0.6, 1.0
), nrow = 2, ncol = 2)

# 2変量正規分布に従う誤差項を生成
errors <- mvrnorm(n = n_obs, mu = c(0, 0), Sigma = sigma)
u1 <- errors[, 1]
u2 <- errors[, 2]

# 2変数の時系列データをVAR(1)プロセスで生成
y1 <- numeric(n_obs)
y2 <- numeric(n_obs)

for (t in 2:n_obs) {
  y1[t] <- 0.7 * y1[t - 1] + 0.2 * y2[t - 1] + u1[t]
  y2[t] <- 0.1 * y1[t - 1] + 0.5 * y2[t - 1] + u2[t]
}

# データフレームにまとめる
common_data <- data.frame(y1 = y1, y2 = y2)

# ケース1: 順序 (y1, y2)
data_order1 <- common_data[, c("y1", "y2")]
var_model1 <- VAR(data_order1, p = 1, type = "none")

# ケース2: 順序 (y2, y1)
data_order2 <- common_data[, c("y2", "y1")]
var_model2 <- VAR(data_order2, p = 1, type = "none")

# OIRFの計算 (ortho = TRUE)
oirf1 <- irf(var_model1, impulse = "y1", response = "y2", n.ahead = 20, ortho = TRUE, boot = FALSE)
oirf2 <- irf(var_model2, impulse = "y1", response = "y2", n.ahead = 20, ortho = TRUE, boot = FALSE)

# 結果のデータフレームを作成して比較
oirf_results <- data.frame(
  period = 1:21,
  order_y1_y2 = oirf1$irf$y1[, 1],
  order_y2_y1 = oirf2$irf$y1[, 1]
)

cat("--- 直交化インパルス応答(OIRF)の結果比較 ---\n")
print(head(oirf_results))

# 2つの結果が異なることを確認
cat("\n2つの結果は同じか？ ->", all.equal(oirf_results$order_y1_y2, oirf_results$order_y2_y1))

# GIRFの計算 (ortho = FALSE)
girf1 <- irf(var_model1, impulse = "y1", response = "y2", n.ahead = 20, ortho = FALSE, boot = FALSE)
girf2 <- irf(var_model2, impulse = "y1", response = "y2", n.ahead = 20, ortho = FALSE, boot = FALSE)

# 結果のデータフレームを作成して比較
girf_results <- data.frame(
  period = 1:21,
  order_y1_y2 = girf1$irf$y1[, 1],
  order_y2_y1 = girf2$irf$y1[, 1]
)

cat("\n--- 一般化インパルス応答(GIRF)の結果比較 ---\n")
print(head(girf_results))

# 2つの結果が同じであることを確認
cat("\n2つの結果は同じか？ ->", all.equal(girf_results$order_y1_y2, girf_results$order_y2_y1))

# OIRFとGIRFの全結果を一つのデータフレームにまとめる
# 信頼区間も取得するためにboot = TRUEで再計算
oirf1_plot <- irf(var_model1, n.ahead = 20, ortho = TRUE, boot = TRUE, seed = 123)
oirf2_plot <- irf(var_model2, n.ahead = 20, ortho = TRUE, boot = TRUE, seed = 123)
girf1_plot <- irf(var_model1, n.ahead = 20, ortho = FALSE, boot = TRUE, seed = 123)
girf2_plot <- irf(var_model2, n.ahead = 20, ortho = FALSE, boot = TRUE, seed = 123)

# プロット用の関数
create_plot_df <- function(irf_obj, order_name) {
  plot_df <- tibble()
  for (imp in names(irf_obj$irf)) {
    for (res in colnames(irf_obj$irf[[imp]])) {
      df <- tibble(
        impulse = imp,
        response = res,
        period = 0:20,
        value = irf_obj$irf[[imp]][, res],
        lower = irf_obj$Lower[[imp]][, res],
        upper = irf_obj$Upper[[imp]][, res],
        order = order_name
      )
      plot_df <- bind_rows(plot_df, df)
    }
  }
  return(plot_df)
}

# データフレームの作成
oirf_plot_data <- bind_rows(
  create_plot_df(oirf1_plot, "Order: y1, y2"),
  create_plot_df(oirf2_plot, "Order: y2, y1")
)

girf_plot_data <- bind_rows(
  create_plot_df(girf1_plot, "Order: y1, y2"),
  create_plot_df(girf2_plot, "Order: y2, y1")
)