Rによるインパルス応答のシミュレーション

# VARモデルの作成 (ラグ次数 p=2)
# -------------------------------------------------
# 関数VARの引数は以下の通り
# -第1引数: 時系列データ (データフレームまたは行列)
# -p: ラグ次数
# -type: 定数項やトレンド項を含めるかどうかを指定
#   "const": 定数項のみ
#   "trend": トレンド項のみ
#   "both": 定数項とトレンド項の両方
#   "none": 定数項もトレンド項もなし
library(vars)
var_model <- VAR(ts_data, p = 2, type = "const")

# VARモデルのラグ次数
p <- var_model$p # p = 2

print("作成したVARモデル:")
var_model

# VARモデルからy1の式の推定係数を取得
# -------------------------------------------------
equation_y1_coeffs <- coef(var_model$varresult$y1)

# 係数を個別の変数に格納
coef_y1_lag1 <- equation_y1_coeffs["y1.l1"]
coef_y1_lag2 <- equation_y1_coeffs["y1.l2"]
coef_y2_lag1 <- equation_y1_coeffs["y2.l1"]
coef_y2_lag2 <- equation_y1_coeffs["y2.l2"]
coef_const_y1 <- equation_y1_coeffs["const"]

# y1の予測値を計算
# -------------------------------------------------
# 予測値を格納するベクトル
y1_predicted <- numeric(n_obs)
y1_predicted[1:p] <- NA # 最初のp時点はラグデータがないためNA

# 元のデータ
y1_actual <- ts_data$y1
y2_actual <- ts_data$y2

# ループで予測値を計算 (t = p+1 から n_obs まで)
for (t in (p + 1):n_obs) {
  y1_predicted[t] <- coef_const_y1 +
    coef_y1_lag1 * y1_actual[t - 1] +
    coef_y1_lag2 * y1_actual[t - 2] +
    coef_y2_lag1 * y2_actual[t - 1] +
    coef_y2_lag2 * y2_actual[t - 2]
}

# y1の残差を計算
# -------------------------------------------------
# 残差を格納するベクトル
residuals_y1 <- numeric(n_obs)
residuals_y1[1:p] <- NA # 最初のp時点はNA

# ループで残差を計算 (t = p+1 から n_obs まで)
for (t in (p + 1):n_obs) {
  residuals_y1[t] <- y1_actual[t] - y1_predicted[t]
}

# NAを除いた残差ベクトル
valid_residuals_y1 <- residuals_y1[!is.na(residuals_y1)]

# y1の残差の標準偏差を計算
# -------------------------------------------------
# 残差標準誤差 = sqrt(RSS / 自由度)
# RSS (Residual Sum of Squares)
rss_y1 <- sum(valid_residuals_y1^2)

# 自由度 (degrees of freedom)
# 有効な観測点の数 = n_obs - p
n_effective_obs <- n_obs - p

# y1の式に含まれる推定パラメータの数 (y1.l1, y1.l2, y2.l1, y2.l2, const の5つ)
num_parameters_y1 <- length(equation_y1_coeffs)

df_residuals_y1 <- n_effective_obs - num_parameters_y1

# y1の残差の分散
variance_residuals_y1 <- rss_y1 / df_residuals_y1

# y1の残差の標準偏差
sd_residuals_y1 <- sqrt(variance_residuals_y1)

print(paste("y1の誤差項の標準偏差:", sd_residuals_y1))

summary_var_model <- summary(var_model)
sd_residuals_y1_vars_pkg <- summary_var_model$varresult$y1$sigma
sd_residuals_y1_vars_pkg_cov <- sqrt(diag(summary_var_model$covres))[1]
print(paste("summaryから得られるy1の残差標準偏差 (直接的に抽出):", sd_residuals_y1_vars_pkg))
print(paste("summaryから得られるy1の残差標準偏差 (残差の分散共分散行列から算出):", sd_residuals_y1_vars_pkg_cov))

# VARモデルからラグ1の係数行列 A_1 を作成
# -----------------------------------------------------------
# A_1 行列は (K x K) の行列で、
# A_1[i,j] は「j番目の変数の1期前の値が、i番目の現在の値に与える影響」を示す係数
# Y_t = A_1 * Y_{t-1} + ...
#
# (y1_t) = (a11  a12) * (y1_{t-1})
# (y2_t)   (a21  a22)   (y2_{t-1})
#
# a11 は y1 の式の y1.l1 の係数
# a12 は y1 の式の y2.l1 の係数
# a21 は y2 の式の y1.l1 の係数
# a22 は y2 の式の y2.l1 の係数

# 変数の数 (K) を取得
K_vars <- var_model$K

# 各方程式の係数を取得
coeffs_eq1 <- coef(var_model$varresult[[1]]) # y1 の方程式の係数
coeffs_eq2 <- coef(var_model$varresult[[2]]) # y2 の方程式の係数

# ts_dataの列名を取得
var_names <- colnames(var_model$y)

# ラグ1の係数名を作成
# 例: "y1.l1", "y2.l1"
lag1_coef_names <- paste0(var_names, ".l1")

# A_1 行列の要素を抽出
# y1 の方程式から
a11 <- coeffs_eq1[lag1_coef_names[1]] # y1.l1 の係数
a12 <- coeffs_eq1[lag1_coef_names[2]] # y2.l1 の係数
# y2 の方程式から
a21 <- coeffs_eq2[lag1_coef_names[1]] # y1.l1 の係数
a22 <- coeffs_eq2[lag1_coef_names[2]] # y2.l1 の係数

# A_1 行列を構成
# A_1 = | a11  a12 |
#       | a21  a22 |
A1_matrix <- matrix(c(a11, a12, a21, a22), nrow = K_vars, ncol = K_vars, byrow = TRUE)

# 誤差項の分散共分散行列
summary_var_model <- summary(var_model)
Sigma_U <- summary_var_model$covres

cat("--- 誤差項の分散共分散行列 Sigma_U ---\n")
print(Sigma_U)
cat("\n")

cat("--- 誤差項の分散共分散行列が対称行列であることを確認 ---\n")
isSymmetric(Sigma_U)
cat("\n")

# Cholesky分解: P P' = Sigma_U (Pは下三角行列)
# Rのchol()は上三角行列R (R'R = Sigma_U) を返すので、P = t(chol(Sigma_U))

cat("--- 固有値をチェックして正定値性(正定値行列)を確認 ---\n")
eigen_values <- eigen(Sigma_U, only.values = TRUE)$values
all(eigen_values > .Machine$double.eps^0.5)
cat("\n")

P_cholesky <- t(chol(Sigma_U))
cat("--- Cholesky分解の下三角行列 P ---\n")
print(P_cholesky)
cat("\n")

# y1への1標準偏差ショックに対する0期目の応答ベクトル
# Pの最初の列ベクトルが、y1へのショックに対する0期目の応答
# (var_modelの変数の順序が y1, y2 である前提)
Y0_response <- P_cholesky[, 1, drop = FALSE] # K_vars x 1 の列行列として保持

cat("--- 0期目応答 (Y0_response) ---\n")
# 行名を設定
var_names <- colnames(var_model$y)
rownames(Y0_response) <- var_names
print(Y0_response)
cat("\n")

# 1期先の応答 (Y_1) の計算
# A1_matrix は (K_vars x K_vars)
# Y0_response は (K_vars x 1)
Y1_response <- A1_matrix %*% Y0_response

cat("--- 1期先の応答 (Y1_response) ---\n")
rownames(Y1_response) <- var_names
print(Y1_response)

irf(var_model,
  impulse = "y1", # y1 の誤差項にショック
  response = c("y1", "y2"), # y1 と y2 の応答を見る
  n.ahead = 1, # 1期先まで計算 (結果は0期目と1期目を含む)
  ortho = TRUE, # 直交化インパルス応答 (ショックは1標準偏差)
  boot = FALSE, # ブートストラップ信頼区間は計算しない
  cumulative = FALSE
) # 各期の応答 (非累積)

# VARモデルからラグ2の係数行列 A_2 を作成
# -----------------------------------------------------------
# A_2 行列は (K x K) の行列で、
# A_2[i,j] は「j番目の変数の2期前の値が、i番目の現在の値に与える影響」を示す係数
# Y_t = ... + A_2 * Y_{t-2} + ...
#
# (y1_t) = ... + (a11_l2  a12_l2) * (y1_{t-2})
# (y2_t)         (a21_l2  a22_l2)   (y2_{t-2})
#
# a11_l2 は y1 の式の y1.l2 の係数
# a12_l2 は y1 の式の y2.l2 の係数
# a21_l2 は y2 の式の y1.l2 の係数
# a22_l2 は y2 の式の y2.l2 の係数

K_vars <- var_model$K
var_names <- colnames(var_model$y)

# 改めて各方程式の係数を取得
coeffs_eq1 <- coef(var_model$varresult[[var_names[1]]]) # y1 の方程式の係数
coeffs_eq2 <- coef(var_model$varresult[[var_names[2]]]) # y2 の方程式の係数

# ラグ2の係数名を作成
lag2_coef_names <- paste0(var_names, ".l2")

# A_2 行列の要素を抽出
# y1 の方程式から
a11_l2 <- coeffs_eq1[lag2_coef_names[1]] # y1.l2 の係数
a12_l2 <- coeffs_eq1[lag2_coef_names[2]] # y2.l2 の係数
# y2 の方程式から
a21_l2 <- coeffs_eq2[lag2_coef_names[1]] # y1.l2 の係数
a22_l2 <- coeffs_eq2[lag2_coef_names[2]] # y2.l2 の係数

# A_2 行列を構成
# A_2 = | a11_l2  a12_l2 |
#       | a21_l2  a22_l2 |
A2_matrix <- matrix(c(a11_l2, a12_l2, a21_l2, a22_l2), nrow = K_vars, ncol = K_vars, byrow = TRUE)

cat("--- 作成した A_2 行列 ---\n")
print(A2_matrix)
cat("\n")

# 2期先の応答 (Y_2) の計算
# -----------------------------------------------------------
# Y_2 = A_1 * Y_1 + A_2 * Y_0
# Y_1 は Y1_response
# Y_0 は Y0_response

# A1_matrix は (K_vars x K_vars)
# Y1_response は (K_vars x 1)
term1_Y2 <- A1_matrix %*% Y1_response

# A2_matrix は (K_vars x K_vars)
# Y0_response は (K_vars x 1)
term2_Y2 <- A2_matrix %*% Y0_response

Y2_response <- term1_Y2 + term2_Y2

cat("--- 2期先の応答 (Y2_response) ---\n")
rownames(Y2_response) <- var_names
print(Y2_response)

irf(var_model,
  impulse = var_names[1], # y1 の誤差項にショック
  response = var_names, # y1 と y2 の応答を見る
  n.ahead = 2, # 2期先まで計算 (結果は0, 1, 2期目を含む)
  ortho = TRUE, # 直交化インパルス応答
  boot = FALSE, # ブートストラップ信頼区間は計算しない
  cumulative = FALSE
)