北京大学金融时间序列分析讲义第24章: 向量自回归模型

24.1 VAR(1)模型

多个资产收益率的联合模型中最常用的是向量自回归
(Vector Autoregression, VAR)模型。
k元时间序列rt服从一个一阶向量自回归模型
即VAR(1),若有

rt=ϕ0+Φrt1+at(24.1)

其中ϕ0k维常数向量,
Φk阶常数方阵,
{at}是序列不相关的弱平稳列,
Eat=0,
Var(at)=Σ>0(对称正定矩阵)。

文献中经常假定{at}服从多元正态分布
Nk(0,Σ)
则这时at是独立同分布随机向量序列。

at=(a1t,,akt)T
ϕ0=(ϕ10,,ϕk0)T,
Φ=(ϕij)k×k

24.2 模型结构和格兰杰因果性

考虑k=2的情形。
模型变成

{r1t=ϕ10+ϕ11r1,t1+ϕ12r2,t1+a1tr2t=ϕ20+ϕ21r1,t1+ϕ22r2,t1+a2t(24.2)

如果ϕ12=ϕ21=0,
r1tr2t两个序列是分离的:

r1t=r2t=ϕ10+ϕ11r1,t1+a1tϕ20+ϕ22r2,t1+a2t

这时两个序列分别服从单独的一元AR(1)模型;
如果a1ta2t不相关,
{r1t}序列与{r2s}序列不相关。
称这样的分离的序列为非耦合的。

如果ϕ120, ϕ210,
这时两个序列之间有相互反馈关系。

如果ϕ12=0, ϕ210,
模型变成

{r1t=ϕ10+ϕ11r1,t1+a1tr2t=ϕ20+ϕ21r1,t1+ϕ22r2,t1+a2t(24.3)

这时r1t不受到r2t的过去值的影响,
r2t受到r1,t1的过去值的影响。
可以将r1t作为输入变量,
r2t作为输出变量。

在统计学文献中,
如果a1ta2t也不相关,
称式(24.3)r1tr2t
之间有传递函数(transfer function)关系。
传递函数模型是VARMA的一种特殊形式,
在控制工程中非常有用,
可以通过调整r1t的值来影响r2t
在计量经济学文献中,
此模型意味着两个序列之间存在格兰杰因果关系,
r1t的过去值影响r2t
r2t的过去值不影响r1t

(Granger 1969)提出了因果关系的概念,
适用于解释VAR模型的结果。
考虑一个二元序列rt=(r1t,r2t)T
的超前l步预测问题。
可以分别用VAR模型和每个分量的一元模型来进行预测。
如果r2t的二元预测比它的一元预测更准确,
则称r1tr2t的格兰格原因,
这是因为r1t的信息提高了r2t的预测精度。
这里预测精度用预测的均方误差来度量。
当然,r2t也可以同时是r1t的格兰格原因。

Ft表示截止到时刻t的可用信息,
Ft包含{rt,rt1,rt2,}
Fi,t表示Ft中去掉第i个分量的信息的集合。
考虑(24.2)的二元VAR(1)模型。
Ft包含{r1t,r2t,r1,t1,r2,t1,},
F1,t包含{r2t,r2,t1,}
考虑基于Ft的超前l步预测,
其预测误差为

et(l)=(e1t(l),e2t(l))T=rt+lE(rt+l|Ft)

其中对r2,t+l的预测误差为

e2t(l)=r2,t+lE(r2,t+l|Ft)

考虑基于F1,tr2,t+l的预测,
其预测误差为

ẽ 2t(l)=r2,t+lE(r2,t+l|F1,t)=r2,t+lE(r2,t+l|r2,t,r2,t1,)

如果

Ee22t(l)=Var(r2,t+l|Ft)<Eẽ 22t(l)=Var(r2,t+l|F1,t)

即使用截止到t时刻为止的全部信息预测r2,t+l
均方误差比不利用第一个分量序列的均方误差小,
则称r1tr2t的格兰格原因。

回到二元VAR(1)模型(24.2)
如果ϕ12=0ϕ210
r2,t+1依赖于r1,t
预测r2,t+1时需要利用r1t的信息,
序列r1t是序列r2t的格兰格原因;
r1,t+1不依赖于r2,t
预测r1,t+1就不需要用到r2t
所以序列r2t不是序列r1t的格兰格原因。

事实上,对(24.3)
不妨设{at}是独立的多元弱平稳列,
这时
E(at+1|Ft)=0,
E(a2,t+1|F1,t)=0,
Var(a2,t+1|Ft)=Var(a2,t+1|F1,t)=σ22

E(r2,t+1|Ft)=Var(r2,t+1|Ft)==ϕ20+ϕ21r1,t+ϕ22r2,tE[(r2,t+1E(r2,t+1|Ft))2|Ft]E(a22,t+1|Ft)=σ22


E(r2,t+1|F1,t)=Var(r2,t+1|F1,t)===ϕ20+ϕ21E(r1,t|F1,t)+ϕ22r2,tE[(r2,t+1E(r2,t+1|F1,t))2|F1,t]E[{ϕ21(r1,tE(r1,t|F1,t))+a2,t+1}2|F1,t]ϕ221Var(r1,t|F1,t)+σ22>σ22

r1tr2t的格兰格原因。

另一方面,

E(r1,t+1|Ft)=Var(r1,t+1|Ft)==E(r1,t+1|F2,t)=Var(r1,t+1|F2,t)==ϕ10+ϕ11r1,tE[(r1,t+1E(r1,t+1|Ft))2|Ft]E(a21,t+1|Ft)=Var(a1,t+1)=σ11ϕ10+ϕ11r1,tE[(r1,t+1E(r1,t+1|F2,t))2|F2,t]E(a21,t+1|F2,t)=Var(a1,t+1)=σ11

所以r2t不是r1t的格兰格原因。

如果ϕ21=0ϕ120
则序列r2t是序列r1t的格兰格原因,
而序列r1t不是序列r2t的格兰格原因。

如果新息at=(a1t,a2t)T的协方差阵Σ不是对角阵,
r1tr2t之间存在同步相关性,
这时两个序列存在即期(同期,瞬时)格兰格因果关系,
这种即期因果关系是双向的。

注意,
这样定义的格兰格因果性依赖于VAR(1)模型(24.1)这样的模型形式(称为简化形式),
如果改变模型的形式(比如改为下面的结构形式),
可能就不会在ϕ12ϕ21这样的简单数值上反映出格兰格因果性。

更多分量的VAR(1)可以有更灵活的格兰格因果关系。
比如三元的模型,
如果Φ是下三角阵:

ϕ11ϕ21ϕ31ϕ22ϕ32ϕ33

则:

  • r2tr3t都不是r1t的格兰格原因;
  • r1tr2t的格兰格原因,
    r3t不是r2t的格兰格原因;
  • r1tr2t都是r3t的格兰格原因。

如果Φ有其他稀疏形式,
关系可能多种多样。

24.3 VAR的简化形式和结构形式

(24.1)的系数矩阵Φ
度量了rt的动态相依性,
r1tr2t之间的同步相依性可以通过
at的协方差阵Σ的非对角元素
σ12来反映。
如果σ12=0
则两个分量r1tr2t之间没有同步线性关系。
在计量经济文献中式(24.1)中的VAR(1)
模型称为简化形式(reduced form)的模型,
因为该模型没有清楚地表现出分量序列之间的同步线性相依性。

可以利用矩阵变换对(24.1)进行变换,
得到显式地表现同步关系的模型。
Σ>0(正定),
存在Cholesky分解
Σ=LGLT,
其中G是对角线元素都为正数的k阶对角方阵,
L为对角线元素都等于1的下三角阵,
定义

bt=L1at=(b1t,,bkt)T


Ebt=Var(bt)=0L1Var(at)LT=G

其中LT=(LT)1=(L1)T
上式表明bt的各分量不相关。
将式(24.1)两边同时左乘L1

L1rt==L1ϕ0+L1Φrt1+L1atϕ0+Φrt1+bt(24.4)

ϕ0=(ϕ10,,ϕk0)T,
Φ=(ϕij)k×k
L1的最后一行为
(wk1,,wk,k1,1)
则模型(24.4)的最后一个方程(第k个方程)为

rkt+i=1k1wkirit=ϕk0+i=1kϕkiri,t1+bkt(24.5)

i<k
因为bktbki不相关,
所以式(24.5)明确表示了
rktrit的同步依赖性。
在计量经济文献中式(24.5)称为rkt
的一个结构方程(structural form)。

rt的任何其它分量rjt(j<k),
可以对VAR模型各个方程的次序进行重排,
比如,
将第j个方程与第k个方程对换位置,
再用前面的方法得到关于rjt的结构方程。
因此,
(24.1)的简化形式等价于式(24.5)的结构形式。

在时间序列分析中通常使用简化形式,因为

  • 简化形式更容易估计;
  • 在预测时,同步形式无法使用。

例2.1

为了说明从简化形式到结构方程的变换,
考虑如下的二元VAR(1)模型:

(r1tr2t)=Var(at)=(0.20.4)+(0.20.60.31.1)(r1,t1r2,t1)+(a1ta2t)Σ=(2111)

Σ作Cholesky分解得
Σ=LGLT,
其中

G=(2000.5)L1=(10.501)

在简化形式的VAR(1)模型两边同时左乘L1

(10.501)(r1tr2t)=Var(bt)=(0.20.3)+(0.20.70.30.95)(r1,t1r2,t1)+(b1tb2t)G=(2000.5)

其中第二个方程为

r2t=0.3+0.5r1t0.7r1,t1+0.95r2,t1+b2t

此方程明确表现了r2tr1t的同步依赖(系数+0.5),
方程中也有r2tr1,t1的交叉依赖,
r2,t1的自身依赖。

为了给出r1t的结构方程,
将简化形式的两个方程对换位置,
得(注意Φ需要行互换、列互换)

(r2tr1t)=Σ=(0.40.2)+(1.10.30.60.2)(r2,t1r1,t1)+(a2ta1t)(1112)

Σ的Cholesky分解有

G=(1001)L1=(1101)

在调换次序的VAR(1)模型的两边同时左乘L1,得

(1101)(r2tr1t)=Var(ct)=(0.40.2)+(1.10.80.60.8)(r2,t1r1,t1)+(c1tc2t)G=(1001)

其中第二个方程为

r1t=0.2+1.0r2t0.8r2,t1+0.8r1,t1+c2t

这个方程明确地表现了r1tr2t的同步线性依赖。

○○○○○

利用R计算LGLT分解的L1G的程序:

ldl.decomp <- function(A){
  R <- chol(A)
  D <- diag(diag(R))
  Linv <- solve(t(solve(D, R)))
  G <- D^2
  
  list(Linv=Linv, G=G)
}

例如

Sigma <- rbind(c(2,1), c(1,1)); Sigma
##      [,1] [,2]
## [1,]    2    1
## [2,]    1    1
## $Linv
##      [,1] [,2]
## [1,]  1.0    0
## [2,] -0.5    1
## 
## $G
##      [,1] [,2]
## [1,]    2  0.0
## [2,]    0  0.5
Sigma <- rbind(c(1,1), c(1,2)); Sigma
##      [,1] [,2]
## [1,]    1    1
## [2,]    1    2
## $Linv
##      [,1] [,2]
## [1,]    1    0
## [2,]   -1    1
## 
## $G
##      [,1] [,2]
## [1,]    1    0
## [2,]    0    1

24.4 VAR(1)模型的平稳性条件和矩

(24.1)式中的VAR(1)模型的rt是弱平稳的,
在方程两边取期望得

Ert=ϕ0+ΦErt1

Ert=Ert1
IΦ满秩(也叫做非奇异)时

μ=Ert=(IΦ)1ϕ0

在式(24.1)中代入
ϕ0=(IΦ)μ,
则模型(24.1)可以写成

rtμ=Φ(rt1μ)+at

记中心化的rtr̃ t=rtμ
则VAR(1)模型(24.1)变成

r̃ t=Φr̃ t1+at(24.6)

(24.6)递归得

r̃ t=at+Φat1+Φ2at2+=j=0Φjatj(24.7)

其中Φ0=Ik
这相当于一元时间序列的Wold表示。

(24.7)可以推出:

  • 因为{at}序列不相关,
    所以Cov(at,rtl)=0(l>0时)。
    因此可以将at称为序列{rt}t
    时刻的一个扰动或者新息(innovation)。
    一般定义新息为at=rtE(rt|Ft1)
    at=rtL(rt|Ht1)
    L(|)表示均方误差最小线性预测,
    Ht1表示由{rt1,rt2,}张成的Hilbert空间。
    平稳的VAR(1)中ait是用
    rt1,rt2,rit作最优线性预测的预测误差。
  • (24.7)两边同时右乘aTt后取期望,
    Cov(rt,at)=Σ
  • 对VAR(1)模型(24.1)导出的Wold表示(24.7)
    rt按系数矩阵Φj依赖于过去的新息atj
    为了使得(24.7)收敛,
    需要jΦj0
    这当且仅当Φ的所有特征值的模都小于1。
    如果存在模大于或等于1的特征值,
    jΦj发散或者不收敛到0矩阵。
    仿照一元情形的特征多项式的讨论,
    因为特征值的模小于1当且仅当
    |λIΦ|的根的模都小于1,
    |λIΦ|=λk|IΦλ1|
    z=λ1
    称关于z的多项式|P(z)|=|IΦz|为模型的(逆序)特征多项式
    则式(24.7)收敛的充分必要条件是特征多项式|P(z)|的根都在单位圆外。
    这个条件称为VAR(1)模型的平稳性条件
    当且仅当平稳性条件成立时,
    模型(24.1)存在满足
    atrt1,rt2,不相关条件的平稳解。
  • 利用推移算子和特征多项式可以将模型(24.1)写成
    P(B)rt=ϕ0+at
  • 利用Wold表示(24.7)式还可得

Var(rt)=Σ+ΦΣΦT+Φ2Σ(Φ2)T+=j=0ΦjΣ(Φj)T

考虑VAR(1)模型的互协方差阵和互相关阵。
将式(24.6)两边乘以r̃ Ttl后取期望,
得互协方差阵

Γl===E(r̃ tr̃ Ttl)=Cov(rt,rtl)ΦE(r̃ t1r̃ Ttl)+E(atr̃ tl)ΦΓl1,(l>0)(24.8)

这个结果是一元时间序列AR(1)模型性质的推广。

通过递推可得

Γl=ΦlΓ0,l>0

其中Γ0=Var(rt)

DΓ0的对角阵,
则VAR(1)模型的互相关阵为

ρl===D12ΓlD12=D12ΦΓl1D12D12ΦD12D12Γl1D12D12ΦD12ρl1=Υρl1

其中Υ=D12ΦD12
(Υ读作Upsilon)。
递推得

ρl=Υlρ0,l=1,2,

例2.2

考虑二元VAR(1)模型

(r1tr2t)=Σ=(53)+(0.20.60.31.1)(r1,t1r2,t1)+(a1ta2t)(10.80.82)

计算Φ的特征值的绝对值:

abs(eigen(rbind(c(0.2, 0.3), c(-0.6, 1.1)))$value)
## [1] 0.8 0.5

特征值的绝对值都小于1。
模型是平稳的。

注意ϕ22=1.1>1
但不影响平稳性。
多元弱平稳一定也是一元弱平稳的。

24.5 VAR(1)模型的边缘模型

对VAR(1)模型(24.1)的如上的Wold表示,
同时也将rt的每个分量,
表示成了一元的无穷阶MA,但其中的新息是多元的。
利用这种思想,将分量的一元模型可以写成一元ARMA模型。

设矩阵A是可逆方阵,
C=(cij)
其中cijA删除第i行和第j列后的子矩阵的行列式乘以(1)i+j
cij称为A(i,j)元素的代数余子式,
CT称为A的伴随矩阵,
记为adj(A)

A1=1|A|adj(A),Aadj(A)=adj(A)A=|A|I

VAR(1)的Wold表示相当于

rt=(IΦB)1(ϕ0+at)

利用伴随矩阵可得

|IΦB|rt=adj(IΦB)(ϕ0+at)

其中|IΦz|zk阶多项式;
adj(IΦz)的每个元素是zk1阶多项式。
这样,分量rit服从一个ARMA(k,q)模型,
模型中的MA部分涉及到多元的新息,
但是可以设法转换成一元的新息。

例2.3

考虑例2.2中的模型。

(r1tr2t)=Σ=(53)+(0.20.60.31.1)(r1,t1r2,t1)+(a1ta2t)(10.80.82)

这时

Φ=IΦz=adj(IΦz)=|IΦz|=(0.20.60.31.1)(10.2z0.6z0.3z11.1z)(11.1z0.6z0.3z10.2z)11.3z+0.4z2

于是

=(11.3B+0.4B2)(r1tr2t)(11.10.60.310.2)(53)+(11.1B0.6B0.3B10.2B)(a1ta2t)

所以分量r1t的一元模型为

r1t1.3r1,t1+0.4r1,t2=5.4+a1t+0.3a2,t11.1a1,t1

ηt=a1t+0.3a2,t11.1a1,t1
易见{ηt}是弱平稳列,
且其ACF在q=1之后截尾(即ρη(j)=0,j>1)。
所以{ηt}是MA(1)序列,
从而分量r1t服从一个一元的ARMA(2,1)模型。

一般地,
k元的VAR(1)模型的分量服从一元的ARMA(k,k1)模型。
由于方程左右可能有公共根,
所以模型阶数可能会降低。

24.6 VAR(p)模型

k元时间序列rt服从一个VAR(p)模型,
如果

rt=ϕ0+Φ1rt1++Φprtp+at(24.9)

其中ϕ0{at}同VAR(1)的规定,
Φk阶方阵(=1,2,,p),
(i,j)元素记为ϕij()
利用向后推移算子(滞后算子)B可以将模型写成

(IΦ1BΦpBp)rt=ϕ0+at


P(z)=IΦ1zΦpzp(24.10)

这是一个从复数zk阶方阵P(z)的变换,
P(z)的每个元素是关于z的阶数不超过p的多项式。
称一元多项式函数det(P(z))(或记为|P(z)|)为模型的(逆序)特征多项式。

利用矩阵推移函数P(B)可以将模型简记为

P(B)rt=ϕ0+at

如果rt弱平稳,
则在下式的逆矩阵存在的条件下

μ=Ert=(IΦ1Φp)1ϕ0=(P(1))1ϕ0

中心化rtr̃ t=rtμ
这时模型(24.9)变成

r̃ t=Φ1r̃ t1++Φpr̃ tp+at(24.11)

称VAR(p)满足平稳性条件,
如果特征多项式det(P(z))的根都在单位圆外。
当VAR(p)满足平稳性条件时,
有如下的Wold表示的平稳解:

rt=μ+j=0Ψjatj,

其中jΨj的元素以负指数速度趋于0,
Ψj有如下的递推公式:

Ψ0=Ψj=I,s=0pΦsΨjs,

s<0Ψs=0

由上述Wold表示可知,
对平稳解rt

  • Cov(rt,at)=Σ=Var(at);
  • Cov(rtl,at)=0(l>0)。

互协方差阵满足如下递推关系

Γl=Φ1Γl1++ΦpΓlp, l>0

这个性质称为VAR(p)的矩方程
是一元AR(p)模型的Yule-Walker方程的多元形式。
用CCM表示为

ρl=Υ1ρl1++Υpρlp, l>0

其中Υj=D12ΦjD12
DΣ的对角阵。

事实上,VAR(p)模型可以改写成一个VAR(1)模型,
然后可以用VAR(1)模型的结果去研究VAR(p)模型的性质。

xt=Var(bt)=Φ=r̃ tr̃ t1r̃ tp+1bt=at0k×10k×1(Σ0k(p1)×k0k×k(p1)0k(p1)×k(p1))Φ1I00Φ20I0Φp100IΦp000


xt=Φxt1+bt(24.12)

是关于kp元时间序列xt的VAR(1)模型。

rt弱平稳当且仅当xt弱平稳,
其充分必要条件是Φ的所有特征值模小于1,
计算可知这相当于要求用行列式表示的z的多项式|P(z)|的根都在单位圆外,
P(z)定义见(24.10)
但是,
因为计算det(IΦ1zΦpzp)这种用行列式表示的多项式的根还需要求多项式的各个系数,
而求Φ的所有特征值或模最大的特征值则有比较成熟的数值方法,
所以为了判断一个VAR(p)模型是否满足平稳性条件,
可以生成Φ并判断其模最大的特征值是否模小于1。

如果det(P(1))=0
Φ有等于1的特征值,
即有单位根,
rt的各个分量中至少有一个是单位根过程,
各个分量之间还有可能存在协整关系,
可以用VECM(向量误差修正模型)建模,
参见25

对VAR(p)模型,
系数矩阵Φl=(ϕij(l))k×k的值也表现了
rt各个分量之间的先导、滞后关系。
比如k=2时,
如果所有ϕ12(l)=0
但存在l>0使得ϕ21(l)0
r1t不依赖于r2t的过去值,
r2t依赖于r1t的过去值。
r1tr2t的格兰杰原因,
r2t不是r1t的格兰杰原因。

24.7 估计和定阶

VAR模型建模也基本遵循定阶、模型估计和模型检验这样的反复尝试过程。
一元的PACF可以推广到多元情形用以辅助定阶。

考虑如下的阶数递进的VAR模型:

rt=rt=rt=ϕ0+Φ1rt1+atϕ0+Φ1rt1+Φ2rt2+atϕ0+Φ1rt1++Φprtp+at(24.13)

模型参数可以对每个方程分别用OLS(最小二乘)方法估计,
在多元统计分析文献中称为多元线性估计,
参见(Johnson and Wichern 1998)

对于(24.13)的第i个方程,
Φ̂ (i)jΦj的最小二乘估计,
ϕ̂ (i)0ϕ0的的最小二乘估计。
残差为

â (i)t=rtΦ̂ (i)1rt1Φ̂ (i)irti

i=0时定义â (0)t=rtr¯
其中r¯rt的样本均值。

残差的方差阵估计为

Σ̂ i=1T(k+1)i1t=i+1Tâ (i)t[â (i)t]T

(参见(Tsay 2014) 节2.5.1 P.34。)

为了定阶,
可以对l=1,2,逐一地检验
H0:Φl=0Ha:Φl0

l=1,要检验
H0:Φ1=0Ha:Φ10,
使用检验统计量

M(1)=(Tk52)ln|Σ̂ 1||Σ̂ 0|

在一些正则性条件下当H0成立时M(1)渐近服从χ2(k2)分布,参见(Tiao and Box 1981)
计算右侧p值。

一般地,
可以通过(24.13)的第l1个方程和第l个方程检验
H0:Φl=0Ha:Φl0
取检验统计量

M(l)=(Tlkl32)ln|Σ̂ l||Σ̂ l1|

H0M(l)渐近服从χ2(k2)分布,
计算右侧p值。

另一种定阶方法是利用AIC或其它类似的信息准则。
{at}是独立同多元正态分布的多元白噪声,
可以用最大似然(ML)方法估计模型,
对于VAR模型,
最小二乘估计ϕ̂ 0Φ̂ j
与条件最大似然估计等价,
但误差项方差阵Σ的ML估计为

Σ̃ i=1Tt=i+1Tâ (i)t[â (i)t]T

VAR(i)模型在正态假定下的AIC值定义为

AIC(i)=ln|Σ̃ i|+2k2iT

取最终的阶p使得
AIC(p)=min0ip0AIC(i)
其中p0是预先规定的可取的最大阶数。

类似的信息准则还有

BIC(i)=HQ(i)=ln|Σ̃ i|+k2ilnTTln|Σ̃ i|+k2ilnlnTT

HQ准则是Hannan和Quinn(1979)提出的。

这些准则的选择结果不受量纲的影响,
对数据乘以常数λ
结果会使得ln|Σ̃ i|固定地增加ln(λ2k)

24.7.1 英、加、美GDP的VAR建模

考虑英国、加拿大、美国GDP的季度增长率,
时间从1980年第二季度到2011年第二季度。
数据经过季节调整,
来自圣路易斯联邦储备银行的数据库。
GDP是以本国货币百亿为单位,
增长率表示为GDP的对数的差分。

读入数据,计算对数增长率:

da <- read_table("q-gdp-ukcaus.txt")
## 
## -- Column specification --------------------------------------------------------
## cols(
##   year = col_double(),
##   mon = col_double(),
##   uk = col_double(),
##   ca = col_double(),
##   us = col_double()
## )
ts.gdp3 <- ts(as.matrix(da[,c("uk", "ca", "us")]), 
              start=c(1980,1), frequency=4)
ts.gdp3r <- diff(log(ts.gdp3))*100

三个国家的季度GDP对数增长率的时间序列图:

plot(as.xts(ts.gdp3r), type="l", 
     multi.panel=TRUE, theme="white",
     main="英国、加拿大、美国GDP的季度增长率(%)",
     major.ticks="years",
     grid.ticks.on = "years")

英国、加拿大、美国GDP的季度增长率

图24.1: 英国、加拿大、美国GDP的季度增长率

利用MTS包的VAR()函数估计VAR(1)模型:

library(MTS, quietly = TRUE)
Z <- coredata(as.xts(ts.gdp3r))
m1.gdp3r <- MTS::VAR(Z, 1)
## Constant term: 
## Estimates:  0.1713324 0.1182869 0.2785892 
## Std.Error:  0.06790162 0.07193106 0.07877173 
## AR coefficient matrix 
## AR( 1 )-matrix 
##       [,1]  [,2]   [,3]
## [1,] 0.434 0.189 0.0373
## [2,] 0.185 0.245 0.3917
## [3,] 0.322 0.182 0.1674
## standard error 
##        [,1]   [,2]   [,3]
## [1,] 0.0811 0.0827 0.0872
## [2,] 0.0859 0.0877 0.0923
## [3,] 0.0940 0.0960 0.1011
##   
## Residuals cov-mtx: 
##            [,1]       [,2]       [,3]
## [1,] 0.28933472 0.01965508 0.06619853
## [2,] 0.01965508 0.32469319 0.16862723
## [3,] 0.06619853 0.16862723 0.38938665
##   
## det(SSE) =  0.02721916 
## AIC =  -3.459834 
## BIC =  -3.256196 
## HQ  =  -3.377107

估计VAR(2)模型:

Z <- coredata(as.xts(ts.gdp3r))
m2.gdp3r <- MTS::VAR(Z, 2)
## Constant term: 
## Estimates:  0.1258163 0.1231581 0.2895581 
## Std.Error:  0.07266338 0.07382941 0.0816888 
## AR coefficient matrix 
## AR( 1 )-matrix 
##       [,1]  [,2]   [,3]
## [1,] 0.393 0.103 0.0521
## [2,] 0.351 0.338 0.4691
## [3,] 0.491 0.240 0.2356
## standard error 
##        [,1]   [,2]   [,3]
## [1,] 0.0934 0.0984 0.0911
## [2,] 0.0949 0.1000 0.0926
## [3,] 0.1050 0.1106 0.1024
## AR( 2 )-matrix 
##         [,1]   [,2]     [,3]
## [1,]  0.0566  0.106  0.01889
## [2,] -0.1914 -0.175 -0.00868
## [3,] -0.3120 -0.131  0.08531
## standard error 
##        [,1]   [,2]   [,3]
## [1,] 0.0924 0.0876 0.0938
## [2,] 0.0939 0.0890 0.0953
## [3,] 0.1038 0.0984 0.1055
##   
## Residuals cov-mtx: 
##            [,1]       [,2]       [,3]
## [1,] 0.28244420 0.02654091 0.07435286
## [2,] 0.02654091 0.29158166 0.13948786
## [3,] 0.07435286 0.13948786 0.35696571
##   
## det(SSE) =  0.02258974 
## AIC =  -3.502259 
## BIC =  -3.094982 
## HQ  =  -3.336804

VAR(1)的AIC为3.46,VAR(2)的AIC为3.50,
VAR(2)占优。
rt的三个分量分别是英国、加拿大、美国的GDP对数增长率(单位为百分之一),
得到的模型VAR(2)模型可以写成

rt=Var(at)=0.130.120.29+0.390.350.490.100.340.240.050.470.24rt1+0.060.190.310.110.180.130.020.010.09rt2+at0.280.030.070.030.290.140.070.140.36

结果中也包括系数的标准误差,可以看出有些系数是不显著的。
比如,三个序列的次序为英国、加拿大、美国,
ϕ12(l), ϕ13(l)(l=1,2)都不显著,
说明英国的GDP季度增长率不受加拿大和美国的过去值的影响,
或者粗略地说,
加拿大和美国的GDP季度增长率不是英国的GDP季度增长率的格兰格原因。
更严格的分析应该进行假设检验。

利用MTS包的VARorder函数可以计算VAR定阶的M(i)统计量检验和各种信息准则:

library(MTS, quietly = TRUE)
Z <- coredata(as.xts(ts.gdp3r))
m3.gdp3r <- VARorder(Z/100)
## selected order: aic =  2 
## selected order: bic =  1 
## selected order: hq =  1 
## Summary table:  
##        p      AIC      BIC       HQ     M(p) p-value
##  [1,]  0 -30.9560 -30.9560 -30.9560   0.0000  0.0000
##  [2,]  1 -31.8830 -31.6794 -31.8003 115.1329  0.0000
##  [3,]  2 -31.9643 -31.5570 -31.7988  23.5389  0.0051
##  [4,]  3 -31.9236 -31.3127 -31.6754  10.4864  0.3126
##  [5,]  4 -31.8971 -31.0826 -31.5662  11.5767  0.2382
##  [6,]  5 -31.7818 -30.7636 -31.3682   2.7406  0.9737
##  [7,]  6 -31.7112 -30.4893 -31.2148   6.7822  0.6598
##  [8,]  7 -31.6180 -30.1925 -31.0389   4.5469  0.8719
##  [9,]  8 -31.7570 -30.1279 -31.0952  24.4833  0.0036
## [10,]  9 -31.6897 -29.8569 -30.9451   6.4007  0.6992
## [11,] 10 -31.5994 -29.5630 -30.7721   4.3226  0.8889
## [12,] 11 -31.6036 -29.3636 -30.6936  11.4922  0.2435
## [13,] 12 -31.6183 -29.1746 -30.6255  11.8168  0.2238
## [14,] 13 -31.6718 -29.0245 -30.5964  14.1266  0.1179

从AIC比较来看,应该取p=2
从检验来看,
i=3阶开始Φi就不显著了,
Φ1Φ2显著,
所以应该取p=2阶。

○○○○○

还可以用vars包的VAR()函数估计。
估计VAR(1):

da1 <- coredata(as.xts(ts.gdp3r))
var1 <- vars::VAR(da1, p=1)
summary(var1)
## 
## VAR Estimation Results:
## ========================= 
## Endogenous variables: uk, ca, us 
## Deterministic variables: const 
## Sample size: 124 
## Log Likelihood: -304.407 
## Roots of the characteristic polynomial:
## 0.7091 0.08735 0.05004
## Call:
## vars::VAR(y = da1, p = 1)
## 
## 
## Estimation results for equation uk: 
## =================================== 
## uk = uk.l1 + ca.l1 + us.l1 + const 
## 
##       Estimate Std. Error t value Pr(>|t|)    
## uk.l1  0.43435    0.08106   5.358 4.12e-07 ***
## ca.l1  0.18888    0.08275   2.282   0.0242 *  
## us.l1  0.03727    0.08716   0.428   0.6697    
## const  0.17133    0.06790   2.523   0.0129 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 0.5468 on 120 degrees of freedom
## Multiple R-Squared: 0.3687,  Adjusted R-squared: 0.3529 
## F-statistic: 23.36 on 3 and 120 DF,  p-value: 5.596e-12 
## 
## 
## Estimation results for equation ca: 
## =================================== 
## ca = uk.l1 + ca.l1 + us.l1 + const 
## 
##       Estimate Std. Error t value Pr(>|t|)    
## uk.l1  0.18499    0.08587   2.154   0.0332 *  
## ca.l1  0.24475    0.08766   2.792   0.0061 ** 
## us.l1  0.39166    0.09233   4.242 4.38e-05 ***
## const  0.11829    0.07193   1.644   0.1027    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 0.5792 on 120 degrees of freedom
## Multiple R-Squared: 0.4685,  Adjusted R-squared: 0.4552 
## F-statistic: 35.26 on 3 and 120 DF,  p-value: < 2.2e-16 
## 
## 
## Estimation results for equation us: 
## =================================== 
## us = uk.l1 + ca.l1 + us.l1 + const 
## 
##       Estimate Std. Error t value Pr(>|t|)    
## uk.l1  0.32153    0.09404   3.419 0.000859 ***
## ca.l1  0.18196    0.09600   1.895 0.060438 .  
## us.l1  0.16740    0.10111   1.656 0.100410    
## const  0.27859    0.07877   3.537 0.000577 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 0.6343 on 120 degrees of freedom
## Multiple R-Squared: 0.3044,  Adjusted R-squared: 0.287 
## F-statistic:  17.5 on 3 and 120 DF,  p-value: 1.725e-09 
## 
## 
## 
## Covariance matrix of residuals:
##         uk      ca      us
## uk 0.29898 0.02031 0.06841
## ca 0.02031 0.33552 0.17425
## us 0.06841 0.17425 0.40237
## 
## Correlation matrix of residuals:
##         uk      ca     us
## uk 1.00000 0.06413 0.1972
## ca 0.06413 1.00000 0.4742
## us 0.19722 0.47424 1.0000

MTS::VAR()的输出格式不同,
vars::VAR()对每个分量单独输出其方程,
好处是用我们熟悉的格式显示估计值、标准误差、检验统计量和p值。
得到的模型与MTS::VAR()结果相同,
白噪声方差估计略有差别。

用AIC自动选择阶,指定ic="AIC"lag.max=5:

da1 <- coredata(as.xts(ts.gdp3r))
var2 <- vars::VAR(da1, ic="AIC", lag.max=5)
summary(var2)
## 
## VAR Estimation Results:
## ========================= 
## Endogenous variables: uk, ca, us 
## Deterministic variables: const 
## Sample size: 121 
## Log Likelihood: -252.647 
## Roots of the characteristic polynomial:
## 0.785 0.7516 0.7516 0.7336 0.7336 0.6144 0.6144 0.5679 0.5679 0.5165 0.5122 0.5122
## Call:
## vars::VAR(y = da1, lag.max = 5, ic = "AIC")
## 
## 
## Estimation results for equation uk: 
## =================================== 
## uk = uk.l1 + ca.l1 + us.l1 + uk.l2 + ca.l2 + us.l2 + uk.l3 + ca.l3 + us.l3 + uk.l4 + ca.l4 + us.l4 + const 
## 
##        Estimate Std. Error t value Pr(>|t|)    
## uk.l1  0.515685   0.095298   5.411  3.8e-07 ***
## ca.l1  0.071863   0.094459   0.761  0.44844    
## us.l1  0.063904   0.088668   0.721  0.47264    
## uk.l2 -0.050370   0.101228  -0.498  0.61979    
## ca.l2  0.160043   0.098620   1.623  0.10754    
## us.l2 -0.001984   0.095495  -0.021  0.98346    
## uk.l3  0.052363   0.102949   0.509  0.61205    
## ca.l3 -0.278830   0.098341  -2.835  0.00547 ** 
## us.l3  0.141145   0.093728   1.506  0.13501    
## uk.l4  0.040130   0.091050   0.441  0.66028    
## ca.l4  0.261731   0.086947   3.010  0.00325 ** 
## us.l4 -0.246485   0.088234  -2.794  0.00617 ** 
## const  0.147957   0.074786   1.978  0.05043 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 0.5013 on 108 degrees of freedom
## Multiple R-Squared: 0.478,   Adjusted R-squared:  0.42 
## F-statistic: 8.242 on 12 and 108 DF,  p-value: 7.84e-11 
## 
## 
## Estimation results for equation ca: 
## =================================== 
## ca = uk.l1 + ca.l1 + us.l1 + uk.l2 + ca.l2 + us.l2 + uk.l3 + ca.l3 + us.l3 + uk.l4 + ca.l4 + us.l4 + const 
## 
##       Estimate Std. Error t value Pr(>|t|)    
## uk.l1  0.37782    0.10199   3.705 0.000336 ***
## ca.l1  0.31603    0.10109   3.126 0.002276 ** 
## us.l1  0.40963    0.09489   4.317 3.52e-05 ***
## uk.l2 -0.17399    0.10834  -1.606 0.111193    
## ca.l2 -0.25407    0.10554  -2.407 0.017771 *  
## us.l2  0.06295    0.10220   0.616 0.539247    
## uk.l3  0.09615    0.11018   0.873 0.384757    
## ca.l3  0.12035    0.10525   1.143 0.255362    
## us.l3  0.01366    0.10031   0.136 0.891925    
## uk.l4  0.07470    0.09744   0.767 0.445006    
## ca.l4 -0.09031    0.09305  -0.971 0.333959    
## us.l4 -0.09781    0.09443  -1.036 0.302622    
## const  0.07757    0.08004   0.969 0.334593    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 0.5365 on 108 degrees of freedom
## Multiple R-Squared: 0.5661,  Adjusted R-squared: 0.5179 
## F-statistic: 11.74 on 12 and 108 DF,  p-value: 8.189e-15 
## 
## 
## Estimation results for equation us: 
## =================================== 
## us = uk.l1 + ca.l1 + us.l1 + uk.l2 + ca.l2 + us.l2 + uk.l3 + ca.l3 + us.l3 + uk.l4 + ca.l4 + us.l4 + const 
## 
##       Estimate Std. Error t value Pr(>|t|)    
## uk.l1  0.51905    0.10938   4.745 6.41e-06 ***
## ca.l1  0.17304    0.10841   1.596   0.1134    
## us.l1  0.15039    0.10177   1.478   0.1424    
## uk.l2 -0.21784    0.11618  -1.875   0.0635 .  
## ca.l2 -0.15931    0.11319  -1.407   0.1622    
## us.l2  0.22561    0.10960   2.058   0.0420 *  
## uk.l3  0.04783    0.11816   0.405   0.6864    
## ca.l3 -0.07856    0.11287  -0.696   0.4879    
## us.l3  0.07384    0.10758   0.686   0.4939    
## uk.l4  0.15409    0.10450   1.475   0.1433    
## ca.l4 -0.15182    0.09979  -1.521   0.1311    
## us.l4 -0.05350    0.10127  -0.528   0.5984    
## const  0.23868    0.08584   2.781   0.0064 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 0.5754 on 108 degrees of freedom
## Multiple R-Squared: 0.4531,  Adjusted R-squared: 0.3923 
## F-statistic: 7.457 on 12 and 108 DF,  p-value: 7.522e-10 
## 
## 
## 
## Covariance matrix of residuals:
##         uk      ca     us
## uk 0.25130 0.04912 0.1001
## ca 0.04912 0.28783 0.1084
## us 0.10009 0.10840 0.3310
## 
## Correlation matrix of residuals:
##        uk     ca     us
## uk 1.0000 0.1826 0.3470
## ca 0.1826 1.0000 0.3512
## us 0.3470 0.3512 1.0000

结果选择了4阶模型,
MTS::VARorder()结果不同。

vars包的SVAR()函数可以用来估计VAR的结构形式,
但其定义不同于我们教材的定义。

○○○○○

24.8 模型检验

可以计算模型残差,
对残差进行多元白噪声检验(多元混成检验)。
残差的多元混成检验因为使用了估计的参数,
所以统计量的自由度会减少k2p
这是系数矩阵Φj, j=1,2,,p中的参数个数。
如果系数矩阵中某些参数固定为0,
应按无约束的参数个数计算要扣除的自由度。
在MTS包的mq()函数中用adj=指定需要减少的自由度。

例3.2
对例3.1的三个国家的GDP季度增速建模,
VAR(2)模型的残差的多元混成检验程序如下:

resi <- m2.gdp3r$residuals
MTS::mq(resi, adj=3^2 * 2)
## Ljung-Box Statistics:  
##           m       Q(m)     df    p-value
##  [1,]   1.000     0.816  -9.000     1.00
##  [2,]   2.000     3.978   0.000     1.00
##  [3,]   3.000    16.665   9.000     0.05
##  [4,]   4.000    35.122  18.000     0.01
##  [5,]   5.000    38.189  27.000     0.07
##  [6,]   6.000    41.239  36.000     0.25
##  [7,]   7.000    47.621  45.000     0.37
##  [8,]   8.000    61.677  54.000     0.22
##  [9,]   9.000    67.366  63.000     0.33
## [10,]  10.000    76.930  72.000     0.32
## [11,]  11.000    81.567  81.000     0.46
## [12,]  12.000    93.112  90.000     0.39
## [13,]  13.000   105.327  99.000     0.31
## [14,]  14.000   116.279 108.000     0.28
## [15,]  15.000   128.974 117.000     0.21
## [16,]  16.000   134.704 126.000     0.28
## [17,]  17.000   138.552 135.000     0.40
## [18,]  18.000   146.256 144.000     0.43
## [19,]  19.000   162.418 153.000     0.29
## [20,]  20.000   171.948 162.000     0.28
## [21,]  21.000   174.913 171.000     0.40
## [22,]  22.000   182.056 180.000     0.44
## [23,]  23.000   190.276 189.000     0.46
## [24,]  24.000   202.141 198.000     0.41

北京大学金融时间序列分析讲义第24章: 向量自回归模型

检验结果只有在滞后4显著,
基本可以认为模型是充分的。

对残差还可以进行异方差等检验。

24.9 模型简化

当VAR中分量个数k较大时,
模型有许多参数,
系数矩阵中参数个数为k2p个。
如果没有先验知识要求参数非零,
可以将不显著的参数约束为零再估计。

模型简化没有公认的最优做法。
一种办法是计算参数估计值与标准误差的比值,
称为t比值,
将t比值绝对值小于1.645(相当于0.10检验水平)或者小于1.96(相当于0.05检验水平)的系数设置为零。
MTS包的refVAR()函数输入无约束的VAR建模结果,
以及thres=1.645thres=1.96这样的t比值界限,
生成设置部分系数为零的约束估计结果,
得到精简的模型。
可以通过AIC比较完全模型和约束模型。

例3.3
对三个国家的GDP对数增长率的VAR(2)模型进行化简。

library(MTS, quietly = TRUE)
Z <- coredata(as.xts(ts.gdp3r))
cat("================ Full model:\n")
## ================ Full model:
## Constant term: 
## Estimates:  0.1258163 0.1231581 0.2895581 
## Std.Error:  0.07266338 0.07382941 0.0816888 
## AR coefficient matrix 
## AR( 1 )-matrix 
##       [,1]  [,2]   [,3]
## [1,] 0.393 0.103 0.0521
## [2,] 0.351 0.338 0.4691
## [3,] 0.491 0.240 0.2356
## standard error 
##        [,1]   [,2]   [,3]
## [1,] 0.0934 0.0984 0.0911
## [2,] 0.0949 0.1000 0.0926
## [3,] 0.1050 0.1106 0.1024
## AR( 2 )-matrix 
##         [,1]   [,2]     [,3]
## [1,]  0.0566  0.106  0.01889
## [2,] -0.1914 -0.175 -0.00868
## [3,] -0.3120 -0.131  0.08531
## standard error 
##        [,1]   [,2]   [,3]
## [1,] 0.0924 0.0876 0.0938
## [2,] 0.0939 0.0890 0.0953
## [3,] 0.1038 0.0984 0.1055
##   
## Residuals cov-mtx: 
##            [,1]       [,2]       [,3]
## [1,] 0.28244420 0.02654091 0.07435286
## [2,] 0.02654091 0.29158166 0.13948786
## [3,] 0.07435286 0.13948786 0.35696571
##   
## det(SSE) =  0.02258974 
## AIC =  -3.502259 
## BIC =  -3.094982 
## HQ  =  -3.336804
cat("\n\n================ Restricted model:\n")
## 
## 
## ================ Restricted model:
mods2.gdp3r <- refVAR(mods1.gdp3r, thres=1.96)
## Constant term: 
## Estimates:  0.1628247 0 0.2827525 
## Std.Error:  0.06814101 0 0.07972864 
## AR coefficient matrix 
## AR( 1 )-matrix 
##       [,1]  [,2]  [,3]
## [1,] 0.467 0.207 0.000
## [2,] 0.334 0.270 0.496
## [3,] 0.468 0.225 0.232
## standard error 
##        [,1]   [,2]   [,3]
## [1,] 0.0790 0.0686 0.0000
## [2,] 0.0921 0.0875 0.0913
## [3,] 0.1027 0.0963 0.1023
## AR( 2 )-matrix 
##        [,1] [,2] [,3]
## [1,]  0.000    0    0
## [2,] -0.197    0    0
## [3,] -0.301    0    0
## standard error 
##        [,1] [,2] [,3]
## [1,] 0.0000    0    0
## [2,] 0.0921    0    0
## [3,] 0.1008    0    0
##   
## Residuals cov-mtx: 
##            [,1]       [,2]       [,3]
## [1,] 0.29003669 0.01803456 0.07055856
## [2,] 0.01803456 0.30802503 0.14598345
## [3,] 0.07055856 0.14598345 0.36268779
##   
## det(SSE) =  0.02494104 
## AIC =  -3.531241 
## BIC =  -3.304976 
## HQ  =  -3.439321

无约束的VAR(2)的AIC值为3.50
约束10个参数为零的VAR(2)的AIC值为3.53
所以约束模型较优。
约束模型为

rt=Var(at)=0.1600.29+0.470.330.470.210.270.2300.500.23rt1+00.200.30000000rt2+at0.290.020.070.020.310.150.070.150.36

对这个简化的模型的检验,
可以提取残差后用MTS::mq()函数检验,
这时自由度缩减个数由原来的32×2=18个,
减少到9个,
因为模型中约束了9个参数等于零。
MTS包还提供了一个MTSdiag()函数,
输入模型结果和adj=自由度缩减个数,
作残差的CCM估计表、图和残差的多元混成检验:

library(MTS, quietly = TRUE)
MTSdiag(mods2.gdp3r, adj=9)
## [1] "Covariance matrix:"
##        uk     ca     us
## uk 0.2924 0.0182 0.0711
## ca 0.0182 0.3084 0.1472
## us 0.0711 0.1472 0.3657
## CCM at lag:  0 
##        [,1]   [,2]  [,3]
## [1,] 1.0000 0.0605 0.218
## [2,] 0.0605 1.0000 0.438
## [3,] 0.2175 0.4382 1.000
## Simplified matrix: 
## CCM at lag:  1 
## . . . 
## . . . 
## . . . 
## CCM at lag:  2 
## . . . 
## . . . 
## . . . 
## CCM at lag:  3 
## . . . 
## . . . 
## . . . 
## CCM at lag:  4 
## . . - 
## . . . 
## . . . 
## CCM at lag:  5 
## . . . 
## . . . 
## . . . 
## CCM at lag:  6 
## . . . 
## . . . 
## . . . 
## CCM at lag:  7 
## . . . 
## . . . 
## . . . 
## CCM at lag:  8 
## . . . 
## . . . 
## . . . 
## CCM at lag:  9 
## . . . 
## . . . 
## . . . 
## CCM at lag:  10 
## . . . 
## . . . 
## . . . 
## CCM at lag:  11 
## . . . 
## . . . 
## . . . 
## CCM at lag:  12 
## . . . 
## . . . 
## . . . 
## CCM at lag:  13 
## . - . 
## . . . 
## . . . 
## CCM at lag:  14 
## . - . 
## . . . 
## . . . 
## CCM at lag:  15 
## . + . 
## . . . 
## . . . 
## CCM at lag:  16 
## . . . 
## . . . 
## . . . 
## CCM at lag:  17 
## . . . 
## . . . 
## . . . 
## CCM at lag:  18 
## . . . 
## . . . 
## . . . 
## CCM at lag:  19 
## . . . 
## . . + 
## . . . 
## CCM at lag:  20 
## . . . 
## . . . 
## . . . 
## CCM at lag:  21 
## . . . 
## . . . 
## . . . 
## CCM at lag:  22 
## . . . 
## . . . 
## . . . 
## CCM at lag:  23 
## . . . 
## . . . 
## . . . 
## CCM at lag:  24 
## . . . 
## . . . 
## . . .

北京大学金融时间序列分析讲义第24章: 向量自回归模型

## Hit Enter for p-value plot of individual ccm:

北京大学金融时间序列分析讲义第24章: 向量自回归模型

## Hit Enter to compute MQ-statistics: 
## 
## Ljung-Box Statistics:  
##          m       Q(m)     df    p-value
##  [1,]   1.00      1.78    0.00     1.00
##  [2,]   2.00     12.41    9.00     0.19
##  [3,]   3.00     22.60   18.00     0.21
##  [4,]   4.00     37.71   27.00     0.08
##  [5,]   5.00     41.65   36.00     0.24
##  [6,]   6.00     44.95   45.00     0.47
##  [7,]   7.00     51.50   54.00     0.57
##  [8,]   8.00     64.87   63.00     0.41
##  [9,]   9.00     72.50   72.00     0.46
## [10,]  10.00     81.58   81.00     0.46
## [11,]  11.00     86.12   90.00     0.60
## [12,]  12.00     98.08   99.00     0.51
## [13,]  13.00    112.31  108.00     0.37
## [14,]  14.00    121.89  117.00     0.36
## [15,]  15.00    134.58  126.00     0.28
## [16,]  16.00    139.16  135.00     0.39
## [17,]  17.00    145.85  144.00     0.44
## [18,]  18.00    152.56  153.00     0.49
## [19,]  19.00    165.91  162.00     0.40
## [20,]  20.00    175.22  171.00     0.40
## [21,]  21.00    180.56  180.00     0.47
## [22,]  22.00    187.40  189.00     0.52
## [23,]  23.00    193.78  198.00     0.57
## [24,]  24.00    204.65  207.00     0.53

北京大学金融时间序列分析讲义第24章: 向量自回归模型

## Hit Enter to obtain residual plots:

北京大学金融时间序列分析讲义第24章: 向量自回归模型

从残差的CCM和多元混成检验结果来看,
约束的模型是充分的。

注:蔡教授的多元金融时间序列专著中,
残差混成检验时调整的自由度是k2p还是k2p+k
没有讲得很清楚。
在§2.7.2(P.52)说明要调整的自由度是k2p
但是在§2.7.3(P.57)的例子中,
调整后的自由度为9m12
这里扣除的12个自由度加上约束的9个零参数,
是21个参数,包括了ϕ0部分。

从简化模型可以得到三国的GDP季度对数增长率服从如下的模型:

:::r1t=0.16+0.47r1,t1+0.21r2,t1+a1tr2t=0.33r1,t1+0.27r2,t1+0.50r3,t10.20r1,t2+a2tr3t=0.28+0.47r1,t1+0.23r2,t1+0.23r3,t10.30r1,t2+a3t

残差的相关阵为:

ρ̂ 0=1.000.060.220.061.000.440.220.441.00

残差的相关阵计算程序:

cor(mods2.gdp3r$residuals)
##            [,1]       [,2]      [,3]
## [1,] 1.00000000 0.06054285 0.2175489
## [2,] 0.06054285 1.00000000 0.4382489
## [3,] 0.21754885 0.43824891 1.0000000

24.10 基于VAR模型的格兰杰因果性检验

如果模型可以简化为某些代表格兰杰因果性的系数等于零,
则可以据此进行格兰杰因果性的检验。
在二元的VAR(1)模型中,
如果约束ϕ12(1)=0后的模型与无约束模型没有显著差异,
r2t不是r1t的格兰杰原因。
p阶以及k元的情形类似。

为了比较无约束与约束的模型,
使用对数似然比检验,
得到的统计量在约束参数等于零的零假设下渐近服从卡方分布。

MTS包中GrangerTest()函数执行格兰杰因果性检验,
默认第一个分量为单向的格兰杰原因,
可以用locInput=序号指定哪一个或者哪几个分量作为原因。
用基于VAR的方法检验格兰杰因果性,
局限是各分量也必须平稳,
不支持协整模型。

例3.4

检验美国的GDP季度增长率是不是英国和加拿大的增长率的单向格兰杰原因。

library(MTS, quietly = TRUE, warn.conflicts=FALSE)
Z <- coredata(as.xts(ts.gdp3r))
GrangerTest(Z, p=2, locInput=3)
## Number of targeted zero parameters:  4 
## Chi-square test for Granger Causality and p-value:  27.2262 1.789152e-05

美国是第3个分量,
零假设为

H0:ϕ31(1)=ϕ32(1)=ϕ31(2)=ϕ32(2)=0,

即预报x3t时不需要用到x1,t1, x2,t1,
x1,t2, x2,t2
结果p值小于0.05,
拒绝H0
说明美国的GDP季度增长率不是英国和加拿大的增长率的单向格兰杰原因,
即英国和加拿大也是美国的格兰杰原因。

对加拿大进行检验:

GrangerTest(Z, p=2, locInput=2)
## Number of targeted zero parameters:  4 
## Chi-square test for Granger Causality and p-value:  48.83871 6.309173e-10

检验结果显著,
加拿大不是美国、英国的单向格兰格原因,
即美国和英国是加拿大的格兰杰原因。

对英国进行检验:

GrangerTest(Z, p=2, locInput=1)
## Number of targeted zero parameters:  4 
## Chi-square test for Granger Causality and p-value:  8.948851 0.06239076 
## Constant term: 
## Estimates:  0.2104448 0.1231581 0.2895581 
## Std.Error:  0.06685632 0.07382941 0.0816888 
## AR coefficient matrix 
## AR( 1 )-matrix 
##       [,1]  [,2]  [,3]
## [1,] 0.473 0.000 0.000
## [2,] 0.351 0.338 0.469
## [3,] 0.491 0.240 0.236
## standard error 
##        [,1]  [,2]   [,3]
## [1,] 0.0899 0.000 0.0000
## [2,] 0.0949 0.100 0.0926
## [3,] 0.1050 0.111 0.1024
## AR( 2 )-matrix 
##        [,1]   [,2]     [,3]
## [1,]  0.151  0.000  0.00000
## [2,] -0.191 -0.175 -0.00868
## [3,] -0.312 -0.131  0.08531
## standard error 
##        [,1]   [,2]   [,3]
## [1,] 0.0859 0.0000 0.0000
## [2,] 0.0939 0.0890 0.0953
## [3,] 0.1038 0.0984 0.1055
##   
## Residuals cov-mtx: 
##            [,1]       [,2]       [,3]
## [1,] 0.30423344 0.02654091 0.07435286
## [2,] 0.02654091 0.29158166 0.13948786
## [3,] 0.07435286 0.13948786 0.35696571
##   
## det(SSE) =  0.02443371 
## AIC =  -3.487791 
## BIC =  -3.17102 
## HQ  =  -3.359104

综合以上三个结果,
第三个检验说明在0.05水平下加拿大和美国不是英国的格兰杰原因,
而前两个检验说明英国是美国和加拿大的格兰杰原因,
所以可以认为英国是美国和加拿大单向的格兰杰原因。
第三个检验结果中还给出了约束系数等于零的模型估计结果。

24.11 预测

若VAR(p)模型已知,
满足平稳性条件,
{at}是独立的弱平稳时间序列。
Ft表示截止到t时刻为止的rs,st所包含的信息,
E(at|Ft1)=0
基于t时刻的信息进行超前l步预测,
预测为

rt(l)=E(rt+l|Ft)

l=1

rt(1)=ϕ0+Φ1rt++Φprt+1p

l=2

rt(2)==E(rt+2|Ft)ϕ0+Φ1E(rt+1|Ft)+Φ2rt++Φprt+2p

若记

rt(l)={E(rt+l|Ft),rt+l,l>0l0

则超前l步预报可以写成

rt(l)=E(rt+l|Ft)=ϕ0+j=1pΦjrt(lj)

可见超前多步预测可以递推地计算。

对于满足平稳性条件的VAR(p)模型,
可以证明

limlrt(l)=μ=Ert

即VAR(p)的预测具有均值回归性。
均值回归速度由特征多项式|P(z)|的最接近单位圆的根的模与1的接近程度决定,
根越接近单位圆,均值回归越慢。

l步超前预测误差为

et(l)=rt+lrt(l)=rt+lE(rt+l|Ft)

为了研究预测误差方差阵,
最好利用VAR的无穷阶MA表示。

类似一元情形,
VAR(p)模型可以写成一个无穷阶MA形式

rt=μ+l=0Ψlatl(24.14)

其中Ψ0=I,
μ=[P(1)]1(假定[P(1)]1存在),
系数矩阵Ψl可以通过下式使得方程两边z的同次项相等求解:

(IΦ1zΦpzp)(I+Ψ1z+Ψz2+)=I

易见

et(l)=at+l+Ψ1at+l1++Ψl1at+1

从而

Var(et(1))=Var(et(l))==ΣΣ+j=1l1ΨjΣΨTjVar(et(l1))+Ψl1ΣΨTl1, l=2,3,

l

Var(et(l))j=0ΨjΣΨTj=Var(rt)

这一点与均值回归性质吻合。

如果模型是从数据中估计得到的,
则点预测仍利用相同的公式,
将估计参数代入即可:

rt(1)=rt(l)=ϕ̂ 0+Φ̂ 1rt++Φ̂ prt+1pϕ̂ 0+Φ̂ 1rt(l1)++Φ̂ prt(lp)

预测误差的估计还要考虑到参数估计带来的误差,
比较复杂,参见(Tsay 2014) §2.9.2。

MTS包的VARpred()函数可以从VAR的建模结果计算点预测值,
不考虑参数估计误差的预测标准误差(Standard errors of predictions),
考虑参数估计误差的预测标准误差(Root mean squared errors of predictions)。

例3.5

例3.1中建立的三个国家的GDP增速的VAR(2)模型是基于1980年第二季度到2011年第二季度的数据,
用建立的模型进行超前1到8步预测。
第一个预测对应2011年第三季度,
最后一个预测对应2013年第二季度。

library(MTS, quietly = TRUE)
VARpred(m2.gdp3r, 8)
## orig  125 
## Forecasts at origin:  125 
##          uk      ca     us
## [1,] 0.3129 0.05166 0.1660
## [2,] 0.2647 0.31687 0.4889
## [3,] 0.3143 0.48231 0.5205
## [4,] 0.3839 0.53053 0.5998
## [5,] 0.4412 0.56978 0.6297
## [6,] 0.4799 0.59478 0.6530
## [7,] 0.5068 0.60967 0.6630
## [8,] 0.5247 0.61689 0.6688
## Standard Errors of predictions:  
##        [,1]   [,2]   [,3]
## [1,] 0.5315 0.5400 0.5975
## [2,] 0.5804 0.7165 0.7077
## [3,] 0.6202 0.7672 0.7345
## [4,] 0.6484 0.7785 0.7442
## [5,] 0.6629 0.7824 0.7475
## [6,] 0.6692 0.7838 0.7484
## [7,] 0.6719 0.7842 0.7486
## [8,] 0.6729 0.7843 0.7487
## Root mean square errors of predictions:  
##        [,1]   [,2]   [,3]
## [1,] 0.5461 0.5549 0.6140
## [2,] 0.6001 0.7799 0.7499
## [3,] 0.6365 0.7879 0.7456
## [4,] 0.6601 0.7832 0.7484
## [5,] 0.6689 0.7841 0.7488
## [6,] 0.6719 0.7844 0.7488
## [7,] 0.6730 0.7844 0.7487
## [8,] 0.6734 0.7844 0.7487

多步预测会趋近到序列均值,计算序列均值:

Z <- coredata(as.xts(ts.gdp3r))
colMeans(Z)
##        uk        ca        us 
## 0.5223092 0.6153672 0.6473996

可以看出在超前8步时的点预测值很接近于序列的均值。

不管是不考虑参数估计误差的标准误差(Standard errors of predictions)
还是考虑参数估计误差的标准误差(Root mean squared errors of predictions),
超前l步预报当l时都应该趋于序列的标准差。
计算序列的样本标准差:

Z <- coredata(as.xts(ts.gdp3r))
apply(Z, 2, sd)
##        uk        ca        us 
## 0.7086442 0.7851955 0.7872912

输出中的标准误差或者根均方误差可以用来计算预测区间。
比如,
计算美国GDP季度对数增长率超前2步预测区间,
即2011年第四季度的预测区间,
在95%置信度下可计算为
0.4889±1.96×0.7077
0.4889±1.96×0.7499
后一个区间用到的标准误差包含了参数估计误差带来的额外误差估计。

24.12 脉冲响应函数

回忆一元ARMA模型的Wold表示:

Xt=j=0ψjεj,

其中{εt}是白噪声列,
称为{Xt}新息
{εt}为独立白噪声时有

εt=XtE(Xt|Xt1,Xt2,).

{ψj}称为{Xt}Wold系数
又称为{Xt}脉冲响应函数(IRF, impulse response function),
这里函数自变量为滞后j
易见

Xtεtj=ψj,

所以ψj代表了过去的输入项εtj的一个单位的变化对j期后的Xt的影响大小。
“脉冲响应”这个名词是假想取输入{εt}ε0=1
其它的εt=0时,
输入序列是一个脉冲,
这个脉冲在t>=0时刻的输出为Xt=ψt
所以{ψj,j=0,1,2,}序列是脉冲输入对应的输出。

前面给出了VAR模型的无穷阶MA表示(24.14)

rt=μ+l=0Ψlatl

Ψl是过去的信息对rt的影响的系数,
Ψl的元素为时间序列rt脉冲响应函数的系数。
等价地,
Ψl的元素也是at对未来的rt+l的线性影响的系数。

k=2的序列rt=(r1t,r2t)T为例,
Ψl=(ψij(l))k×k
考虑r1t的一个变化对r2,t+j的影响。
为此,设μ=0
at=0(t0时),
a0=(1,0)T,
at=0(t0时),
这样来考虑a1,0的变化引起的后续变化,
也可以认为是r1,0的一个变化引起的后续变化。
这时

r0=r1=r2=a0=(10)Ψ1a0=(ψ11(1)ψ21(1))Ψ2a0=(ψ11(2)ψ21(2))

所以Ψl(i,j)元可以看成是rjt的一个变动对ri,t+l造成的影响。


Ψn=l=0nΨl

Ψn表示对rt的单位冲击的累积响应
Ψn的元素为第n个短期乘数,
全部的响应的累积值

Ψ=l=0Ψl

称为总乘数或长期效应。

可以作ψij(l)沿l=0,1,2,变化的折线图,
称为脉冲响应函数图;
以及nl=0ψij(l)沿n=0,1,2,变化的折线图,
称为累积响应函数图。

例3.6

对例3.2中建立的简化的VAR(2)模型作脉冲响应函数图和累积响应图。
三个分量两两的图形共3×3×2幅。

library(MTS, quietly = TRUE)
VARMAirf(Phi = mods2.gdp3r$Phi, 
  Sigma = mods2.gdp3r$Sig, 
  orth=FALSE)

北京大学金融时间序列分析讲义第24章: 向量自回归模型

## Press return to continue

北京大学金融时间序列分析讲义第24章: 向量自回归模型

在脉冲响应的9幅图中,
右上角的图是ψ13(l)的图形,
代表美国(分量3)的数据的一个冲击对英国(分量1)的延迟的影响。

○○○○○

这样的脉冲响应有一些问题。
由于新息at存在分量间的同步相关性,
改变a1t的值不可能不同时影响到a2t的值。
所以原始的脉冲响应Ψl不够实际。
用数学解释,
改变a1t后对rt+l的影响,
可以看成是求
rt+la1t
利用无穷阶MA表示(24.14)
因为{at}序列不相关,
所以(以k=2为例)

rt+la1t=Ψlata1t=Ψl(1a2ta1t)

其中a2t不是a1t的函数,
其关系可以用线性回归

a2t=σ21σ11a1t+e

表示,
因此a2ta1t可估计为σ21σ11

rt+la1t=Ψl(1σ21σ11)

所以a1trt+l的影响不仅涉及Ψl的第一列元素的影响,
还涉及Σ的第一列元素的影响以及Ψl的第二列元素的影响。

可以利用Cholesky分解使得新息正交化(分量间不相关)。
at的方差阵Σ有Cholesky分解
Σ=LLT
其中L是对角线元素都为正数的下三角阵。
定义正交化的新息bt=L1at
Var(bt)=I是单位阵,
从而bt各分量不相关。
Ψl=ΨlL=(ψij(l))k×k
(24.14)可以用正交化的新息bt表示为

rt==μ+l=0ΨlLL1atμ+l=0Ψlbt

因为bt的分量不相关,
可得

ri,t+lbjt=ψij(l)

ψij(l)bjt对未来的观测ri,t+l的影响,
称系数矩阵Ψl
rt的带正交新息bt的脉冲响应系数,
称沿l=0,1,2,变化的ψij(l)rt
的带正交新息bt的脉冲响应函数。

正交化以后,
ψij(0)L1的第(i,j)元素,
可以表示rjt的一个冲击对rit的即期影响。

这样的三角分解的缺点是结果依赖于rt的分量次序,
对脉冲函数的理解与rt的分量次序有关。
特别地,
a1t=σ11b1t
a1t仅做了方差标准化。

例3.7

对例3.2中建立的简化的VAR(2)模型作正交化的脉冲响应函数图和累积响应图。

library(MTS, quietly = TRUE)
VARMAirf(Phi = mods2.gdp3r$Phi, 
  Sigma = mods2.gdp3r$Sig, 
  orth=TRUE)

北京大学金融时间序列分析讲义第24章: 向量自回归模型

## Press return to continue

北京大学金融时间序列分析讲义第24章: 向量自回归模型

这三个分量的新息的即期相关不大,
所以正交化的脉冲响应与原始的表现相近。

○○○○○

韭菜热线原创版权所有,发布者:风生水起,转载请注明出处:https://www.9crx.com/74791.html

(0)
打赏
风生水起的头像风生水起普通用户
上一篇 2023年8月5日 22:37
下一篇 2023年8月7日 01:35

相关推荐

  • 总回报过山车

    这是一组有趣的图表,特别会引起我们这些关注经济和市场周期的人的共鸣。想象一下,五年前,您在标准普尔 500 指数中投资了 10,000 美元。如果将股息进行再投资,但根据通货膨胀进行调整,那么今天它值多少钱? 您的投资购买力已增加至 13,653 美元,年化实际回报率为 6.24%。 如果我们在 2009 年 3 月提出同样的问题,答案将是令人沮丧的 6,6…

    2023年10月9日
    9000
  • 市场效率较低=阿尔法较高?

    介绍 学生经常向我寻求职业建议。这并不是一次特别令人满意的经历。一方面,这些人通常都是非常聪明和勤奋的人,拥有牛津或剑桥化学工程、天体物理学或其他一些具有挑战性学科的博士学位。我希望他们能够坚持科学,为我们的文明创造一些有意义的东西,而不是试图每年产生一些多余的基点。 另一方面,一些学生很早就决定从事金融职业并进行相应的学习。告诉他们建造更好的化肥或火箭飞船…

    2023年7月21日
    10400
  • 规划未来支出变化

    退休后的持续规划包括定期评估是否可以增加支出或必须减少支出以保持正常运转。德里克·塔普博士在最近的 Kitces.com 文章《如何沟通护栏提款策略可以改善客户体验并减轻压力》中说: “然而,蒙特卡洛模拟的结果通常没有考虑到在此过程中可能做出的潜在调整(例如,如果市场回报较弱且成功概率下降,则减少提款,反之亦然),这使得它们对于持续的规划工作不太有用,因为顾…

    2024年5月17日
    4200
  • 在高度分散的世界中释放阿尔法机遇

    作者: Rich Mathieson,Christopher DiPrimio,24 年 3 月 6 日 关键点 离散度上升:从长期停滞状态向通货再膨胀状态的转变,不仅导致全球盈利增长扩大,而且导致公司业绩和业绩的离散度显着提高。 阿尔法优于贝塔:除了加剧分散性之外,新制度可能会抑制近几十年来的整体市场表现。这反映出机会集有利于阿尔法回报来源而不是贝塔回报来…

    2024年3月15日
    5800
  • 北京大学R语言教程(李东风)第14章: 工作空间和变量赋值

    工作空间 R把在命令行定义的变量都保存到工作空间中,在退出R时可以选择是否保存工作空间。这也是R与其他如C、Java这样的语言的区别之一。 用ls()命令可以查看工作空间中的内容。 随着多次在命令行使用R,工作空间的变量越来越多,使得重名的可能性越来越大,而且工作空间中变量太多也让我们不容易查看其内容。在命令行定义的变量称为“全局变量”,在编程实践中,全局变…

    2023年11月2日
    12300

发表回复

登录后才能评论
客服
客服
关注订阅号
关注订阅号
分享本页
返回顶部