Rでベクトル誤差修正モデル

# シミュレーションのパラメータ設定
seed <- 20250605
set.seed(seed)

n_obs <- 1000 # 生成するデータ点の数

# VECMのパラメータ
# 誤差修正項: ECT_t = Y_t - beta * X_t - c_const
alpha_y <- -0.3 # Yの調整係数 (ECT > 0 のとき Y を減少させる)
alpha_x <- 0.15 # Xの調整係数 (ECT > 0 のとき X を増加させ、βXを増やし差を縮める)
beta <- 0.8 # 共和分ベクトルの係数 (Y_t ≈ beta * X_t + c_const)
c_const <- 10 # 共和分関係の定数項

# 差分ラグ項の係数 (p-1 = 1 の場合)
# ΔY_t = ... + gamma11_1*ΔY_{t-1} + gamma12_1*ΔX_{t-1} + ...
# ΔX_t = ... + gamma21_1*ΔY_{t-1} + gamma22_1*ΔX_{t-1} + ...
gamma11_1 <- 0.4 # ΔY_{t-1} が ΔY_t に与える影響
gamma12_1 <- 0.05 # ΔX_{t-1} が ΔY_t に与える影響
gamma21_1 <- 0.1 # ΔY_{t-1} が ΔX_t に与える影響
gamma22_1 <- 0.2 # ΔX_{t-1} が ΔX_t に与える影響

# 誤差項(イノベーション)の標準偏差
sd_eps_y <- 0.1
sd_eps_x <- 0.1

# 初期値の設定
# YとXの初期値
Y_vec <- numeric(n_obs)
X_vec <- numeric(n_obs)
# Y_vec[1] <- 100
# X_vec[1] <- (Y_vec[1] - c_const) / beta # 均衡点に近い値 X = (100-10)/0.8 = 112.5
Y_vec[1] <- 90 # 少し均衡からずらして開始
X_vec[1] <- 100

# ΔYとΔXの初期値 (t=1時点ではラグがないため0とする)
dY_vec <- numeric(n_obs)
dX_vec <- numeric(n_obs)
dY_vec[1] <- 0
dX_vec[1] <- 0

# 誤差修正項を格納するベクトル
ECT_vec <- numeric(n_obs)

# 時系列データの生成ループ (t = 2 から n_obs まで)
for (t in 2:n_obs) {
  # 1. 前期の誤差修正項 (ECT_{t-1}) を計算
  ECT_current <- Y_vec[t - 1] - beta * X_vec[t - 1] - c_const
  ECT_vec[t - 1] <- ECT_current # 保存用

  # 2. 誤差項 (ε_yt, ε_xt) を生成
  eps_y <- rnorm(1, mean = 0, sd = sd_eps_y)
  eps_x <- rnorm(1, mean = 0, sd = sd_eps_x)

  # 3. ΔY_t と ΔX_t を計算
  dY_vec[t] <- alpha_y * ECT_current +
    gamma11_1 * dY_vec[t - 1] +
    gamma12_1 * dX_vec[t - 1] +
    eps_y

  dX_vec[t] <- alpha_x * ECT_current +
    gamma21_1 * dY_vec[t - 1] +
    gamma22_1 * dX_vec[t - 1] +
    eps_x

  # 4. Y_t と X_t を計算 (レベルデータ)
  Y_vec[t] <- Y_vec[t - 1] + dY_vec[t]
  X_vec[t] <- X_vec[t - 1] + dX_vec[t]
}
# 最後のECTを計算
ECT_vec[n_obs] <- Y_vec[n_obs] - beta * X_vec[n_obs] - c_const


# 生成されたデータをデータフレームにまとめる
sim_data <- tidyr::tibble(
  Time = 1:n_obs,
  Y = Y_vec,
  X = X_vec,
  Delta_Y = dY_vec,
  Delta_X = dX_vec,
  ECT = ECT_vec
)

library(tseries)
adf.test(sim_data$X)
adf.test(sim_data$Y)

adf.test(sim_data$Delta_X)
adf.test(sim_data$Delta_Y)

library(urca)
jo_test <- ca.jo(sim_data[, c("X", "Y")], type = "trace", ecdet = "const", K = 2)
summary(jo_test)

# ect は Error Correction Term(誤差修正項)
vecm_model <- cajorls(jo_test, r = 1)
summary(vecm_model$rlm)

library(dplyr)
# p_level_var = nrow(sim_data) - nrow(vecm_model$rlm$model)
p_level_var <- 2
cat(paste0("VECMの基礎となるレベルVARのラグ次数 (p_level_var): ", p_level_var, "\n"))

# ΔYの当てはめ値を取得
delta_Y_fitted_vecm <- vecm_model$rlm$fitted.values[, "Y.d"]
cat(paste0("VECMによるΔYの当てはめ値の数: ", length(delta_Y_fitted_vecm), "\n"))

# Yのレベルの当てはめ値を計算
# delta_Y_fitted_vecm[j] は、元の時系列の ΔY_{j + p_level_var} の当てはめ。
# Y_fitted_vecm[j] は、元の Y_{j + p_level_var} の当てはめ。
# Y_{j + p_level_var} (fit) = Y_{j + p_level_var - 1} (act) + ΔY_{j + p_level_var} (fit)
Y_fitted_vecm <- numeric(length(delta_Y_fitted_vecm))
for (j in 1:length(delta_Y_fitted_vecm)) {
  Y_fitted_vecm[j] <- sim_data$Y[j + p_level_var - 1] + delta_Y_fitted_vecm[j]
}
cat(paste0("VECMによるYのレベルの当てはめ値の数: ", length(Y_fitted_vecm), "\n"))

# 比較のため、実際のYも同じ期間に合わせる
# VECMの推定には最初の p_level_var 個の観測値がラグとして使われるため、
# 当てはめ/予測は (p_level_var + 1) 番目の観測から始まる。
actual_Y_for_vecm <- sim_data$Y[(p_level_var + 1):nrow(sim_data)]
cat(paste0("VECM比較用の実測Yのデータ数: ", length(actual_Y_for_vecm), "\n"))

# 差分データを準備
diff_data <- sim_data %>%
  mutate(
    Delta_X_actual = c(NA, diff(X)),
    Delta_Y_actual = c(NA, diff(Y))
  ) %>%
  select(Delta_X_actual, Delta_Y_actual) %>%
  na.omit() # 最初のNA行を削除
cat(paste0("差分データの行数: ", nrow(diff_data), "\n"))

# 差分VARのラグ次数 p_diff_var は、VECMの p_level_var - 1 に相当
p_diff_var <- p_level_var - 1
cat(paste0("差分VARモデルのラグ次数 (p_diff_var): ", p_diff_var, "\n"))

# 差分VARモデルを推定
dvar_model <- vars::VAR(diff_data, p = p_diff_var, type = "const")

# 差分VARモデルでΔYを予測
fitted_dvar_output <- fitted(dvar_model)
delta_Y_fitted_dvar <- fitted_dvar_output[, "Delta_Y_actual"]
cat(paste0("差分VARによるΔYの当てはめ値の数: ", length(delta_Y_fitted_dvar), "\n"))

# Yのレベルの当てはめ値を計算
# delta_Y_fitted_dvar[j] は、元の時系列の ΔY_{j + p_diff_var + 1} の当てはめ。
# Y_fitted_dvar[j] は、元の Y_{j + p_diff_var + 1} の当てはめ。
# Y_{j + p_diff_var + 1} (fit) = Y_{j + p_diff_var} (act) + ΔY_{j + p_diff_var + 1} (fit)
Y_fitted_dvar <- numeric(length(delta_Y_fitted_dvar))
for (j in 1:length(delta_Y_fitted_dvar)) {
  Y_fitted_dvar[j] <- sim_data$Y[j + p_diff_var] + delta_Y_fitted_dvar[j]
}
cat(paste0("差分VARによるYのレベルの当てはめ値の数: ", length(Y_fitted_dvar), "\n"))

# 比較のため、実際のYも同じ期間に合わせる
# 差分系列のVAR(p_diff_var)では、
# 1. 差分を取るために最初の1観測が失われる
# 2. VARのラグとして最初のp_diff_var観測が失われる
# よって、当てはめ/予測は元の時系列の (1 + p_diff_var + 1) = (p_diff_var + 2) 番目の観測から始まる。
actual_Y_for_dvar <- sim_data$Y[(p_diff_var + 2):nrow(sim_data)]
cat(paste0("差分VAR比較用の実測Yのデータ数: ", length(actual_Y_for_dvar), "\n"))

cat("\n--- モデル比較 (Yの当てはめ精度) ---\n")

# VECMの評価指標
mse_vecm <- Metrics::mse(actual_Y_for_vecm, Y_fitted_vecm)
mae_vecm <- Metrics::mae(actual_Y_for_vecm, Y_fitted_vecm)
cat("VECM:\n")
cat(paste0("  平均二乗誤差 (MSE): ", round(mse_vecm, 6), "\n"))
cat(paste0("  平均絶対誤差 (MAE): ", round(mae_vecm, 6), "\n"))
cat(paste0("  評価に使用した観測点の数: ", length(actual_Y_for_vecm), "\n"))

# 差分VARの評価指標
mse_dvar <- Metrics::mse(actual_Y_for_dvar, Y_fitted_dvar)
mae_dvar <- Metrics::mae(actual_Y_for_dvar, Y_fitted_dvar)
cat("\n差分VAR (dVAR):\n")
cat(paste0("  平均二乗誤差 (MSE): ", round(mse_dvar, 6), "\n"))
cat(paste0("  平均絶対誤差 (MAE): ", round(mae_dvar, 6), "\n"))
cat(paste0("  評価に使用した観測点の数: ", length(actual_Y_for_dvar), "\n"))