VAR模型分析联合内生变量的动态关系

一、实验介绍

1.1 实验内容

VAR模型是向量自回归模型的简称，是基于数据的统计性质建立的一种常用的计量经济模型，它把系统中每一个内生变量作为系统中所有内生变量的滞后值的函数来构造模型，从而将单变量自回归模型推广到由多元时间序列变量组成的“向量”自回归模型。本实验运用 R 语言来建立两变量的向量自回归模型，首先是检验两变量序列的平稳性，然后进行协整检验，确定待拟合模型的滞后阶数，再拟合VAR模型和做脉冲响应分析，最终对拟合的VAR模型进行分析于预测。通过本实验学会用VAR模型处理多个相关经济指标的分析与预测。

1.2 实验知识点

平稳性检验
协整检验
滞后阶数的确定
VAR 模型的拟合
脉冲响应分析
VAR 模型的预测

1.3 实验环境

R version 3.4.1
Xfce 终端

1.4 适合人群

本课程难度为较难，属于中级级别课程，适合具有 R 语言基础并且有一定经管背景的用户，熟悉 R 语言基础知识并加深巩固，学会用 R 语言建立计量经济模型分析经济数据。

1.5 实验准备

打开Xfce 终端，下载实验所需的数据，并启动 R ：

$ wget http://labfile.oss.aliyuncs.com/courses/910/monthdata.csv
$ sudo R

二、实验原理

以下实验原理介绍来自 MBA 智库百科：向量自回归模型

向量自回归模型简称VAR模型，是一种常用的计量经济模型，1980年由克里斯托弗·西姆斯（Christopher Sims）提出。VAR模型是用模型中所有当期变量对所有变量的若干滞后变量进行回归。VAR模型用来估计联合内生变量的动态关系，而不带有任何事先约束条件。它是AR模型的推广，此模型目前已得到广泛应用。

向量自回归(VAR)是基于数据的统计性质建立模型，VAR模型把系统中每一个内生变量作为系统中所有内生变量的滞后值的函数来构造模型，从而将单变量自回归模型推广到由多元时间序列变量组成的“向量”自回归模型。VAR模型是处理多个相关经济指标的分析与预测最容易操作的模型之一，并且在一定的条件下，多元MA和ARMA模型也可转化成VAR模型，因此近年来VAR模型受到越来越多的经济工作者的重视。

三、实验步骤

本实验利用 2005 年 8 月汇率改革后消费者价格指数（CPI）和人民币名义有效汇率（NEER）的月度数据建立两变量 VAR 模型来分析人民币汇率变动与CPI之间的动态关系。

3.1 实验准备

> data<-read.csv("monthdata.csv")  #读取数据
> CPI<-data[,1]    #将 data 数据的第一列赋值给 CPI 数据框
> NEER<-data[,2]
> length(CPI)  #查看 CPI 变量的数据长度
> length(NEER)

可以看到数据长度为 137，即包含了 137 个月度 CPI 和 NEER 数据。

下面我们绘制两个变量的时序图，观察消费者价格指数和名义汇率的波动情况。

> CPI.ts<-ts(CPI,start=c(2005,8),end=c(2016,12),freq=12)   #将数据转化为时间序列
> NEER.ts<-ts(NEER,start=c(2005,8),end=c(2016,12),freq=12)
> par(mfrow=c(2,1))
> plot(CPI.ts,type="l",xlab="Date", ylab="CPI")
> plot(NEER.ts,type="l",xlab="Date", ylab="NEER")

可以看到消费者价格指数（CPI）的波动幅度较大，而名义有效汇率波动较小但明显呈现一个波动上升的趋势，两个变量序列大致是不平稳的。

3.2 平稳性检验

在拟合 VAR 模型之前，需要对变量进行平稳性检验，如果要拟合的内生变量都是平稳的或者同阶单整的才可进行 VAR 模型的拟合。但在平稳性检验之前需要先项两变量取对数，以消除时间序列的异方差的影响。

> lncpi<-log(CPI)
> lnneer<-log(NEER)

进行平稳性检验本实验选择 urca 包中的 ur.df 函数进行单位根检验，若存在单位根则序列不平稳；反之，不存在单位根则序列平稳。

> install.packages("urca")    #镜像选择 China (Lanzhou) [https]
> library(urca)
> urt.lncpi<-ur.df(lncpi,type='trend',selectlags='AIC')
> urt.lnneer<-ur.df(lnneer,type='trend',selectlags='AIC')
> summary(urt.lncpi)    ############################################### # Augmented Dickey-Fuller Test Unit Root Test # ############################################### Test regression trend Call:lm(formula = z.diff ~ z.lag.1 + 1 + tt + z.diff.lag)Residuals:Min        1Q    Median        3Q       Max -0.028676 -0.003311  0.000381  0.003035  0.043118 Coefficients:Estimate Std. Error t value Pr(>|t|)  (Intercept)  3.047e-01  1.477e-01   2.063   0.0411 *z.lag.1     -6.555e-02  3.185e-02  -2.058   0.0416 *tt          -1.443e-05  1.678e-05  -0.860   0.3912  z.diff.lag  -1.178e-01  8.627e-02  -1.365   0.1745  ---Signif. codes:  0 \u2018***\u2019 0.001 \u2018**\u2019 0.01 \u2018*\u2019 0.05 \u2018.\u2019 0.1 \u2018 \u2019 1Residual standard error: 0.007505 on 131 degrees of freedomMultiple R-squared:  0.05455,    Adjusted R-squared:  0.0329 F-statistic: 2.519 on 3 and 131 DF,  p-value: 0.06084Value of test-statistic is: -2.0578 1.5223 2.2719 Critical values for test statistics: 1pct  5pct 10pcttau3 -3.99 -3.43 -3.13phi2  6.22  4.75  4.07phi3  8.43  6.49  5.47

> summary(urt.lnneer)   ############################################### # Augmented Dickey-Fuller Test Unit Root Test # ############################################### Test regression trend Call:lm(formula = z.diff ~ z.lag.1 + 1 + tt + z.diff.lag)Residuals:Min     1Q    Median        3Q       Max -0.038574 -0.007164 -0.000990  0.009182  0.033482 Coefficients:Estimate Std. Error t value Pr(>|t|)    (Intercept)  0.4429099  0.1380821   3.208  0.00168 ** z.lag.1     -0.0995380  0.0312008  -3.190  0.00178 ** tt           0.0003157  0.0001058   2.983  0.00341 ** z.diff.lag   0.4099249  0.0810914   5.055 1.42e-06 ***---Signif. codes:  0 \u2018***\u2019 0.001 \u2018**\u2019 0.01 \u2018*\u2019 0.05 \u2018.\u2019 0.1 \u2018 \u2019 1Residual standard error: 0.01267 on 131 degrees of freedomMultiple R-squared:  0.1877,    Adjusted R-squared:  0.1691 F-statistic: 10.09 on 3 and 131 DF,  p-value: 5.01e-06Value of test-statistic is: -3.1902 4.2783 5.1509 Critical values for test statistics: 1pct  5pct 10pcttau3 -3.99 -3.43 -3.13phi2  6.22  4.75  4.07phi3  8.43  6.49  5.47

判断是否平稳主要看详细拟合结果的最后两个部分，即：

1.Value of test-statistic is:

2.Critical values for test statistics:

1 是检验统计量的值，2 是对应的显著性水平下检验统计量的临界值。

单位根检验的原假设是序列存在单位根。

对于 lncpi 和 lnneer 的单位根检验结果:

lncpi	lnneer
Value of test-statistic is: -2.0578 1.5223 2.2719	Value of test-statistic is: -3.1902 4.2783 5.1509
Critical values for test statistics:	Critical values for test statistics:
1pct 5pct 10pct	1pct 5pct 10pct
tau3 -3.99 -3.43 -3.13	tau3 -3.99 -3.43 -3.13
phi2 6.22 4.75 4.07	phi2 6.22 4.75 4.07
phi3 8.43 6.49 5.47	phi3 8.43 6.49 5.47

lncpi的检验统计量的值 -2.0578 在 1%、5%、10% 的显著性水平下都大于临界值(-3.99 -3.43 -3.13)，则不能拒绝原价设而接受存在单位根的假设，说明 lncpi 序列是不平稳的。（注意，统计量的值和临界值为负，统计量的值大于临界值是接受原假设；若统计量的值和临界值为正值，统计量的值大于临界值是拒绝原假设。）
同理，lnneer 的检验统计量的值 -3.1902 在 1% 、5% 的显著性水平下都大于对应的临界值(-3.99 -3.43 -3.13)，不能拒绝原假设，lnneer 序列也是不平稳的，存在单位根。

由于两个变量都存在单位根，对序列差分后检验序列的平稳性。

> dlncpi<-diff(lncpi)   #取一阶差分
> dlnneer<-diff(lnneer)
> urt.dlncpi<-ur.df(dlncpi,type='trend',selectlags='AIC')
> urt.dlnneer<-ur.df(dlnneer,type='trend',selectlags='AIC')
> summary(urt.dlncpi)
> summary(urt.dlnneer)

由于检验的详细结果太长，截取需要的部分一阶差分平稳性检验结果：

lncpi 的一阶差分平稳性检验

由于 dlncpi 的检验统计量的值为 -8.1895，在 1%、5%、10% 的显著性水平上都小于对应的临界值，因此拒绝存在单位根的原假设，即 lncpi 是一阶差分平稳的。
lnneer 的一阶差分平稳性检验

由于 dlnneer 的检验统计量的值为 -7.0409，在 1%、5%、10% 的显著性水平上都小于对应的临界值，因此拒绝存在单位根的原假设，即 lncpi 是一阶差分平稳的。

我们再来看一阶差分后的时间序列图：

> par(mfrow=c(2,1))
> plot(dlncpi,type="l",xlab="Date", ylab="diff.CPI")
> plot(dlnneer,type="l",xlab="Date", ylab="diff.NEER")

可以看到差分后的两个变量序列没有明显的波动聚集或则上升下降的趋势，比较平稳。

两个内生变量都是一阶差分平稳的，即都是一阶单整的，是不平稳的时间序列，不能做 Granger 因果检验，只能做协整检验。

3.3 协整检验

协整检验主要针对非平稳的单个序列，但它们的线性组合可能是平稳的。几个变量之间可能存在的一种长期均衡关系进行检验，表现为存在某个协整方程。

由于所有变量都是一阶单整的，是非平稳时间序列，因此各变量之间可能存在协整关系，如果要对所选择的内生变量进行VAR模型的构建，需要进行协整检验，以判断各个变量之间是否存在长期稳定的协整关系，处理各变量之间的是否存在伪回归问题。在本实验我们运用E-G两步法协整检验，需要用到 lmtest 包中的 dwtest 函数。

E-G两步法：第一步:回归方程的估计

> fit<-lm(lncpi~lnneer)
> summary(fit)

由于拟合的模型变量和截距项的 t 检验在 5% 的显著性水平上都是显著的，且模型的 F 检验在5%的显著性水平上也是显著的，因此拟合的线性模型是合理的。

协整回归方程为：lncpi=4.80296-0.03666 lnneer+εt

> install.packages("lmtest")
> library(zoo)
> library(lmtest)
> dwtest(fit)    #检验序列的自相关性

由于 Durbin-Watson test 检验的 p-value 的值几乎为 0，该数值小于 0.05 ，说明在 5% 的显著性水平上残差序列不独立，具有自相关性。

检验残差序列的平稳性：
```
> error<-residuals(fit)  #提取残差序列
> urt.residuals<-ur.df(error,type='none',selectlags='AIC')
> summary(urt.residuals)
```
根据残差序列的平稳性检验结果，在 5% 的显著性水平上拒绝残差序列存在单位根的原假设，即残差序列是平稳的，说明 CPI 和名义有效汇率两个序列之间存在协整关系，意味着我国的消费者价格指数和名义有效汇率之间具有长期均衡关系，增长或者减少具有协同效应。

。
E-G两步法：第二步：误差修正模型的建立

> error.lag<-error[-c(137,138)]
> ecm.fit<-lm(dlncpi~error.lag+dlnneer)   #拟合误差修正模型
> summary(ecm.fit)
> dwtest(ecm.fit)

协整回归方程：

误差修正项的系数为负，符合误差修正机制，反映了上一期偏离长期均衡的数量将在下一期得到反向修正，这也符合之前证明的协整关系。

3.4 滞后阶数的确定

在拟合 VAR 模型之前还需要确定拟合几阶 VAR 模型，也就是确定滞后阶数。确定滞后阶数需要用到 vars 包中的 VARselect 函数。

> install.packages("vars")
> library(MASS)
> library(sandwich)
> library(strucchange)
> library(vars)
> data.new<-data.frame(lncpi,lnneer)  #合并数据
> VARselect(data.new,lag.max=10,type="const")  #在10以内选择最优滞后阶数

根据结果，不同的信息准则有不同的滞后阶数，选择 2 阶或者 4 阶都是可以的，一般来说选择在相同条件下更加简洁的模型，因此选择 2 阶滞后。

3.5 VAR模型的拟合和预测

3.5.1.VAR 模型的拟合

在确定好最优滞后阶数以后我们就可以拟合模型：

> var<-VAR(data.new,lag.max=2,ic="AIC")
> summary(var)VAR Estimation Results:   #VAR(2) 模型拟合结果========================= Endogenous variables: lncpi, lnneer Deterministic variables: const Sample size: 135 Log Likelihood: 873.227 Roots of the characteristic polynomial:0.9621 0.9621 0.3829 0.1847Call:VAR(y = data.new, lag.max = 2, ic = "AIC")Estimation results for equation lncpi:     #模型估计结果一====================================== lncpi = lncpi.l1 + lnneer.l1 + lncpi.l2 + lnneer.l2 + const Estimate Std. Error t value Pr(>|t|)    lncpi.l1   0.80917    0.08578   9.433   <2e-16 ***lnneer.l1 -0.09248    0.04705  -1.965   0.0515 .  lncpi.l2   0.12907    0.08483   1.521   0.1306    lnneer.l2  0.08381    0.04712   1.779   0.0777 .  const      0.32677    0.15855   2.061   0.0413 *  ---Signif. codes:  0 \u2018***\u2019 0.001 \u2018**\u2019 0.01 \u2018*\u2019 0.05 \u2018.\u2019 0.1 \u2018 \u2019 1Residual standard error: 0.007369 on 130 degrees of freedomMultiple R-Squared: 0.8791,    Adjusted R-squared: 0.8754 F-statistic: 236.3 on 4 and 130 DF,  p-value: < 2.2e-16 Estimation results for equation lnneer:        #模型估计结果二======================================= lnneer = lncpi.l1 + lnneer.l1 + lncpi.l2 + lnneer.l2 + const Estimate Std. Error t value Pr(>|t|)    lncpi.l1  -0.18430    0.14954  -1.232 0.220030    lnneer.l1  1.31275    0.08203  16.003  < 2e-16 ***lncpi.l2   0.28850    0.14789   1.951 0.053238 .  lnneer.l2 -0.31998    0.08215  -3.895 0.000156 ***const     -0.44726    0.27641  -1.618 0.108065    ---Signif. codes:  0 \u2018***\u2019 0.001 \u2018**\u2019 0.01 \u2018*\u2019 0.05 \u2018.\u2019 0.1 \u2018 \u2019 1Residual standard error: 0.01285 on 130 degrees of freedomMultiple R-Squared: 0.9908,    Adjusted R-squared: 0.9905 F-statistic:  3504 on 4 and 130 DF,  p-value: < 2.2e-16 Covariance matrix of residuals:       #残差的协方差矩阵lncpi    lnneerlncpi  5.430e-05 7.761e-06lnneer 7.761e-06 1.650e-04Correlation matrix of residuals:      #残差的相关阵lncpi  lnneerlncpi  1.00000 0.08198lnneer 0.08198 1.00000

根据模型拟合结果，所得的 VAR(2) 的模型方程（保留两位小数）为：

也可以用 coef 函数来查看模型估计的简明结果。

> coef(var)   #查看模型估计的简明结果

会拟合出一致的 VAR(2) 方程（系数保留两位小数）。

> plot(var)       #画出每个变量的时序图、残差图、ACF 图、PACF 图
按<Return>键来看下一个图      #也就是<Enter>键

变量 lncpi 的时序图、残差图、ACF 图、PACF 图
变量 lnneer 的时序图、残差图、ACF 图、PACF 图

3.5.2 脉冲响应分析

在查看完拟合结果的图形之后，我们来绘制拟合结果的脉冲响应图，需要用到 var 包中的 irf 函数。然后通过 plot 函数画出图形。

> var.irf<-irf(var)
> plot(var.irf)
按<Return>键来看下一个图：

根据图形 lncpi 自身以及 lnneer 的波动对 lncpi 有正向的冲击。 lncpi 对自身的影响没有滞后期，并且自身波动的影响随着时期的增加会越来越小。lnneer 波动对 lncpi 的影响在第一期以前是逐渐减少的，但在第一期减少为 0 之后随着时间的增加影响越来越大。

lncpi 波动对 lnneer 有一个负向的冲击，并且随着时期数的增加负向的影响会越来越大。lnneer 波动对自身有一个正向的冲击，这个冲击从开始先增加，在第二期达到最大值以后又逐渐减少。

3.5.3 VAR(2)模型的预测

模型预测通常需要用到 predict 函数，具体用法如下：

> var.predict<-predict(var,n.ahead=10,ci=0.95)
> var.predict

这样我们就得到了两变量（lncpi、lnneer）的 VAR(2)模型的滞后 10 期预测结果：

四、实验总结

通过本次实验深入了解向量自回归模型，学会运用 R 语言来进行 VAR模型建模。学习如何检验一个序列的平稳性，以及如何运用协整检验来分析非平稳时间序列的关系。学会运用 var 包中的包括 VARselect() 函数、VAR() 函数、coef() 函数以及 irf() 函数各种函数来进行 VAR 模型的建模，选择最优滞后阶数，拟合 VAR 模型，对模型结果进行脉冲响应分析以及模型的预测。希望同过实验能够独立的运用 R 语言进行 VAR 模型建模于分析，巩固 R 的同时加深自己的金融知识。