一、自回归模型（AR模型）

1.1 概念

自回归模型（英语：Autoregressive model，简称AR模型），是统计上一种处理时间序列的方法，用同一变量例如X的之前各期，亦即至X_{t-1}来预测本期X_{t}的表现，并假设它们为线性关系。因为这是从回归分析中的线性回归发展而来，只是不用X预测Y，而是用X预测X（自己）；所以叫做自回归。自回归模型被广泛运用在经济学、资讯学、自然现象的预测上

1.2 定义

AR(p)模型的数学描述如下：

X_{t}=c+\sum_{i=1}^{p} \varphi_{i} X_{t-i}+\varepsilon_{t} \\

其中：c是常数项，\varphi_{i}是X_t与X_{t-i}的相关系数。\varepsilon_{t}是误差项，满足：

\begin{array}{l}\mathrm{E}\left(\varepsilon_{t}\right)=0 \rightarrow \text { 误差项的均值为0 } \\\mathrm{E}\left(\varepsilon_t^{2}\right)=\sigma^2\rightarrow \text { 误差项的方差为$\sigma^2$ } \\\mathrm{E}\left(\varepsilon_{t}\varepsilon_{t-k}\right)=0(k \ne 0)\rightarrow \text { 误差项不存在自相关 }\end{array} \\

用文字可以总结为：X的当期值等于一个或数个前期值的线性组合，加常数项，加随机误差。

二、向量自回归模型（VAR模型）

2.1 概念

向量自回归模型（英语：Vector Autoregression model，简称VAR模型）是一种常用的计量经济模型。它扩充了只能使用一个变量的自回归模型（简称：AR模型），使容纳大于1个变量，因此经常用在多变量时间序列模型的分析上。

2.2 定义

VAR(p)模型的数学描述如下：

y_{t}=c+A_{1} y_{t-1}+A_{2} y_{t-2}+\cdots+A_{p} y_{t-p}+e_{t} \\

其中：c是n \times 1常数向量，A_i(1\leq i \leq p)是n × n矩阵。 e_t是n \times 1误差向量，满足：

\begin{array}{l} \mathrm{E}\left(e_{t}\right)=0 \rightarrow \text { 误差向量的均值为0向量 } \\ \mathrm{E}\left(e_{t} e_{t}^{\prime}\right)=\Omega\rightarrow \text { 误差向量的协方差矩阵为$\Omega$（一个$n \times n$正定矩阵） } \\ \mathrm{E}\left(e_{t} e_{t-k}^{\prime}\right)=0(k \ne 0)\rightarrow \text { 误差向量不存在自相关 } \end{array} \\

2.3 VAR模型的步骤

1、选择合适的变量；

2、检验序列平稳性，看序列是否平稳，或者同阶单整；

3、多元混成检验；

4、模型定阶，根据AIC、BIC等准则，选择VAR模型的滞后阶数；

5、模型诊断；

6、协整检验，看变量间是否存在协整关系；

7、格兰杰因果检验；

8、脉冲响应，看变量对外界冲击的反馈；

9、方差分解；

10、结果预测。

三、VAR模型的R语言实现

下面以一个四元的时间序列数据为例，构建VAR模型并进行分析和预测。

3.1 变量选择

根据实际问题或理论，得到感兴趣的变量数据。

3.2 检验序列的平稳性

3.2.1 时序图

一元时间序列的时序图可以用最基本的plot函数，直接传入用ts函数构造的ts类型对象。

data.ts <- ts(X, start, end, frequency)
plot(data.ts, main, ylab, xlab)

若是多元时间序列，建议用xtsExtra包的plot.xts函数。xts 是一种时间序列数据类型，既可以保存等间隔时间序列数据，也可以保存不等间隔的时间序列数据。

首先用as.xts函数将数据转为xts类型。as.xts传入两个参数，第一个参数是不包含时间的数据，第二个参数传入时间数据。

接着用plot.xts函数画出时间序列图。这里介绍下例子中出现的参数，首先传入一个xts类型对象，col参数控制序列线条的颜色，auto.legend用于控制是否自动加图例，main是图像的标题，legend.loc用于控制图例的位置，major.ticks用于控制主要刻度线，auto.grid是控制画网格线的。

data.xts <- xtsExtra::as.xts(data[,-1], as.Date(data[,1]))
plot(data.xts, col =c(1,2,3,4),auto.legend = T, main="时间序列图",
legend.loc = "topright",major.ticks="years",auto.grid=T)

从时间序列图大概可以看出，B序列和C序列是非平稳序列，考虑差分。

3.2.2 单位根检验

用tseries包中的adf.test函数进行更精确的平稳性检验，备择假设为序列平稳。检验结果如图。

tseries::adf.test(data.xts[,1])

因为我们只关心p值，因此可以通过调用检验的p.value属性，并加多个序列的p值放到一个data frame中，便于阅读查看。

library(tseries)
p1 <- round(adf.test(data.xts[,1])$p.value,3)
p2 <- round(adf.test(data.xts[,2])$p.value,3)
p3 <- round(adf.test(data.xts[,3])$p.value,3)
p4 <- round(adf.test(data.xts[,4])$p.value,3)
p_df <- data.frame(matrix(c(p1,p2,p3,p4),nrow=1))
colnames(p_df2) <- c("A","B","C","D")
rownames(p_df2) <- c("单位根检验的p值")
p_df

检验结果和之前时序图得出的结论一致，B序列和C序列不平稳。

3.3 多元混成检验

将一元的 Ljung-Box 白噪声检验推广到了多元的情形即为多元混成检验。对一个多元时间序列，检验原假设：

H_{0}: \rho_{1}=\cdots=\rho_{m}=0 \\

备择假设是不全为零矩阵。这可以检验多元时间序列 r_t为宽白噪声的零假设，即 r_t为弱平稳列且无序列自相关，可以有同步的分量间相关。

用 MTS::mq() 计算多元混成检验，返回不同滞后阶数的统计量计算值及对应的p值，同时返回图像。直接看图像，不超过蓝色虚线即为拒绝原假设。因此，序列不是白噪声序列，有分析价值。

mq(data.xts, lag=12)

3.4 模型定阶

用coredata函数提取出xts类型的核心数据，也就是除了时间数据的数据后，用MTS包的VARorder函数即可返回根据AIC、BIC及HQ信息准则选择的模型阶数。

例子中根据AIC准则选择5阶，根据BIC准则选择1阶，根据HQ准则选择1阶。综合考虑，同时依据尽量低阶的原则，选择建立1阶模型。

Z <- xtsExtra::coredata(data.xts)
m <- MTS::VARorder(Z)

接下来用vars包中的VAR函数可以进行VAR模型的估计。MTS包中也有用来构建VAR模型的VAR函数，但该函数的返回结果是矩阵，不如vars::VAR函数直观。因此，这里选择后者。

var1 <- vars::VAR(Z,p=1)
var1$varresult$'A'

3.5 模型诊断

Eviews中一般看VAR模型根是否在单位圆内。若VAR模型的根都在单位圆内，才可继续后续的分析。

在R中暂时没有找到这样的操作语句，只发现通过残差累积和检验判断的方法。

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

diagnostic_t <- stability(var1, type = c("OLS-CUSUM"), h = 0.15, dynamic = FALSE, rescale = T)
par(mfrow=c(2,2))
plot(diagnostic_t$stability$"A",main='OLS-base CUSUM test：A')
plot(diagnostic_t$stability$"D",main='OLS-base CUSUM test：D')
plot(diagnostic_t$stability$"diff_B",main='OLS-base CUSUM test：diff_B')
plot(diagnostic_t$stability$"diff_C",main='OLS-base CUSUM test：diff_C')

残差检验和检验显示，模型结果稳定。

3.6 协整检验

本例中，变量并未存在同阶单整的关系，因此不考虑协整检验。

3.7 格兰杰因果检验

格兰杰因果检验的原假设是关心的自变量不是因变量的格兰杰原因，因此我们希望拒绝原假设。需要注意格兰杰因果检验只是统计上的因果关系，经济意义上的因果关系还需进行进一步分析。

vars::causality(var1,cause=c('diff_B', 'diff_C', 'D'))

得到的p值很小，拒绝原假设，格兰杰因果性成立。

3.8 脉冲响应

脉冲响应分析，就是某一内生变量对于残差冲击的反应。具体而言，他描述的是在随机误差项上施加一个标准差大小的某一内生变量冲击后，对感兴趣的内生变量的当期值和未来值所产生的影响。

name <- names(train_set)

par(mfrow=c(2,2))
for(i in 1:1)
{
p=vars::VAR(Z, p=1)
svec.irf <- irf(p,response=name[i], n.ahead = 150, boot = TRUE)
##利用循环语句，得出每一个脉冲响应的图，整齐地排列到画布上
for(j in 1:4)
{
p=as.vector(svec.irf$irf[[j]])
q=as.vector(svec.irf$Upper[[j]])
k=as.vector(svec.irf$Lower[[j]])
low=min(k)-0.1
high=max(q)+0.1
plot(p,type='l',main=paste(paste(svec.irf$response,'from'),svec.irf$impulse[j]),
ylim=c(low,high))
lines(q,type='l',lty=2,col='red')
lines(k,type='l',lty=2,col='red')
abline(h=0,col='red')
}
}

从上图可以得知： ① 正向的D冲击，使得A在第25天左右达到正的最大值，随后持续下降到第120期时趋于0； ② 正向的diff_B冲击几乎不会引起A变动； ③ 正向的diff_C冲击效果，使得A在第25天左右达到负的最大值，随后持续下上升到第100期时趋于0。

3.9 方差分解

VAR模型的应用，还可以采用方差分解方法研究模型的动态特征。

方差分解是进一步评价各内生变量对预测方差的贡献度。方差分解是分析预测残差的标准差由不同新息的冲击影响的比例，亦即对应内生变量对标准差的贡献比例。

fevd1<-fevd(var1, n.ahead = 100)
df <- data.frame(fevd1$'A')
df['lag'] <- seq(1,100)

y <- df[,1]
x1 <- df[,2]
x2 <- df[,3]
x3 <- df[,4]
lag <- df[,5]

par(mfrow=c(2,2))
plot(lag,y,type='l',col=1,ylab='预测方差贡献率',main="A")
plot(lag,x1,type='l',col=1,,ylab='预测方差贡献率',main='D')
plot(lag,x2,type='l',col=1,ylab='预测方差贡献率',main='diff_B')
plot(lag,x3,type='l',col=1,ylab='预测方差贡献率',main='diff_C')

短期内，A的变动由自身主导，但随着时间推移，其对自身变动的贡献率逐渐降低，D对A变动的贡献率逐渐增加。到了第70期左右时，A对自身变动的贡献率稳定在33%左右，D对A变动的贡献率稳定在67%左右。而diff_B和diff_C对A变动的影响很小，可忽略不计。

3.10 结果预测

若想通过构建的VAR模型进行预测，使用vars::predict方法即可。

若想评价模型预测效果，我们就需要分隔数据集。这里以4：1的比例将原数据集分割为训练集和测试集。

n <- dim(data.xts)[1]
#xts类型用first、last方法取子集
train_set <- first(data.xts,round(0.8*n))
test_set <- last(data.xts,n-round(0.8*n))
h <- dim(test_set)[1]

之后以训练集训练得到模型var1，并以该模型为基础进行预测。

var1.predict <- predict(var1,n.ahead=h,ci=0.95)
pred <- var1.predict$fcst$'A'[,1]
lower <- var1.predict$fcst$'A'[,2]
higher <- var1.predict$fcst$'A'[,3]
plot(1:length(pred),pred,type='l',ylim=c(36,110))
lines(1:length(pred),lower,col=2,type='l',lty=2)
lines(1:length(pred),higher,col=2,type='l',lty=2)
lines(1:length(pred),coredata(test_set)[,1],col=4,type='l')

从结果可以看到，VAR模型的预测效果与一元的AR模型类似，较平滑，且后期趋于一定值。

四、参考资料

1、wikipedia：https://www.wikipedia.org/

2、CSDN博客：https://blog.csdn.net/Imliao/article/details/80352158?utm_medium=distribute.pc_relevant_t0.none-task-blog-BlogCommendFromMachineLearnPai2-1.compare&depth_1-utm_source=distribute.pc_relevant_t0.none-task-blog-BlogCommendFromMachineLearnPai2-1.compare

3、金融时间序列分析讲义（李东风）：http://www.math.pku.edu.cn/teachers/lidf/course/fts/ftsnotes/html/_ftsnotes/index.html