平稳VAR模型一
Chang
2017年11月17日
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变形以及单元序列变形, 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## function (x, ...) \nUseMethod("mean")
## z1 11.957522
## z2 6.368985
## z3 13.859946zt = 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"## [1] "mean throung simulation"meansm## [,1] [,2] [,3]
## [1,] 11.96229 6.363558 13.86553P1 = diag(ncol(p1))
P2 = P1 * 0
PR1 = cbind(p1, p2)
PR2 = cbind(P1, P2)
P = rbind(PR1, PR2)
P## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 0.47 0.21 0.00 0.00 0.00 0
## [2,] -0.35 0.34 0.47 -0.19 0.18 0
## [3,] 0.47 0.23 0.23 0.30 0.00 0
## [4,] 1.00 0.00 0.00 0.00 0.00 0
## [5,] 0.00 1.00 0.00 0.00 0.00 0
## [6,] 0.00 0.00 1.00 0.00 0.00 0G1 = matrix(0, nrow(sig),ncol(sig))
G2 = G1
G3 = G1
GR1 = cbind(sig, G1)
GR2 = cbind(G2, G3)
G = rbind(GR1, GR2)
G## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 0.285 0.026 0.069 0 0 0
## [2,] 0.026 0.287 0.137 0 0 0
## [3,] 0.069 0.137 0.357 0 0 0
## [4,] 0.000 0.000 0.000 0 0 0
## [5,] 0.000 0.000 0.000 0 0 0
## [6,] 0.000 0.000 0.000 0 0 0k = 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## L10 L20 L30 L11 L21 L31 L12
## z1 0.4324548 0.1209413 0.3030421 0.22865141 0.1904219 0.2134289 0.11230096
## z2 0.1209413 0.6360930 0.3380911 0.02302285 0.4005740 0.3732044 0.04100451
## z3 0.3030421 0.3380911 0.8142372 0.36936535 0.2878116 0.5182749 0.32745188
## L22 L32
## z1 0.1736188 0.1786845
## z2 0.2963368 0.2990570
## z3 0.2841094 0.3962645sm = 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## L10 L20 L30 L11 L21 L31 L12
## z1 1.0000000 0.2305918 0.5106899 0.52872907 0.3630666 0.3596727 0.25968256
## z2 0.2305918 1.0000000 0.4697838 0.04389636 0.6297412 0.5185743 0.07818097
## z3 0.5106899 0.4697838 1.0000000 0.62245842 0.3999194 0.6365159 0.55182539
## L22 L32
## z1 0.3310291 0.3011211
## z2 0.4658703 0.4155451
## z3 0.3947751 0.4866696zt = 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 differs somehow from the parameters's one if T is small"## [1] "check throung simulation, outcomes differs somehow from the parameters's one if T is small"Rohsim ## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 1.0000000 0.2388284 0.5225252 0.53972625 0.3668379 0.3730606
## [2,] 0.2388284 1.0000000 0.4821114 0.05368308 0.6451983 0.5298558
## [3,] 0.5225252 0.4821114 1.0000000 0.62714000 0.4144837 0.6497877
## [,7] [,8] [,9]
## [1,] 0.27232618 0.3355396 0.3131459
## [2,] 0.08489265 0.4779336 0.4336809
## [3,] 0.55778361 0.4065202 0.5027382每一个VAR(p)都可以转化成一个VAR(1)模型，VAR(p)的平稳性分析都可以转换成VAR(1)的平稳性分析， 同时每一个VAR(1)中的中的单独序列都可以转成单元的ARMA序列，于是逻辑上对于每一个VAR(P)模型也是适用。总体上，转换为VAR(1)使用一下公式:
\[ \begin{bmatrix} z_t = \begin{bmatrix}z_{t,10}& \cdots&z_{t,k0}\\\vdots & \vdots& \vdots \end{bmatrix} '\\ z_{t-1}=\begin{bmatrix}z_{t-1,10}& \cdots&z_{t-1,k0}\\\vdots & \vdots& \vdots \end{bmatrix} ' \end{bmatrix} = \Phi_0 + \Phi_1\begin{bmatrix} z_{t-1}\begin{bmatrix}z_{t-1,10}& \cdots&z_{t-1,k0}\\\vdots & \vdots& \vdots \end{bmatrix}' \\ z_{t-2}\begin{bmatrix}z_{t-2,10}& \cdots&z_{t-2,k0}\\\vdots & \vdots& \vdots \end{bmatrix} ' \end{bmatrix} +b_t \quad{\rm with:} \\ \Phi_0 = \begin{bmatrix}\phi_0 \\0_{k*1} \end{bmatrix}, \quad b_t = \begin{bmatrix}a_t \\0_{k*1} \end{bmatrix} \qquad {\rm and} \quad \Phi_1 = \Phi {\rm \quad in \quad Moments \quad Equaition}\]
于是通过系数判断模型是否平稳，取决于\(\Phi_1\) 的 Eingenvalue，当说有Eingenvalue的Modulus绝对值小于一时候模型平稳， 此判断与characteristic roots的Modulus绝对值大于1的要求一致。上面例子的VAR(2)通过计算表面是平稳的， 模拟出来的图像也很明显验证了这一点。
require(MTS)
ev = eigen(P)$values
Mod(ev)## [1] 8.092769e-01 4.304480e-01 4.304480e-01 3.884366e-01 2.553013e-17
## [6] 0.000000e+00abs(ev) < 1## [1] TRUE TRUE TRUE TRUE TRUE TRUEMTSplot(zt)
下面开始进入模型测量。常用方法不外乎 各种LS MLE 和 bayesian 测量。首先看LS Estimate, 公式如普通回归 \[ \hat \beta = (X'X)^{-1}X'Y \\ \hat \sigma_a = \frac{1}{T-(k+1)p -1}\hat A'\hat A, \quad A = Y-X\hat\beta\\ {\rm Cov[ vec} (\hat \beta)] = \widehat{\Sigma_a}\otimes(X'X)^{-1} \]
X = cbind(rep(1, nrow(zt3)), zt2, zt1)
phi_lse = solve(t(X) %*% X) %*% t(X) %*% zt3
rownames(phi_lse) = paste("phi", rep(c(0, 1, 2), c(1,3,3)),
c(0,1:3,1:3) ,sep = "")
colnames(phi_lse) = rownames(mu)
A = zt3 - X %*% phi_lse
sig_lse = crossprod(A, A)/ (nrow(zt) -(1 + ncol(phi_lse))*2 -1)
COV = kronecker(sig_lse, solve(t(X) %*% X))
se = sqrt(diag(COV))
para = cbind(c(phi_lse), se, c(phi_lse)/se )
rownames(para) = paste(rep(colnames(phi_lse), each= 7), rep(rownames(phi_lse), 3), sep = "")
colnames(para) = c("coef", "se", "t-ratio")
coef_given = as.numeric(rbind(t(p0), t(p1), t(p2)))
para = cbind(para, coef_given)
para#LS Estimation## coef se t-ratio coef_given
## z1phi00 4.9174261981 0.066686999 73.73890389 5.00
## z1phi11 0.4802149627 0.005915176 81.18355271 0.47
## z1phi12 0.2060405541 0.006026628 34.18836232 0.21
## z1phi13 0.0008996660 0.005561312 0.16177227 0.00
## z1phi21 -0.0007231090 0.007330365 -0.09864571 0.00
## z1phi22 -0.0047889678 0.005237393 -0.91437999 0.00
## z1phi23 0.0011486804 0.005066705 0.22671150 0.00
## z2phi00 2.9039737957 0.067366546 43.10706077 3.00
## z2phi11 -0.3580733104 0.005975452 -59.92405691 -0.35
## z2phi12 0.3501726948 0.006088040 57.51813007 0.34
## z2phi13 0.4807659822 0.005617982 85.57627845 0.47
## z2phi21 -0.1962219452 0.007405062 -26.49835389 -0.19
## z2phi22 0.1697023011 0.005290763 32.07520343 0.18
## z2phi23 0.0083603348 0.005118336 1.63340882 0.00
## z3phi00 -0.0057094266 0.074534389 -0.07660124 0.00
## z3phi11 0.4638615620 0.006611244 70.16252821 0.47
## z3phi12 0.2325694355 0.006735812 34.52730699 0.23
## z3phi13 0.2393347602 0.006215739 38.50463683 0.23
## z3phi21 0.2965456056 0.008192965 36.19515153 0.30
## z3phi22 -0.0023146123 0.005853704 -0.39540990 0.00
## z3phi23 -0.0006267898 0.005662930 -0.11068295 0.00sigv = cbind(c(sig_lse), c(sig))
colnames(sigv) = c("sig_est", "sig_given")
sigv#comparison of sigma## sig_est sig_given
## [1,] 0.28262655 0.285
## [2,] 0.02707394 0.026
## [3,] 0.07018320 0.069
## [4,] 0.02707394 0.026
## [5,] 0.28841588 0.287
## [6,] 0.13740144 0.137
## [7,] 0.07018320 0.069
## [8,] 0.13740144 0.137
## [9,] 0.35305633 0.357结果现实LS所测量出来的系数跟“捏造”的系数很接近，尽管有些不在95的confidence level下有部分不显著，残差的的协方差矩阵也几乎一致。以为下使用lm函数或者MTS包中的VAR函数得到的结果。
require(MTS)
lm(zt3 ~ zt2 + zt1 )##
## Call:
## lm(formula = zt3 ~ zt2 + zt1)
##
## Coefficients:
## [,1] [,2] [,3]
## (Intercept) 4.9174262 2.9039738 -0.0057094
## zt21 0.4802150 -0.3580733 0.4638616
## zt22 0.2060406 0.3501727 0.2325694
## zt23 0.0008997 0.4807660 0.2393348
## zt11 -0.0007231 -0.1962219 0.2965456
## zt12 -0.0047890 0.1697023 -0.0023146
## zt13 0.0011487 0.0083603 -0.0006268VAR(zt, 2)## Constant term:
## Estimates: 4.917426 2.903974 -0.005709427
## Std.Error: 0.066687 0.06736655 0.07453439
## AR coefficient matrix
## AR( 1 )-matrix
## [,1] [,2] [,3]
## [1,] 0.480 0.206 0.0009
## [2,] -0.358 0.350 0.4808
## [3,] 0.464 0.233 0.2393
## standard error
## [,1] [,2] [,3]
## [1,] 0.00592 0.00603 0.00556
## [2,] 0.00598 0.00609 0.00562
## [3,] 0.00661 0.00674 0.00622
## AR( 2 )-matrix
## [,1] [,2] [,3]
## [1,] -0.000723 -0.00479 0.001149
## [2,] -0.196222 0.16970 0.008360
## [3,] 0.296546 -0.00231 -0.000627
## standard error
## [,1] [,2] [,3]
## [1,] 0.00733 0.00524 0.00507
## [2,] 0.00741 0.00529 0.00512
## [3,] 0.00819 0.00585 0.00566
##
## Residuals cov-mtx:
## [,1] [,2] [,3]
## [1,] 0.28256060 0.02706762 0.07016683
## [2,] 0.02706762 0.28834858 0.13736937
## [3,] 0.07016683 0.13736937 0.35297394
##
## det(SSE) = 0.02227041
## AIC = -3.803296
## BIC = -3.798311
## HQ = -3.801697MLE 将在 Vola Model 中 使用， 此处略过。 接着简单介绍下Bayesian Estimator，详细介绍在State Space Model中。Bayesian Estimator 更多的是种概念而不是方法，测量算法于LS一样，不同地，Baysian假定先验地已经知道了残差服从inverted Wishart Distribution \(\Sigma_a \sim W^{-1}(V_0, n_0)\), 根据定义已知道其期望为0， Covriance \(V_0\) 和自由度 \(n_0\) 需要假设已经知，其中为保证\(V_0\)为正，自由度 \(n_0 > K + 1\), K 为序列观测量的个数。假定系数的服从正态分 \({\rm vec}(\beta) \sim N({\rm vec}(\beta_0), \quad \Sigma_a \otimes C^{-1})\) 先验值的选择很主观，有许多成俗的做法， 简单地，对于VAR, Minnesota Prior 常被使用:
\[ V_0 = I_k, \qquad n_0 = k + 2 \\ \beta_0 = 0_{(kp+1)*k}\\ \Sigma_a \otimes C^{-1} = C = \lambda I_{kp+1}, \quad \lambda{\rm \quad is \quad small,\quad such \quad as\quad 0.1 } \]
总体上， Bayesian Estimatior是寻找LS Estimator 和 先验值的均衡。 系数和残差测量公式根据以上，结合LS Estimator为如下 \[ \hat \beta = (X'X + C)^{-1} (X'Y + C\beta_0)\\ A = Y - \hat Y\\\ B = \hat \beta - \beta_0\\ S = A'A + B'CB\\ \Sigma_a = \frac{V_0 + S}{n_0 + T - p - k -1}\\ {\rm Cov}(\hat\beta) = \Sigma_a \otimes (X'X + C)^{-1} \] 仍让以上面的VAR(2)为例子计算如下
phi0 = matrix(0, ncol(zt)*p+1, ncol(zt))
lmd = 0.1
bcv = lmd*diag(nrow(phi0))
acv = diag(ncol(phi0))
x = cbind(rep(1, nrow(zt3)), zt2, zt1)
y = zt3
phi = solve(crossprod(x, x) + bcv) %*% (crossprod(x, y) + bcv %*% phi0)
A = zt3 - x %*% phi
B = phi - phi0
S = crossprod(A, A) + t(B) %*% bcv %*% B
byssig = (acv + S) / ((k+2) + nrow(zt3) - ncol(phi0) -1)
phicv = kronecker(byssig, solve(crossprod(x, x) + bcv))
se = diag(phicv)^0.5
phivec = c(phi)
bysest = cbind(phivec, se, phivec/se)
bysest = cbind(bysest, para[, 4])
dimnames(bysest) = dimnames(para[,1:4])
bysest # bayesian estimate## coef se t-ratio coef_given
## z1phi00 4.9097329625 0.066639071 73.67649166 5.00
## z1phi11 0.4804942018 0.005914738 81.23676725 0.47
## z1phi12 0.2061400042 0.006026881 34.20343154 0.21
## z1phi13 0.0008623303 0.005561633 0.15504983 0.00
## z1phi21 -0.0004531627 0.007330200 -0.06182133 0.00
## z1phi22 -0.0048735546 0.005237637 -0.93048725 0.00
## z1phi23 0.0012590570 0.005066884 0.24848742 0.00
## z2phi00 2.8993806728 0.067311847 43.07385418 3.00
## z2phi11 -0.3578998999 0.005974452 -59.90505710 -0.35
## z2phi12 0.3502329507 0.006087727 57.53098776 0.34
## z2phi13 0.4807367817 0.005617783 85.57411733 0.47
## z2phi21 -0.1960600567 0.007404204 -26.47955804 -0.19
## z2phi22 0.1696498753 0.005290515 32.06679544 0.18
## z2phi23 0.0084272314 0.005118039 1.64657442 0.00
## z3phi00 -0.0056543874 0.074469367 -0.07592904 0.00
## z3phi11 0.4638574648 0.006609738 70.17788973 0.47
## z3phi12 0.2325652909 0.006735058 34.53055603 0.23
## z3phi13 0.2393362106 0.006215142 38.50856474 0.23
## z3phi21 0.2965417230 0.008191521 36.20105772 0.30
## z3phi22 -0.0023119666 0.005853076 -0.39500030 0.00
## z3phi23 -0.0006246360 0.005662259 -0.11031569 0.00sigv[,1] = c(byssig)
sigv # comparsion sigma## sig_est sig_given
## [1,] 0.28266591 0.285
## [2,] 0.02711391 0.026
## [3,] 0.07016530 0.069
## [4,] 0.02711391 0.026
## [5,] 0.28840220 0.287
## [6,] 0.13736465 0.137
## [7,] 0.07016530 0.069
## [8,] 0.13736465 0.137
## [9,] 0.35299689 0.357预测结果和捏造的系数以及残差方差矩阵几乎一致， 系数也大多显著，不显著的都很小，本质上就很难与零区别。同样，测量也可以用MTS中的BVAR函数获得
library(MTS)
BVAR(zt, p = 2, C = bcv, V = acv)## Bayesian estimate:
## Est s.e. t-ratio
## [1,] 4.9097329625 0.066639071 73.67649166
## [2,] 0.4804942018 0.005914738 81.23676725
## [3,] 0.2061400042 0.006026881 34.20343154
## [4,] 0.0008623303 0.005561633 0.15504983
## [5,] -0.0004531627 0.007330200 -0.06182133
## [6,] -0.0048735546 0.005237637 -0.93048725
## [7,] 0.0012590570 0.005066884 0.24848742
## [8,] 2.8993806728 0.067311847 43.07385418
## [9,] -0.3578998999 0.005974452 -59.90505710
## [10,] 0.3502329507 0.006087727 57.53098776
## [11,] 0.4807367817 0.005617783 85.57411733
## [12,] -0.1960600567 0.007404204 -26.47955804
## [13,] 0.1696498753 0.005290515 32.06679544
## [14,] 0.0084272314 0.005118039 1.64657442
## [15,] -0.0056543874 0.074469367 -0.07592904
## [16,] 0.4638574648 0.006609738 70.17788973
## [17,] 0.2325652909 0.006735058 34.53055603
## [18,] 0.2393362106 0.006215142 38.50856474
## [19,] 0.2965417230 0.008191521 36.20105772
## [20,] -0.0023119666 0.005853076 -0.39500030
## [21,] -0.0006246360 0.005662259 -0.11031569
## Covariance matrix:
## [,1] [,2] [,3]
## [1,] 0.28266591 0.02711391 0.0701653
## [2,] 0.02711391 0.28840220 0.1373646
## [3,] 0.07016530 0.13736465 0.3529969但是肉眼比较哪个Estimator或者哪个模型更好，比较不容易。 好的常用标准类似univarite analysis ，In Sample用Information Criteria， Out of Sample可以使用比如MSE。 首先比较In Sample， 第一个方法是 sequential likelihood test, 使用下面的公式寻找 P 值\(P_\ell\)， 如果 \(P_\ell > \alpha\) ,则lag \(p = \ell - 1\) \[ P_\ell = 1 - \chi^2(M(\ell), k^2), \quad {\rm with:}\\ M(\ell) = -(T -\ell -1.5 -k\ell){\ln}\left(\begin{array}{c}\frac{\mid\hat\Sigma_{a,\ell}\ \mid}{\mid\hat\Sigma_{a,\ell}\mid}\end{array}\right)\\ \]
Information Criteria 的公式和单元分析的一样， 选择 lag p 根据最小的信息值。 \[ {\rm AIC}(\ell) = \ln \mid\hat\Sigma_{a,\ell}\mid + \frac{2}{T}\ell k^2\\ {\rm BIC}(\ell) = \ln \mid\hat\Sigma_{a,\ell} \mid + \frac{\ln(T)}{T}\ell k^2\\ {\rm HQ}(\ell) = \ln \mid\hat\Sigma_{a,\ell} \mid + \frac{2\ln[\ln(T)]}{T}\ell k^2 \] 下面使用例子中的数据，使用LS Estimator， 来对比 lag p = 1 和 p = 2。
p = 0# lag P = 1
sig0a = cov(zt)*(nrow(zt)-1)/nrow(zt)
m0 = -(nrow(x) - p - 1.5-ncol(zt)*p)*log(det(sig0a)/1)
pvp0 = 1 - pchisq(m0, ncol(zt)^2)
aci0 = log(det(sig0a)) + 2/nrow(zt)*p*ncol(zt)^2
bic0 = log(det(sig0a)) + log(nrow(zt))/nrow(zt)*p*(ncol(zt))^2
hq0 = log(det(sig0a)) + 2*(log(log(nrow(zt)))) /nrow(zt)*p*ncol(zt)^2
p = 1# lag P = 1
zp1y = zt[(p+1):nrow(zt),]
zp1x = zt[1: (nrow(zt)-p),]
x = cbind(rep(1, nrow(zp1x)), zp1x)
y = zp1y
b1 = solve(crossprod(x, x), crossprod(x,y))
a1 = y - x %*% b1
sig1a = crossprod(a1, a1) / (nrow(zt)-p)
m1 = -(nrow(x) - p - 1.5-ncol(zt)*p)*log(det(sig1a)/det(sig0a))
pvp1 = 1 - pchisq(m1, ncol(zt)^2)
aci1 = log(det(sig1a)) + 2/nrow(zt)*p*ncol(zt)^2
bic1 = log(det(sig1a)) + log(nrow(zt))/nrow(zt)*p*ncol(zt)^2
hq1 = log(det(sig1a)) + 2*(log(log(nrow(zt)))) /nrow(zt)*p*ncol(zt)^2
p = 2
x = cbind(rep(1, nrow(zt3)), zt2, zt1)
y = zt3
b2 = solve(crossprod(x, x), crossprod(x,y))
a2 = y - x %*% b2
sig2a = crossprod(a2, a2) / (nrow(zt)-p)
m2 = -(nrow(x) - p - 1.5-ncol(zt)*p)*log(det(sig2a)/det(sig1a))
pvp2 = 1 - pchisq(m2, ncol(zt)^2)
aci2 = log(det(sig2a)) + 2/nrow(zt)*p*ncol(zt)^2
bic2 = log(det(sig2a)) + log(nrow(zt))/nrow(zt)*p*ncol(zt)^2
hq2 = log(det(sig2a)) + 2*(log(log(nrow(zt)))) /nrow(zt)*p*ncol(zt)^2
odr = matrix(c(aci0, bic0, hq0, pvp0, aci1, bic1, hq1, pvp1,aci2, bic2, hq2, pvp2), 4, 3)
odr## [,1] [,2] [,3]
## [1,] -2.015364 -3.652739 -3.803296
## [2,] -2.015364 -3.650246 -3.798311
## [3,] -2.015364 -3.651939 -3.801697
## [4,] 0.000000 0.000000 0.000000显然VAR(2)比VAR(1)/(0)要好，是否更高的lag依照 likelihood test 和 IC 更合适，可以用MTS中的VARorder函数检验:
VARorder(zt, 5)## selected order: aic = 2
## selected order: bic = 2
## selected order: hq = 2
## Summary table:
## p AIC BIC HQ M(p) p-value
## [1,] 0 -2.0156 -2.0156 -2.0156 0.0000 0.0000
## [2,] 1 -3.6528 -3.6503 -3.6520 49119.3575 0.0000
## [3,] 2 -3.8032 -3.7982 -3.8016 4529.0898 0.0000
## [4,] 3 -3.8029 -3.7955 -3.8005 10.0331 0.3478
## [5,] 4 -3.8029 -3.7929 -3.7997 15.5062 0.0779
## [6,] 5 -3.8028 -3.7903 -3.7988 14.6535 0.1009结果现实，依照IC，所以三项指标都提示lag p = 2, 但是依照likelihood test应该将lag 3 也考虑进去。哪项更合理没有定论，一般地，更多要考虑模型预测即 Out of Sample 的效果， 如上所言，Granger Causality。
Quelle
Tsay, Ruey S(2014): Multivariate Time Series Analysis With R and Financial Applications, Chapter 2.