R语言建立VAR模型分析联合内生变量的动态关系

最近在写向量自回归的论文，无论是百度还是Google，都没能找到特别合适的R环境下中文资料，大都是Eviews做出来的。所以写这么一篇blog来分享下自己的经验。
注：本文着重介绍VAR的R实现，具体学术性质的东西请参阅相关学术论文。
VAR的定义：
自行Google，很详细，也很简单
VAR模型的用途：
主要是预测分析和内生变量间影响状况分析。
VAR的主要步骤：
（个人拙见，不是标准模板）
    选择合适的变量
    Granger因果检验，进一步观察变量间的关联性，最好做双向检验，不过也有人说单向就足够了，这就人之间人智者见智了
    选择VAR模型滞后阶数
    拟合VAR模型
    诊断性检验：包括系统平稳性检验、正态性检验、序列相关误差等
    脉冲响应分析
    方差分解
    预测分析
各个步骤在R中的实现方法：
R中有个叫“vars”的package，主要用来做向量自回归分析，所以先安装并加载该包：
install.packages(vars)
library(vars)
1.选择变量
根据理论分析选择出相关联的变量，不多说。
2.Granger因果检验
vars包里面有个专门做格兰杰因果检验的函数：

causality(x, cause = NULL, vcov.=NULL, boot=FALSE, boot.runs=100)
另外还有一个适用于普通线性回归模型的Granger test的函数：

grangertest(x, y, order = 1, na.action = na.omit, ...)

这两个函数最直接的区别在于，第二个不用拟合VAR模型即可使用，而第一个必须在拟合VAR模型之后使用。
3.选择合适的滞后阶数
没有一个定论，主要是通过不同信息准则选择出合适的结果，且最好选择最简阶数（也就是最低阶数）。
相关函数：

VARselect(y, lag.max = 10, type = c("const", "trend", "both", "none"),
season = NULL, exogen = NULL)

函数会return一个结果，分别是根据AIC、HQ、SC、FPE四个信息准则得出的最优阶数。
4.拟合VAR模型

var(x, y = NULL, na.rm = FALSE, use)

5.诊断性检验
也就是检验模型的有效性。
系统平稳性：

stability(x, type = c("OLS-CUSUM", "Rec-CUSUM", "Rec-MOSUM",
"OLS-MOSUM", "RE", "ME", "Score-CUSUM", "Score-MOSUM",
"fluctuation"), h = 0.15, dynamic = FALSE, rescale = TRUE)

这里使用“OLS-CUSUM”，它给出的是残差累积和，在该检验生成的曲线图中，残差累积和曲线以时间为横坐标，图中绘出两条临界线，如果累积和超出了这两条临界线，则说明参数不具有稳定性。

结果如下图：

说明系统稳定。
正态性检验：

normality.test(x, multivariate.only = TRUE)

序列相关误差检验：

serial.test(x, lags.pt = 16, lags.bg = 5, type = c("PT.asymptotic",
"PT.adjusted", "BG", "ES") )
6.脉冲响应分析
脉冲响应分析，直白的来说就是对于某一内生变量对于残差冲击的反应。具体而言，他描述的是在随机误差项上施加一个标准差大小的冲击后对内生变量的当期值和未来值所产生的影响。

irf(x, impulse = NULL, response = NULL, n.ahead = 10,
ortho = TRUE, cumulative = FALSE, boot = TRUE, ci = 0.95,
runs = 100, seed = NULL, ...)
示例：

var<-VAR(timeseries,lag.max=2)
var.irf<-irf(var)
plot(var.irf)

结果：

解读：
标题栏说明，这是Y_ln对各个变量（包括Y_ln自身）的脉冲响应（impulse response），其中可以看出来自Y_ln的正向冲击，来自FDI_ln的正向冲击、来自INDUSTRY_ln的冲击不断减小到负向。其余变量的冲击较小。
7.方差分解
VAR模型的应用，还可以采用方差分解方法研究模型的动态特征。方差分解是进一步评价各内生变量对预测方差的贡献度。方差分解是分析预测残差的标准差由不同新息的冲击影响的比例，亦即对应内生变量对标准差的贡献比例。

fevd(x, n.ahead=10, ...)
示例：

var<-VAR(timeseries,lag.max=2)
fevd1<-fevd(var, n.ahead = 5)$Y_ln
结果：
Y_ln REER_ln M0_ln CPI_ln RETAIL_ln FDI_ln INDUSTRY_ln
[1,] 1.0000000 0.000000000 0.0000000 0.00000000 0.00000000 0.00000000 0.00000000
[2,] 0.5660281 0.004363083 0.3085364 0.01686071 0.01356081 0.06509447 0.02555642
[3,] 0.5411924 0.009721985 0.2755711 0.01899613 0.07313395 0.05837871 0.02300568
[4,] 0.5259530 0.020262020 0.2783238 0.01870045 0.06689414 0.06883620 0.02103032
[5,] 0.5268243 0.036825419 0.2697744 0.01855353 0.06276992 0.06550223 0.01975014
解读：
例子中选取的是Y_ln变量的方差分解结果，如果不加‘$Y_ln’，则会return全部变量的结果。
最左边的是滞后期数，一共5期，结果表明当滞后期为1时，其自身对预测方差的贡献率为100%，用人话讲就是自身其变化。随着滞后期增加，Y_ln的贡献率下降，其他变量逐渐增加。不管怎么变化，每一行（也就是每一期）各个变量的贡献率之和都为1。
8.模型预测
没什么好说的，举例示之。

var.predict<-predict(var,n.ahead=3,ci=0.95)
var.predict
结果：
$Y_ln
fcst lower upper CI
[1,] 8.335729 8.208656 8.462802 0.1270727
[2,] 8.284560 8.076325 8.492795 0.2082349
[3,] 8.299723 8.078930 8.520516 0.2207930
fcst:点估计值
lower:区间估计下界
upper:区间估计上界
CI:置信区间
9.预测结果可视化
除了直接使用plot()函数绘图以外，vars包有一个fanchart()函数可以绘制扇形图，示意图：

总结：以上内容基本上实现了建立向量自回归模型，并进行分析所需的主要功能。至于更细分的点，就需要具体问题具体分析了。如文中有任何错误，请及时留言，谢谢