21向量自回归模型
向量自回归模型简称VAR模型,是一种常用的计量经济模型,1980年由克里斯托弗·西姆斯(Christopher Sims)提出。VAR模型是用模型中所有当期变量对所有变量的若干滞后变量进行回归。VAR模型用来估计联合内生变量的动态关系,而不带有任何事先约束条件。它是AR模型的推广,此模型目前已得到广泛应用。 向量自回归(VAR)是基于数据的统计性质建立模型,VAR模型把系统中每一个内生变量作为系统中所有内生变量的滞后值的函数来构造模型,从而将单变量自回归模型推广到由多元时间序列变量组成的“向量”自回归模型。VAR模型是处理多个相关经济指标的分析与预测最容易操作的模型之一,并且在一定的条件下,多元MA和ARMA模型也可转化成VAR模型,因此近年来VAR模型受到越来越多的经济工作者的重视。

利用 2005 年 8 月汇率改革后消费者价格指数(CPI)和人民币名义有效汇率(NEER)的月度数据建立两变量 VAR 模型来分析人民币汇率变动与CPI之间的动态关系。
实验准备
> data <- readLines("http://labfile.oss.aliyuncs.com/courses/910/monthdata.csv")
> data <- unlist(strsplit(data, split=","))
> data <- as.data.frame(matrix(data, ncol=2, byrow=T))
> colnames(data) <- data[1,]
> data <- data[-1,]
> data <- as.data.frame(sapply(data, as.numeric))
> CPI <- data[,1]
> NEER <- data[,2]
> length(CPI)
[1] 137
> length(NEER)
[1] 137
可以看到数据长度为 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)的波动幅度较大,而名义有效汇率波动较小但明显呈现一个波动上升的趋势,两个变量序列大致是不平稳的。
平稳性检验
在拟合 VAR 模型之前,需要对变量进行平稳性检验,如果要拟合的内生变量都是平稳的或者同阶单整的才可进行 VAR 模型的拟合。但在平稳性检验之前需要先项两变量取对数,以消除时间序列的异方差的影响。
> lncpi <- log(CPI)
> lnneer <- log(NEER)
进行平稳性检验本实验选择 urca 包中的 ur.df 函数进行单位根检验,若存在单位根则序列不平稳;反之,不存在单位根则序列平稳。
> 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 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.007505 on 131 degrees of freedom
Multiple R-squared: 0.05455, Adjusted R-squared: 0.0329
F-statistic: 2.519 on 3 and 131 DF, p-value: 0.06084
Value of test-statistic is: -2.0578 1.5223 2.2719
Critical values for test statistics:
1pct 5pct 10pct
tau3 -3.99 -3.43 -3.13
phi2 6.22 4.75 4.07
phi3 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 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.01267 on 131 degrees of freedom
Multiple R-squared: 0.1877, Adjusted R-squared: 0.1691
F-statistic: 10.09 on 3 and 131 DF, p-value: 5.01e-06
Value of test-statistic is: -3.1902 4.2783 5.1509
Critical values for test statistics:
1pct 5pct 10pct
tau3 -3.99 -3.43 -3.13
phi2 6.22 4.75 4.07
phi3 8.43 6.49 5.47
判断是否平稳主要看详细拟合结果的最后两个部分,即:
1.Value of test-statistic is:
2.Critical values for test statistics:
1 是检验统计量的值,2 是对应的显著性水平下检验统计量的临界值。 单位根检验的原假设是序列存在单位根。
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)
###############################################
# 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.029902 -0.003086 0.000184 0.003073 0.043732
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.399e-04 1.345e-03 0.476 0.635
z.lag.1 -1.086e+00 1.327e-01 -8.189 2.13e-13 ***
tt -8.213e-06 1.710e-05 -0.480 0.632
z.diff.lag -5.347e-02 8.749e-02 -0.611 0.542
---
Signif. codes:
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.007642 on 130 degrees of freedom
Multiple R-squared: 0.5753, Adjusted R-squared: 0.5655
F-statistic: 58.7 on 3 and 130 DF, p-value: < 2.2e-16
Value of test-statistic is: -8.1895 22.3559 33.5337
Critical values for test statistics:
1pct 5pct 10pct
tau3 -3.99 -3.43 -3.13
phi2 6.22 4.75 4.07
phi3 8.43 6.49 5.47
> summary(urt.dlnneer)
###############################################
# 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.038709 -0.008143 0.000447 0.008358 0.032081
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.369e-03 2.342e-03 1.011 0.314
z.lag.1 -7.016e-01 9.964e-02 -7.041 9.91e-11 ***
tt -7.819e-06 2.941e-05 -0.266 0.791
z.diff.lag 8.103e-02 8.737e-02 0.927 0.355
---
Signif. codes:
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.01314 on 130 degrees of freedom
Multiple R-squared: 0.3291, Adjusted R-squared: 0.3136
F-statistic: 21.26 on 3 and 130 DF, p-value: 2.884e-11
Value of test-statistic is: -7.0409 16.5395 24.8091
Critical values for test statistics:
1pct 5pct 10pct
tau3 -3.99 -3.43 -3.13
phi2 6.22 4.75 4.07
phi3 8.43 6.49 5.47
由于 dlncpi 的检验统计量的值为 -8.1895,在 1%、5%、10% 的显著性水平上都小于对应的临界值,因此拒绝存在单位根的原假设,即 lncpi 是一阶差分平稳的。
由于 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 因果检验,只能做协整检验。
协整检验
协整检验主要针对非平稳的单个序列,但它们的线性组合可能是平稳的。几个变量之间可能存在的一种长期均衡关系进行检验,表现为存在某个协整方程。
由于所有变量都是一阶单整的,是非平稳时间序列,因此各变量之间可能存在协整关系,如果要对所选择的内生变量进行VAR模型的构建,需要进行协整检验,以判断各个变量之间是否存在长期稳定的协整关系,处理各变量之间的是否存在伪回归问题。在本实验我们运用E-G两步法协整检验。
第一步:回归方程的估计
> fit <- lm(lncpi~lnneer)
> summary(fit)
Call:
lm(formula = lncpi ~ lnneer)
Residuals:
Min 1Q Median 3Q Max
-0.047468 -0.008233 -0.001576 0.009218 0.052499
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.80296 0.06073 79.085 <2e-16 ***
lnneer -0.03666 0.01305 -2.809 0.0057 **
---
Signif. codes:
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.02031 on 135 degrees of freedom
Multiple R-squared: 0.05523, Adjusted R-squared: 0.04823
F-statistic: 7.891 on 1 and 135 DF, p-value: 0.005705
由于拟合的模型变量和截距项的 t 检验在 5% 的显著性水平上都是显著的,且模型的 F 检验在5%的显著性水平上也是显著的,因此拟合的线性模型是合理的。 协整回归方程为:lncpi=4.80296-0.03666 lnneer+εt
> library(zoo)
> library(lmtest)
> dwtest(fit)
Durbin-Watson test
data: fit
DW = 0.14135, p-value < 2.2e-16
alternative hypothesis: true autocorrelation is greater than 0
由于 Durbin-Watson test 检验的 p-value 的值几乎为 0,该数值小于 0.05 ,说明在 5% 的显著性水平上残差序列不独立,具有自相关性。
检验残差序列的平稳性:
> error <- residuals(fit)
> urt.residuals <- ur.df(error, type='none', selectlags='AIC')
> summary(urt.residuals)
###############################################
# Augmented Dickey-Fuller Test Unit Root Test #
###############################################
Test regression none
Call:
lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)
Residuals:
Min 1Q Median 3Q Max
-0.027777 -0.003141 0.000325 0.003295 0.043617
Coefficients:
Estimate Std. Error t value Pr(>|t|)
z.lag.1 -0.06830 0.03224 -2.118 0.036 *
z.diff.lag -0.13308 0.08540 -1.558 0.122
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.007454 on 133 degrees of freedom
Multiple R-squared: 0.05906, Adjusted R-squared: 0.04491
F-statistic: 4.174 on 2 and 133 DF, p-value: 0.01745
Value of test-statistic is: -2.1184
Critical values for test statistics:
1pct 5pct 10pct
tau1 -2.58 -1.95 -1.62
根据残差序列的平稳性检验结果,在 5% 的显著性水平上拒绝残差序列存在单位根的原假设,即残差序列是平稳的,说明 CPI 和名义有效汇率两个序列之间存在协整关系,意味着我国的消费者价格指数和名义有效汇率之间具有长期均衡关系,增长或者减少具有协同效应。
第二步:误差修正模型的建立
> error.lag <- error[-c(137, 138)]
> ecm.fit <- lm(dlncpi~error.lag+dlnneer)
> summary(ecm.fit)
Call:
lm(formula = dlncpi ~ error.lag + dlnneer)
Residuals:
Min 1Q Median 3Q Max
-0.033309 -0.003180 0.000015 0.003477 0.042589
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -2.299e-05 6.541e-04 -0.035 0.9720
error.lag -8.414e-02 3.236e-02 -2.600 0.0104 *
dlnneer 3.008e-02 4.744e-02 0.634 0.5271
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.00748 on 133 degrees of freedom
Multiple R-squared: 0.04846, Adjusted R-squared: 0.03415
F-statistic: 3.387 on 2 and 133 DF, p-value: 0.03676
> dwtest(ecm.fit)
Durbin-Watson test
data: ecm.fit
DW = 2.1784, p-value = 0.8259
alternative hypothesis: true autocorrelation is greater than 0
协整回归方程:

误差修正项的系数为负,符合误差修正机制,反映了上一期偏离长期均衡的数量将在下一期得到反向修正,这也符合之前证明的协整关系。
滞后阶数的确立
在拟合 VAR 模型之前还需要确定拟合几阶 VAR 模型,也就是确定滞后阶数。
> library(MASS)
> library(sandwich)
> library(strucchange)
> library(vars)
> data.new <- data.frame(lncpi, lnneer)
> VARselect(data.new, lag.max=10, type="const")
$selection
AIC(n) HQ(n) SC(n) FPE(n)
4 2 2 4
$criteria
1 2 3 4 5
AIC(n) -1.831254e+01 -1.844750e+01 -1.841657e+01 -1.851707e+01 -1.848107e+01
HQ(n) -1.825794e+01 -1.835651e+01 -1.828919e+01 -1.835329e+01 -1.828090e+01
SC(n) -1.817817e+01 -1.822355e+01 -1.810304e+01 -1.811396e+01 -1.798838e+01
FPE(n) 1.114228e-08 9.736147e-09 1.004341e-08 9.085414e-09 9.422136e-09
6 7 8 9 10
AIC(n) -1.842969e+01 -1.838650e+01 -1.838871e+01 -1.844871e+01 -1.839595e+01
HQ(n) -1.819311e+01 -1.811354e+01 -1.807935e+01 -1.810295e+01 -1.801380e+01
SC(n) -1.784741e+01 -1.771465e+01 -1.762727e+01 -1.759769e+01 -1.745536e+01
FPE(n) 9.924624e-09 1.037063e-08 1.035831e-08 9.767702e-09 1.031338e-08
根据结果,不同的信息准则有不同的滞后阶数,选择 2 阶或者 4 阶都是可以的,一般来说选择在相同条件下更加简洁的模型,因此选择 2 阶滞后。
VAR模型的拟合和预测
在确定好最优滞后阶数以后我们就可以拟合模型:
> var <- VAR(data.new, lag.max=2, ic="AIC")
> summary(var)
VAR Estimation Results:
=========================
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.1847
Call:
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 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.007369 on 130 degrees of freedom
Multiple 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 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.01285 on 130 degrees of freedom
Multiple 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 lnneer
lncpi 5.430e-05 7.761e-06
lnneer 7.761e-06 1.650e-04
Correlation matrix of residuals:
lncpi lnneer
lncpi 1.00000 0.08198
lnneer 0.08198 1.00000
根据模型拟合结果,所得的 VAR(2) 的模型方程(保留两位小数)为:

也可以用 coef 函数来查看模型估计的简明结果。
> coef(var)
$lncpi
Estimate Std. Error t value Pr(>|t|)
lncpi.l1 0.80917274 0.08578001 9.433115 2.063276e-16
lnneer.l1 -0.09248192 0.04705422 -1.965433 5.149594e-02
lncpi.l2 0.12906563 0.08483246 1.521418 1.305837e-01
lnneer.l2 0.08381139 0.04712463 1.778505 7.765840e-02
const 0.32677229 0.15855122 2.060989 4.129779e-02
$lnneer
Estimate Std. Error t value Pr(>|t|)
lncpi.l1 -0.1842958 0.14954402 -1.232385 2.200295e-01
lnneer.l1 1.3127548 0.08203166 16.003025 1.594736e-32
lncpi.l2 0.2885006 0.14789212 1.950750 5.323813e-02
lnneer.l2 -0.3199781 0.08215442 -3.894837 1.564088e-04
const -0.4472580 0.27640922 -1.618101 1.080647e-01
拟合出一致的 VAR(2) 方程(系数保留两位小数)。
> plot(var)
Hit <Return> to see next plot:
变量 lncpi 的时序图、残差图、ACF 图、PACF 图

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

在查看完拟合结果的图形之后,绘制拟合结果的脉冲响应图
> var.irf <- irf(var)
> plot(var.irf)
Hit <Return> to see next plot:


根据图形 lncpi 自身以及 lnneer 的波动对 lncpi 有正向的冲击。 lncpi 对自身的影响没有滞后期,并且自身波动的影响随着时期的增加会越来越小。lnneer 波动对 lncpi 的影响在第一期以前是逐渐减少的,但在第一期减少为 0 之后随着时间的增加影响越来越大。
lncpi 波动对 lnneer 有一个负向的冲击,并且随着时期数的增加负向的影响会越来越大。lnneer 波动对自身有一个正向的冲击,这个冲击从开始先增加,在第二期达到最大值以后又逐渐减少。
VAR(2)模型的预测
> var.predict <- predict(var, n.ahead=10, ci=0.95)
> var.predict
$lncpi
fcst lower upper CI
[1,] 4.624716 4.610273 4.639159 0.01444318
[2,] 4.623953 4.605348 4.642558 0.01860518
[3,] 4.623361 4.601393 4.645329 0.02196807
[4,] 4.622861 4.598303 4.647420 0.02455823
[5,] 4.622416 4.595772 4.649060 0.02664388
[6,] 4.622011 4.593662 4.650360 0.02834876
[7,] 4.621639 4.591878 4.651401 0.02976171
[8,] 4.621298 4.590354 4.652242 0.03094400
[9,] 4.620983 4.589043 4.652924 0.03194066
[10,] 4.620695 4.587909 4.653481 0.03278571
$lnneer
fcst lower upper CI
[1,] 4.814544 4.789365 4.839724 0.02517943
[2,] 4.816007 4.774543 4.857471 0.04146393
[3,] 4.816488 4.761790 4.871186 0.05469804
[4,] 4.816540 4.750928 4.882152 0.06561178
[5,] 4.816376 4.741493 4.891258 0.07488244
[6,] 4.816081 4.733118 4.899045 0.08296352
[7,] 4.815693 4.725534 4.905853 0.09015978
[8,] 4.815230 4.718554 4.911906 0.09667590
[9,] 4.814702 4.712050 4.917355 0.10265273
[10,] 4.814116 4.705927 4.922306 0.10818965
得到了两变量(lncpi、lnneer)的 VAR(2)模型的滞后 10 期预测结果。