介紹
最近我們被客戶要求撰寫關于向量自回歸的研究報告，包括一些圖形和統計輸出。向量自回歸（VAR）模型的一般缺點是，估計系數的數量與滯后的數量成比例地增加。因此，隨著滯后次數的增加，每個參數可用的信息較少。在貝葉斯VAR文獻中，減輕這種所謂的維數詛咒的一種方法是隨機搜索變量選擇（SSVS）統計分析軟件的選擇，由等人提出（2008）。SSVS的基本思想是將通常使用的先驗方差分配給應包含在模型中的參數，將不相關參數的先驗方差接近零。這樣，通常就可以估算出相關參數，并且無關變量的后驗值接近于零，因此它們對預測和沖激響應沒有顯著影響。這是通過在模型之前添加層次結構來實現的，其中在采樣算法的每個步驟中評估變量的相關性。
視頻：線性回歸中的貝葉斯推斷與R語言預測工人工資數據案例
貝葉斯推斷線性回歸與R語言預測工人工資數據
，時長09:58
視頻：向量自回歸VAR數學原理及R軟件經濟數據脈沖響應分析實例
【視頻】向量自回歸VAR數學原理及R語言軟件經濟數據脈沖響應分析實例
，時長12:01
這篇文章介紹了使用SSVS估計貝葉斯向量自回歸（BVAR）模型。它使用Lü（2007）的數據集E1，其中包含有關至德國固定投資，可支配收入和消費支出的數據。加載數據并生成數據：
# 加載和轉換數據
e1 <- diff(log(e1))
# 生成VAR
data <- gen_var(e1, p = 4, deterministic = "const")
# 獲取數據矩陣
y <- data$Y[, 1:71]
x <- data$Z[, 1:71]估算值
根據等人所述的半自動方法來設置參數的先驗方差（2008）。對于所有變量，先驗包含概率設置為0.5。誤差方差-協方差矩陣的先驗信息不足。
# 重置隨機數提高可重復性
set.seed(1234567)
t <- ncol(y) # 觀察數
k <- nrow(y) # 內生變量數
m <- k * nrow(x) # 估計系數數
# 系數先驗
a_mu_prior <- matrix(0, m) # 先驗均值的向量
# SSVS先驗（半自動方法）
ols <- tcrossprod(y, x) %*% solve(tcrossprod(x)) # OLS估計
sigma_ols <- tcrossprod(y - ols %*% x) / (t - nrow(x)) # OLS誤差協方差矩陣
cov_ols <- kronecker(solve(tcrossprod(x)), sigma_ols)
se_ols <- matrix(sqrt(diag(cov_ols))) # OLS標準誤
# 先驗參數
prob_prior <- matrix(0.5, m)
# 方差-協方差矩陣
u_sigma_df_prior <- 0 # 方差-協方差矩陣
u_sigma_scale_prior <- diag(0, k) # 先驗協方差矩陣
u_sigma_df_post <- t + u_sigma_df_prior # 后驗自由度初始參數值設置為零，這意味著在Gibbs采樣器的第一步中應相對自由地估算所有參數。
可以直接將SSVS添加到VAR模型的標準Gibbs采樣器算法中。在此示例中，常數項從SSVS中排除，這可以通過指定來實現 = 1:36。具有SSVS的Gibbs采樣器的輸出可以用通常的方式進一步分析。因此，可以通過計算參數的繪制方式獲得點估計：
## invest income cons
## invest.1 -0.102 0.011 -0.002
## income.1 0.044 -0.031 0.168
## cons.1 0.074 0.140 -0.287
## invest.2 -0.013 0.002 0.004
## income.2 0.015 0.004 0.315
## cons.2 0.027 -0.001 0.006
## invest.3 0.033 0.000 0.000
## income.3 -0.008 0.021 0.013
## cons.3 -0.043 0.007 0.019
## invest.4 0.250 0.001 -0.005
## income.4 -0.064 -0.010 0.025
## cons.4 -0.023 0.001 0.000
## const 0.014 0.017 0.014
還可以通過計算變量的均值來獲得每個變量的后驗概率。從下面的輸出中可以看出，在VAR（4）模型中似乎只有幾個變量是相關的。常數項的概率為100％，因為它們已從SSVS中排除。
## invest income cons
## invest.1 0.43 0.23 0.10
## income.1 0.10 0.18 0.67
## cons.1 0.11 0.40 0.77
## invest.2 0.11 0.09 0.14
## income.2 0.08 0.07 0.98
## cons.2 0.07 0.06 0.08
## invest.3 0.19 0.07 0.06
## income.3 0.06 0.13 0.10
## cons.3 0.09 0.07 0.12
## invest.4 0.78 0.09 0.16
## income.4 0.13 0.09 0.18
## cons.4 0.09 0.07 0.06
## const 1.00 1.00 1.00
給定這些值，研究人員可以按照常規方式進行操作統計分析軟件的選擇，并根據Gibbs采樣器的輸出獲得預測和脈沖響應。這種方法的優勢在于它不僅考慮了參數不確定性，而且還考慮了模型不確定性。這可以通過系數的直方圖來說明，該直方圖描述了收入的第一個滯后項與消費當前值之間的關系。
hist(draws_a[6,], 通過兩個峰描述模型不確定性，并通過右峰在它們周圍的分布來描述參數不確定性。
但是，如果研究人員不希望使用模型，變量的相關性可能會從采樣算法的一個步驟更改為另一個步驟，那么另一種方法將是僅使用高概率的模型。這可以通過進一步的模擬來完成，在該模擬中，對于不相關的變量使用非常嚴格的先驗，而對于相關參數則使用沒有信息的先驗。
后方抽取的均值類似于Lü（2007，5.2.10節）中的OLS估計值：
## invest income cons
## invest.1 -0.219 0.001 -0.001
## income.1 0.000 0.000 0.262
## cons.1 0.000 0.238 -0.334
## invest.2 0.000 0.000 0.001
## income.2 0.000 0.000 0.329
## cons.2 0.000 0.000 0.000
## invest.3 0.000 0.000 0.000
## income.3 0.000 0.000 0.000
## cons.3 0.000 0.000 0.000
## invest.4 0.328 0.000 -0.001
## income.4 0.000 0.000 0.000
## cons.4 0.000 0.000 0.000
## const 0.015 0.015 0.014
評價
bvar功能可用于將Gibbs采樣器的相關輸出收集到標準化對象中，例如獲得預測或irf進行脈沖響應分析。
hin(bvar_est, thin = 5)
預測
可以使用函數獲得置信區間的預測。
plot(bvar_pred)脈沖響應分析
plot(OIR