平稳VAR模型一

knitr::opts_chunk$set(echo = TRUE, message = F, warning = F)

Tsay书的第一章通过上两篇得到总结，主要目的是清楚多元模型的数据结构和模型目标。第二章用以了解 Stationary VAR 时间序列的基本知识，主要包括:

• 模型成分和结构
• 平稳性特征/条件
• 模型测量
• 模型筛选
• 预测

总体的，一个VAR(p)表达与AR(p)类似:

$$z_t = \phi_0 + \sum_{i = 1}^p \phi_i z_{t - i} + a_t$$

参见上篇详细表达。 值得注意的是，VAR(p) 的 p 取值于最大lag， 如一个观测对象数量为2的VAR(2)模型， 可以由一个AR(1)和一个AR(2)模型组成， 如下, 例子中$z_{2, t-1}$如果不能为z_{1t} 的预测提供有效信息, 依照Granger
Causality $\phi_{1,12} = 0$, 原因是它不能提高预测的准确性。

$$z_{1t} = \phi_{10} + \phi_{1, 11}z_{1, t-1} +\phi_{1,12}z_{2, t-1}+ a_{1t} \\ z_{2t} = \phi_{20} + \phi_{1, 21}z_{1, t-1} + \phi_{1,22}z_{2, t-1}+a_{2t}$$

另外，现实中P的值都很小，以下在尽量不失总体性的情况下以如下三元VAR(2)模型为例子。 对 VAR(2) 模型结构和成分分析主要考虑moment分析，平稳性， MA变形以及单元序列变形三方面。

$$\left(\begin{array}{c}z_{1t}\\z_{2t}\\z_{3t}\end{array}\right) = \left(\begin{array}{c}5\\3\\0\end{array}\right) + \begin{bmatrix}0.47 & 0.21 & 0 \\-0.35 & 0.34& 0.47\\0.47& 0.23& 0.23 \end{bmatrix} \left(\begin{array}{c}z_{1,t-1}\\z_{2,t-1}\\z_{3,t-1}\end{array}\right)\\ + \begin{bmatrix}0 & 0 & 0 \\-0.19 & 0.18& 0\\0.3& 0& 0 \end{bmatrix} \left(\begin{array}{c}z_{1,t-1}\\z_{2,t-1}\\z_{3,t-1}\end{array}\right) + \left(\begin{array}{c}a_{1t}\\a_{2t}\\a_{3t}\end{array}\right) ,\\{\rm with }\quad\Sigma_a= \begin{bmatrix}0.285 & 0.026 & 0.069 \\0.026 & 0.287& 0.137\\0.069& 0.137& 0.357 \end{bmatrix}.$$

本篇仅考虑前两个Monment，即均值和协方差，公式和计算过程如下。

• 均值

$$\mu = (I_k - \phi_1 -\phi_2)^{-1}\phi_0$$

• 方差和前p-1个lag 的协方差矩阵:

$${\rm vec}(\Gamma_0^{*}) = [I_{(kp)^2} - \Phi\otimes\Phi]^{-1}{\rm vec}(\Gamma_b) \qquad {\rm with:}\\ \Phi = \begin{bmatrix}\phi_1 & \phi_2 &\cdots & \phi_{p-1}& \phi_{p}\\I & 0 &\cdots& 0& 0\\0&I&\cdots& 0& 0\\\vdots&\vdots&\ddots&\vdots&\vdots&\\0& &\cdots &I&0\end{bmatrix},\\ {\rm and}\qquad \Gamma_b =\begin{bmatrix}\Gamma_a & 0_k \\0_k & 0_k \end{bmatrix} \\ {\rm generally, \quad} \Gamma_\ell = \sum_{i = 1}^p\phi_i\Gamma_{\ell-i} \quad\forall \ell > 0$$

使用上述VAR(2)数值进行如下计算

library(MTS)
p0 = matrix(c(5, 3, 0), 3, 1)
p1 = matrix(c(.47, -.35, .47, .21, .34, .23, 0, .47, .23), 3, 3)
p2 = matrix(c(0,-.19, .3, 0, .18, 0, 0,0,0),3,3)
sig = matrix(c(.285, .026, .069,.026,.287,.137,.069,.137,.357),3,3)
k =  nrow(sig)
Ik= diag(k)
mu = solve(Ik - p1 - p2) %*% p0
rownames(mu) = c("z1","z2","z3")
colnames(mu) = c(mean)
mu
zt = VARMAsim(30000, arlags = c(1,2), phi = cbind(p1,p2), cnst = c(p0),
              sigma = sig, set.seed(100))
meansm = matrix(apply(zt$series, 2, mean), nrow =1) 
"mean throung simulation"
meansm
P1 = diag(ncol(p1))
P2 = P1 * 0
PR1 = cbind(p1, p2)
PR2 = cbind(P1, P2)
P = rbind(PR1, PR2)
P
G1 = matrix(0, nrow(sig),ncol(sig))
G2 = G1
G3 = G1
GR1 = cbind(sig, G1)
GR2 = cbind(G2, G3)
G = rbind(GR1, GR2)
G
k = 3
p = 2
nVI = (k*p)^2
VI = diag(nVI)
VKP = kronecker(P, P)
G0 = solve(VI - VKP) %*% matrix(as.numeric(G), ncol = 1)
nG = length(G)^0.5
GM = matrix(G0, nrow = nG)
g0 = GM[1:nrow(sig), 1:ncol(sig)]
g1 = GM[1:nrow(sig), (ncol(g0)+1):ncol(GM)]
g2 = p1 %*% g1 + p2 %*% g0
GamaM = cbind(g0, g1, g2)
rownames(GamaM) = c("z1","z2","z3")
colnames(GamaM) = paste("L", rep(c(1,2,3),  3), rep(c(0,1,2),  each =3), sep = "")
GamaM
sm = diag(diag(g0))
se = sm^0.5
s = solve(se)
rho0 = s %*% g0  %*% s
rho1 = s %*% g1  %*% s
rho2 = s %*% g2  %*% s
RohM = cbind(rho0, rho1, rho2)
dimnames(RohM) = dimnames(GamaM)
RohM
zt = zt$series
nr = nrow(zt)
zt1 = zt[1: (nr - 2),]
zt2 = zt[2:(nr -1),]
zt3 = zt[3:nr,]
Rohsim = cbind(cor(zt), cor(zt2,zt1), cor(zt3, zt1))
"check throung simulation, outcomes differ somehow from the parameters's one if T is small"
Rohsim 

Quelle

Rmarkdwon版本
Tsay, Ruey S(2014): *Multivariate Time Series Analysis
With R and Financial Applications*, Chapter 2.