R語言隨機(jī)搜索變量選擇SSVS估計(jì)貝葉斯向量自回歸（BVAR）模型|附代碼數(shù)據(jù)
原文鏈接：http://tecdat.cn/?p=9390
原文出處：拓端數(shù)據(jù)部落公眾號
?最近我們被客戶要求撰寫關(guān)于貝葉斯向量自回歸（BVAR）的研究報(bào)告，包括一些圖形和統(tǒng)計(jì)輸出。
介紹
向量自回歸（VAR）模型的一般缺點(diǎn)是，估計(jì)系數(shù)的數(shù)量與滯后的數(shù)量成比例地增加。因此，隨著滯后次數(shù)的增加，每個(gè)參數(shù)可用的信息較少。在貝葉斯VAR文獻(xiàn)中，減輕這種所謂的維數(shù)詛咒的一種方法是隨機(jī)搜索變量選擇（SSVS），由George等人提出（2008）。SSVS的基本思想是將通常使用的先驗(yàn)方差分配給應(yīng)包含在模型中的參數(shù)，將不相關(guān)參數(shù)的先驗(yàn)方差接近零。這樣，通常就可以估算出相關(guān)參數(shù)，并且無關(guān)變量的后驗(yàn)值接近于零，因此它們對預(yù)測和沖激響應(yīng)沒有顯著影響。這是通過在模型之前添加層次結(jié)構(gòu)來實(shí)現(xiàn)的，其中在采樣算法的每個(gè)步驟中評估變量的相關(guān)性。
這篇文章介紹了使用SSVS估計(jì)貝葉斯向量自回歸（BVAR）模型。它使用Lütkepohl（2007）的數(shù)據(jù)集E1，其中包含有關(guān)1960Q1至1982Q4德國固定投資，可支配收入和消費(fèi)支出的數(shù)據(jù)。加載數(shù)據(jù)并生成數(shù)據(jù)：
r
# 加載和轉(zhuǎn)換數(shù)據(jù)e1 <- diff(log(e1))# 生成VARdata <- gen_var(e1, p = 4, deterministic = "const")# 獲取數(shù)據(jù)矩陣y <- data$Y[, 1:71]x <- data$Z[, 1:71]
估算值
根據(jù)George等人所述的半自動(dòng)方法來設(shè)置參數(shù)的先驗(yàn)方差（2008）。對于所有變量，先驗(yàn)包含概率設(shè)置為0.5。誤差方差-協(xié)方差矩陣的先驗(yàn)信息不足。
r
# 重置隨機(jī)數(shù)提高可重復(fù)性set.seed(1234567)t <- ncol(y) # 觀察數(shù)k <- nrow(y) # 內(nèi)生變量數(shù)m <- k * nrow(x) # 估計(jì)系數(shù)數(shù)# 系數(shù)先驗(yàn)a_mu_prior <- matrix(0, m) # 先驗(yàn)均值的向量# SSVS先驗(yàn)（半自動(dòng)方法）ols <- tcrossprod(y, x) %*% solve(tcrossprod(x)) # OLS估計(jì)sigma_ols <- tcrossprod(y - ols %*% x) / (t - nrow(x)) # OLS誤差協(xié)方差矩陣cov_ols <- kronecker(solve(tcrossprod(x)), sigma_ols)se_ols <- matrix(sqrt(diag(cov_ols))) # OLS標(biāo)準(zhǔn)誤 # 先驗(yàn)參數(shù)prob_prior <- matrix(0.5, m)# ?方差-協(xié)方差矩陣u_sigma_df_prior <- 0 # 方差-協(xié)方差矩陣u_sigma_scale_prior <- diag(0, k) # 先驗(yàn)協(xié)方差矩陣u_sigma_df_post <- t + u_sigma_df_prior # 后驗(yàn)自由度
初始參數(shù)值設(shè)置為零，這意味著在Gibbs采樣器的第一步中應(yīng)相對自由地估算所有參數(shù)。
可以直接將SSVS添加到VAR模型的標(biāo)準(zhǔn)Gibbs采樣器算法中。在此示例中，常數(shù)項(xiàng)從SSVS中排除，這可以通過指定來實(shí)現(xiàn)include = 1:36。具有SSVS的Gibbs采樣器的輸出可以用通常的方式進(jìn)一步分析。因此，可以通過計(jì)算參數(shù)的繪制方式獲得點(diǎn)估計(jì)：
## ? ? ? ? ?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
還可以通過計(jì)算變量的均值來獲得每個(gè)變量的后驗(yàn)概率。從下面的輸出中可以看出，在VAR（4）模型中似乎只有幾個(gè)變量是相關(guān)的。常數(shù)項(xiàng)的概率為100％，因?yàn)樗鼈円褟腟SVS中排除。
## ? ? ? ? ?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
給定這些值，研究人員可以按照常規(guī)方式進(jìn)行操作，并根據(jù)Gibbs采樣器的輸出獲得預(yù)測和脈沖響應(yīng)。這種方法的優(yōu)勢在于它不僅考慮了參數(shù)不確定性，而且還考慮了模型不確定性。這可以通過系數(shù)的直方圖來說明，該直方圖描述了收入的第一個(gè)滯后項(xiàng)與消費(fèi)當(dāng)前值之間的關(guān)系。
r
hist(draws_a[6,],
?
通過兩個(gè)峰描述模型不確定性，并通過右峰在它們周圍的分布來描述參數(shù)不確定性。
但是，如果研究人員不希望使用模型，變量的相關(guān)性可能會(huì)從采樣算法的一個(gè)步驟更改為另一個(gè)步驟，那么另一種方法將是僅使用高概率的模型。這可以通過進(jìn)一步的模擬來完成，在該模擬中，對于不相關(guān)的變量使用非常嚴(yán)格的先驗(yàn)，而對于相關(guān)參數(shù)則使用沒有信息的先驗(yàn)。
后方抽取的均值類似于Lütkepohl（2007，5.2.10節(jié)）中的OLS估計(jì)值：
## ? ? ? ? ?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
評價(jià)
bvar功能可用于將Gibbs采樣器的相關(guān)輸出收集到標(biāo)準(zhǔn)化對象中，例如predict獲得預(yù)測或irf進(jìn)行脈沖響應(yīng)分析。
?
r
hin(bvar_est, thin = 5)
預(yù)測
可以使用函數(shù)獲得置信區(qū)間的預(yù)測predict。
r
plot(bvar_pred)
?
脈沖響應(yīng)分析
r
plot(OIR
最受歡迎的見解
1.matlab使用貝葉斯優(yōu)化的深度學(xué)習(xí)
2.matlab貝葉斯隱馬爾可夫hmm模型實(shí)現(xiàn)
3.R語言Gibbs抽樣的貝葉斯簡單線性回歸仿真
4.R語言中的block Gibbs吉布斯采樣貝葉斯多元線性回歸
5.R語言中的Stan概率編程MCMC采樣的貝葉斯模型
6.Python用PyMC3實(shí)現(xiàn)貝葉斯線性回歸模型
7.R語言使用貝葉斯 層次模型進(jìn)行空間數(shù)據(jù)分析
8.R語言隨機(jī)搜索變量選擇SSVS估計(jì)貝葉斯向量自回歸（BVAR）模型
9.matlab貝葉斯隱馬爾可夫hmm模型實(shí)現(xiàn)
