Rの関数:VAR {vars}

Rの関数から VAR {vars} を確認します。

関数 VAR とは

VAR は、ベクトル自己回帰（Vector Autoregression）モデルを推定するための関数です。

単変量の自己回帰（AR）モデルが「自身の過去の値」のみを用いて現在の値を説明するのに対し、VARモデルは「複数の変数の過去の値」を用いて、変数間の相互依存関係を動的にモデル化します。

この関数は、指定されたラグ次数（過去の期間）に基づいてデータ行列を構築し、方程式ごと（変数ごと）に通常の最小二乗法（OLS）を適用して係数を推定します。

関数 VAR の引数

library(vars)
args(VAR)
function (y, p = 1, type = c("const", "trend", "both", "none"), 
    season = NULL, exogen = NULL, lag.max = NULL, ic = c("AIC", 
        "HQ", "SC", "FPE")) 
NULL
  1. y

    • 分析対象となる多変量時系列データです。
    • データフレームまたは行列（matrix）。各列が変数、各行が観測時点を表します。少なくとも2つの変数が含まれている必要があります。
  2. p

    • VARモデルのラグ次数（Lag Order）です。
    • デフォルト: 1
    • 「何時点過去のデータまでを説明変数として用いるか」を指定します。例えば p=2 ならば、1期前と2期前の値がモデルに含まれます。
    • lag.max が指定されている場合、この引数は無視され（または選択されたラグ数が格納され）ます。
  3. type

    • モデルに含める決定論的項（Deterministic terms）の種類を指定します。
    • 選択肢:

      • "const": 定数項（切片）のみを含めます（デフォルト）。
      • "trend": トレンド項（時間の経過とともに直線的に変化する項）のみを含めます。
      • "both": 定数項とトレンド項の両方を含めます。
      • "none": 決定論的項を含めません。
  4. season

    • 季節ダミー変数を含める場合の、季節の周期（整数）を指定します。
    • 例: 月次データなら 12、四半期データなら 4 を指定します。これにより、季節変動を捉えるための中心化された季節ダミーがモデルに追加されます。
  5. exogen

    • 外生変数（Exogenous variables）のデータ行列です。
    • モデル内の変数には影響を与えるが、モデル内の変数からは影響を受けない（と仮定される）変数を追加する場合に使用します。
  6. lag.max

    • ラグ次数 p を情報量基準に基づいて自動選択する場合の、探索するラグの上限値です。
    • NULL でない場合、内部的に VARselect 関数が呼び出され、指定された情報量基準（ic）を最小化するラグ次数が自動的に p として採用されます。
  7. ic

    • lag.max を使用してラグ次数を自動選択する際に用いる情報量基準です。
    • 選択肢:

      • "AIC": 赤池情報量基準（Akaike Information Criterion）。
      • "HQ": Hannan-Quinn 情報量基準。
      • "SC": Schwarz 基準（BIC: Bayesian Information Criterion）。
      • "FPE": 最終予測誤差（Final Prediction Error）。

シミュレーションコード

以下に、VAR 関数の機能を確認するためのサンプルデータを用いたシミュレーションコードを示します。

このシミュレーションでは、互いに影響し合う2つの変数（「広告費」と「売上」のような関係）を生成します。

  • 変数X: 自身の過去だけでなく、変数Yの過去の影響も受ける。
  • 変数Y: 自身の過去だけでなく、変数Xの過去の影響も受ける。

このような「フィードバック関係」を持つデータを生成し、VAR 関数がその真の係数構造を正しく推定できるかを確認します。

なお、有意水準は5%とします。

シミュレーションデータの生成

# パッケージの読み込み
library(ggplot2)
library(tidyr)

# 乱数シードの固定
seed <- 20260103
set.seed(seed)


# 設定: 2変量VAR(1)モデル
# 変数: y1, y2
# データ長: 200時点
n_obs <- 200

# 係数行列 A の設定 (真のモデル)
# y1(t) = 0.5 * y1(t-1) + 0.3 * y2(t-1) + e1
# y2(t) = -0.4 * y1(t-1) + 0.6 * y2(t-1) + e2
#
# 解説:
# y1 は y2 の過去から正の影響 (0.3) を受ける
# y2 は y1 の過去から負の影響 (-0.4) を受ける
A_true <- matrix(c(0.5, 0.3, -0.4, 0.6), nrow = 2, byrow = TRUE)

# 誤差項の生成 (標準正規分布)
e1 <- rnorm(n_obs)
e2 <- rnorm(n_obs)

# データの格納用行列
data_sim <- matrix(0, nrow = n_obs, ncol = 2)
colnames(data_sim) <- c("Variable_A", "Variable_B")

# ループによる時系列データの生成
# 初期値は0とする
for (t in 2:n_obs) {
  # 前期の値ベクトル
  y_prev <- matrix(data_sim[t - 1, ], ncol = 1)

  # VAR(1) プロセス: y(t) = A %*% y(t-1) + e(t)
  y_curr <- A_true %*% y_prev + c(e1[t], e2[t])

  data_sim[t, ] <- as.vector(y_curr)
}

# 扱いやすいようにデータフレームに変換
df_sim <- as.data.frame(data_sim)
df_sim$Time <- 1:n_obs

cat("データ概要:\n")
cat(sprintf("観測数: %d\n", n_obs))
cat("真のモデル構造:\n")
cat("  Variable_A = 0.5 * A(t-1) + 0.3 * B(t-1) + e\n")
cat("  Variable_B = -0.4 * A(t-1) + 0.6 * B(t-1) + e\n\n")

# データの可視化
# ggplot用にデータを縦持ち（ロング形式）に変換
df_long <- pivot_longer(df_sim,
  cols = c("Variable_A", "Variable_B"),
  names_to = "Variable", values_to = "Value"
)

p1 <- ggplot(df_long, aes(x = Time, y = Value, color = Variable)) +
  geom_line(linewidth = 0.8) +
  labs(
    title = "生成された2変量時系列データ",
    subtitle = "互いに影響し合いながら変動している",
    x = "時間 (Time)",
    y = "値 (Value)",
    color = "変数"
  ) +
  theme_minimal()

print(p1)
データ概要:
観測数: 200
真のモデル構造:
  Variable_A = 0.5 * A(t-1) + 0.3 * B(t-1) + e
  Variable_B = -0.4 * A(t-1) + 0.6 * B(t-1) + e
Figure 1

VARモデルの推定

# VAR関数の実行
# ここでは真のラグ次数 p=1 が分かっているものとして指定します
# 定数項はデータ生成に含めていないため type="const" でも "none" でも良いですが、
# 一般的には定数項を含めて推定することが多いため "const" とします。
# (真の切片は0なので、推定される切片も0に近くなるはずです)

var_model <- VAR(df_sim[, c("Variable_A", "Variable_B")], p = 1, type = "const")

print(summary(var_model))

VAR Estimation Results:
========================= 
Endogenous variables: Variable_A, Variable_B 
Deterministic variables: const 
Sample size: 199 
Log Likelihood: -550.449 
Roots of the characteristic polynomial:
0.6566 0.6566
Call:
VAR(y = df_sim[, c("Variable_A", "Variable_B")], p = 1, type = "const")


Estimation results for equation Variable_A: 
=========================================== 
Variable_A = Variable_A.l1 + Variable_B.l1 + const 

              Estimate Std. Error t value Pr(>|t|)    
Variable_A.l1  0.44906    0.05573   8.058 7.46e-14 ***
Variable_B.l1  0.34416    0.04710   7.307 6.75e-12 ***
const         -0.05708    0.06680  -0.854    0.394    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


Residual standard error: 0.9375 on 196 degrees of freedom
Multiple R-Squared: 0.3948, Adjusted R-squared: 0.3887 
F-statistic: 63.94 on 2 and 196 DF,  p-value: < 2.2e-16 


Estimation results for equation Variable_B: 
=========================================== 
Variable_B = Variable_A.l1 + Variable_B.l1 + const 

              Estimate Std. Error t value Pr(>|t|)    
Variable_A.l1 -0.41945    0.06005  -6.986 4.31e-11 ***
Variable_B.l1  0.63857    0.05075  12.583  < 2e-16 ***
const          0.04334    0.07198   0.602    0.548    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


Residual standard error: 1.01 on 196 degrees of freedom
Multiple R-Squared: 0.4988, Adjusted R-squared: 0.4936 
F-statistic: 97.51 on 2 and 196 DF,  p-value: < 2.2e-16 



Covariance matrix of residuals:
           Variable_A Variable_B
Variable_A    0.87889   -0.06167
Variable_B   -0.06167    1.02032

Correlation matrix of residuals:
           Variable_A Variable_B
Variable_A    1.00000   -0.06512
Variable_B   -0.06512    1.00000
モデル全体の概要 (VAR Estimation Results)
  • Sample size: 199: 元の観測数は200ですが、ラグ次数 \(p=1\) を設定したため、最初の1時点分が初期値として利用され、推定に使用されたサンプルは199個となります。
  • Roots of the characteristic polynomial: 特性多項式の根（0.6566）がすべて1未満（単位円の内側）に収まっています。本結果は、この時系列システムが発散することなく、安定的（定常）であることを示唆しています。
Variable_A の方程式 (Estimation results for equation Variable_A)

ここでは、「Variable_A の現在の値」が「過去の自分」と「過去の Variable_B」からどのような影響を受けているかが推定されています。

  • Variable_A.l1 (1期前の自分):

    • 推定値: 0.44906
    • 真の値: 0.5
    • 解説: 真の値に近い数値が推定されており、p値は設定した有意水準を下回っています。自己相関を正しく捉えているといえます。
  • Variable_B.l1 (1期前の相手):

    • 推定値: 0.34416
    • 真の値: 0.3
    • 解説: こちらも真の値に近く、p値は設定した有意水準を下回っています。「Variable_B の過去が Variable_A にプラスの影響を与える」という因果関係を正しく検出できています。
  • const (定数項):

    • 推定値: -0.05708
    • 解説: p値は 0.394 と設定した有意水準を上回っています。シミュレーションデータは切片ゼロで生成したため、この項がゼロと区別がつかないことは正しい結果です。
Variable_B の方程式 (Estimation results for equation Variable_B)

同様に、「Variable_B の現在の値」に対する影響の推定結果です。

  • Variable_A.l1 (1期前の相手):

    • 推定値: -0.41945
    • 真の値: -0.4
    • 解説: 真の値に近く、設定した有意水準を下回っています。「Variable_A の過去が Variable_B にマイナスの影響を与える」という逆方向の因果関係も特定されています。
  • Variable_B.l1 (1期前の自分):

    • 推定値: 0.63857
    • 真の値: 0.6
    • 解説: こちらも適切に推定されています。
残差の相関 (Correlation matrix of residuals)
  • Variable_A と Variable_B の残差間の相関係数は -0.06512 です。
  • シミュレーションでは互いに独立な（無相関な）誤差項を生成したため、この相関がゼロに近いことは、モデルが構造を適切に説明しており、未説明の部分（残差）に意図しない関係性が残っていないことを示しています。

以上です。