这里分享一下R语言实现VAR和SVAR的整个流程。
主要步骤包括：
1.单位根检验
2.确定滞后阶数
3.格兰杰因果检验
4.模型稳定性检验
5.脉冲响应
6.方差分解
(Johansen协整检验，如果需要的话)
整个过程用到的R语言的扩展包有：
library(zoo)
library(vars)
library(tseries)
首先，数据是下面的样子：
ps：数据是时间序列类型,可以通过下面方法将dataframe转成时间序列类型
data = ts(data)
1.单位根检验
#对data的第一列进行单位根检验
adf.test(data[,1])
2.滞后阶数确定
VARselect函数结果包括AIC、HQ、SC和FPE准则
#参数y为时间序列数据，lag.max为最大滞后阶数
#参数type值包括const截距，trend趋势，both同时包含截距和趋势，none不包含截距和趋势
VARselect(y=data, lag.max = 10, type = c("const"))
3.格兰杰因果检验
格兰杰因果检验有两个方法，第一个是在构造模型之前，第二个是在构造模型之后在模型的基础上进行格兰杰因果检验。
(1)构造模型之前格兰杰因果检验
#函数格式：grangertest(yt~xt)
eg：
grangertest(Value~BCI)
(2)构造模型之后格兰杰因果检验
#函数格式：causality(VARModel,cause)
eg
var = VAR(data ,p = 2, type = "const")
causality(var,cause=c('Count','Value'))
ps:在这里如果想要构建SVAR模型的话，需要根据实际情况构建两个矩阵amat和bmat，然后使用这两个矩阵来构建SVAR模型：
svar = SVAR(var,Amat = amat,Bmat = bmat)
4.模型稳定性检验
#这里使用“OLS-CUSUM”，它给出的是残差累积和，在该检验生成的曲线图中，残差累积和曲线以时间为横坐标，
#图中绘出两条临界线，如果累积和超出了这两条临界线，则说明参数不具有稳定性。
sta = stability(var, type = c("OLS-CUSUM"), h = 0.15, dynamic = FALSE, rescale = TRUE)
plot(sta)##结果稳健
5.脉冲响应
#标题栏说明，这是BCI(或者其他变量)对各个变量（包括BCI自身）的脉冲响应
(1)VAR脉冲响应
var.irf<-irf(var,n.head=10)
plot(var.irf)
(2)SVAR脉冲响应
svar.irf<-irf(svar,n.ahead = 100)
plot(svar.irf)
6.方差分解
#反映了各变量的贡献率
(1)VAR方差分解
fevd1<-fevd(var, n.ahead = 10)
fevd1$Count
(2)SVAR方差分解
fevd2<-fevd(svar, n.ahead = 10)
fevd2$Value
ps：有时候需要进行Johansen协整检验
#Johansen协整检验，
#对r=0（不存在协整关系）的检验统计量大于临界值，表明拒绝原假设
yJoTest = ca.jo(data, type = c("trace"), ecdet = c("none"), K = 2)
summary(yJoTest)
网页链接
lag参数是指的滞后项长度的最大可能值；type参数的设置与DF检验一致，根据是否包含截距项和时间趋势项而定；
由于在进行检验前是不知道真实滞后项长度的，需要使用selectlags参数选择确定“最佳”滞后项的原则，“最佳”滞后长度不大于lag参数所对应的值。
使用R语言进行协整关系检验协整检验是为了检验非平稳序列的因果关系，协整检验是解决伪回归为问题的重要方法。首先回归伪回归例子：伪回归Spurious regression伪回归方程的拟合优度、显著性水平等指标都很好，但是其残差序列是一个非平稳序列，拟合一个伪回归:#调用相关R包library(lmtest)library(tseries)#模拟序列set.seed(123456)e1=rnorm(500)e2=rnorm(500)trd=1:500y1=0.8*trd+cumsum(e1)y2=0.6*trd+cumsum(e2)sr.reg=lm(y1~y2)#提取回归残差error=residuals(sr.reg)#作残差散点图plot(error, main="Plot of error")#对残差进行单位根检验adf.test(error)## Dickey-Fuller = -2.548, Lag order = 7, p-value = 0.3463## alternative hypothesis: stationary#伪回归结果，相关参数都显著summary(sr.reg)## Residuals:## Min 1Q Median 3Q Max## -30.654 -11.526 0.359 11.142 31.006## Coefficients:## Estimate Std. Error t value Pr(>|t|)## (Intercept) -29.32697 1.36716 -21.4 <2e-16 ***## y2 1.44079 0.00752 191.6 <2e-16 ***## Residual standard error: 13.7 on 498 degrees of freedom## Multiple R-squared: 0.987, Adjusted R-squared: 0.987## F-statistic: 3.67e+04 on 1 and 498 DF, p-value: <2e-16dwtest(sr.reg)## DW = 0.0172, p-value <2.2e-16恩格尔-格兰杰检验Engle-Granger第一步：建立两变量（y1,y2）的回归方程，第二部：对该回归方程的残差(resid)进行单位根检验其中，原假设两变量不存在协整关系，备择假设是两变量存在协整关系。利用最小二乘法对回归方程进行估计，从回归方程中提取残差进行检验。set.seed(123456)e1=rnorm(100)e2=rnorm(100)y1=cumsum(e1)y2=0.6*y1+e2# (伪)回归模型lr.reg=lm(y2~y1)error=residuals(lr.reg)adf.test(error)## Dickey-Fuller = -3.988, Lag order = 4, p-value = 0.01262## alternative hypothesis: stationaryerror.lagged=error[-c(99,100)]#建立误差修正模型ECM.REGdy1=diff(y1)dy2=diff(y2)diff.dat=data.frame(embed(cbind(dy1, dy2),2))#emed表示嵌入时间序列dy1,dy2到diff.datcolnames(diff.dat)=c("dy1","dy2","dy1.1","dy2.1")ecm.reg=lm(dy2~error.lagged+dy1.1+dy2.1, data=diff.dat)summary(ecm.reg)## Residuals:## Min 1Q Median 3Q Max## -2.959 -0.544 0.137 0.711 2.307## Coefficients:## Estimate Std. Error t value Pr(>|t|)## (Intercept) 0.0034 0.1036 0.03 0.97## error.lagged -0.9688 0.1585 -6.11 2.2e-08 ***## dy1.1 0.8086 0.1120 7.22 1.4e-10 ***## dy2.1 -1.0589 0.1084 -9.77 5.6e-16 ***## Residual standard error: 1.03 on 94 degrees of freedom## Multiple R-squared: 0.546, Adjusted R-squared: 0.532## F-statistic: 37.7 on 3 and 94 DF, p-value: 4.24e-16par(mfrow=c(2,2))plot(ecm.reg)Johansen-Juselius(JJ)协整检验法，该方法是一种用向量自回归(VAR)模型进行检验的方法，适用于对多重一阶单整I(1)序列进行协整检验。JJ检验有两种：特征值轨迹检验和最大特征值检验。我们可以调用urca包中的ca.jo命令完成这两种检验。其语法：ca.jo(x, type = c("eigen", "trace"), ecdet = c("none", "const", "trend"), K = 2,spec=c("longrun", "transitory"), season = NULL, dumvar = NULL)其中：x为矩阵形式数据框；type用来设置检验方法；ecdet用于设置模型形式：none表示不带截距项，const表示带常数截距项，trend表示带趋势项。K表示自回归序列的滞后阶数；spec表示向量误差修正模型反映的序列间的长期或短期关系；season表示季节效应；dumvar表示哑变量设置。set.seed(12345)e1=rnorm(250,0,0.5)e2=rnorm(250,0,0.5)e3=rnorm(250,0,0.5)#模拟没有移动平均的向量自回归序列；u1.ar1=arima.sim(model=list(ar=0.75), innov=e1, n=250)u2.ar1=arima.sim(model=list(ar=0.3), innov=e2, n=250)y3=cumsum(e3)y1=0.8*y3+u1.ar1y2=-0.3*y3+u2.ar1#合并y1,y2,y3构成进行JJ检验的数据库；y.mat=data.frame(y1, y2, y3)#调用urca包中cajo命令对向量自回归序列进行JJ协整检验vecm=ca.jo(y.mat)jo.results=summary(vecm)#cajorls命令可以得到限制协整阶数的向量误差修正模型的最小二乘法回归结果vecm.r2=cajorls(vecm, r=2)vecm.r2## Call:lm(formula = substitute(form1), data = data.mat)## Coefficients:## y1.d y2.d y3.d## ect1 -0.33129 0.06461 0.01268## ect2 0.09447 -0.70938 -0.00916## constant 0.16837 -0.02702 0.02526## y1.dl1-0.22768 0.02701 0.06816## y2.dl1 0.14445 -0.71561 0.04049## y3.dl1 0.12347 -0.29083 -0.07525## $beta## ect1 ect2## y1.l2 1.000e+00 0.0000## y2.l2 -3.402e-18 1.0000## y3.l2 -7.329e-01 0.2952
