多元时间序列分析基础二
Chang
2017年11月17日
上篇主要内容为:
- 根据均值，协方差矩阵模拟white noise 属性的多元时间序列，
- 根据给定的多元序列计算cross covariance/correlation matrix， 以及
- 多元时间序列的Ljung Box Test的计算。
本篇首先模拟一个初始的多元时间序列， 同样3个观测对象， 500个观测点。首先在假定:
- White Noise属性的残差序列 \(a_t\) 的协方差给定
## [,1] [,2] [,3]
## [1,] 2.0 0.5 0.3
## [2,] 0.5 1.5 0.6
## [3,] 0.3 0.6 3.0- 原始时间\(z_t\)序列可以通过ARMA(1, 2)模型进行模拟和预测
\[\begin{bmatrix}z_{1,3}&z_{2,3}&z_{3,3}\\z_{1,4}&z_{2,4}&z_{3,4}\\\vdots&\vdots&\vdots\\z_{1,502}&z_{2,502}&z_{3,502}\end{bmatrix} =\begin{bmatrix}\phi_{1,0}&\phi_{2,0}&\phi_{3,0}\\\phi_{1,0}&\phi_{2,0}&\phi_{3,0}\\\vdots&\vdots&\vdots\\\phi_{1,0}&\phi_{2,0}&\phi_{3,0}\end{bmatrix} + \begin{bmatrix}z_{1,2}&z_{2,2}&z_{3,2}\\z_{1,3}&z_{2,3}&z_{3,3}\\\vdots&\vdots&\vdots\\z_{1,501}&z_{2,501}&z_{3,501}\end{bmatrix} \begin{bmatrix}\phi_{1,1}&\phi_{1,2}&\phi_{1,3}\\\phi_{2,1}&\phi_{2,2}&\phi_{2,3}\\\phi_{3,1}&\phi_{3,2}&\phi_{3,3}\end{bmatrix}' \\ - \begin{bmatrix}a_{1,2}&a_{2,2}&a_{3,2}\\a_{1,3}&a_{2,3}&a_{3,3}\\\vdots&\vdots&\vdots\\a_{1,501}&a_{2,501}&a_{3,501}\end{bmatrix} \begin{bmatrix}\theta_{1,11}&\theta_{1,12}&\theta_{1,13}\\\theta_{1,21}&\theta_{1,22}&\theta_{1,23}\\\theta_{1,31}&\theta_{1,32}&\theta_{1,33}\end{bmatrix}' - \begin{bmatrix}a_{1,1}&a_{2,1}&a_{3,1}\\a_{1,2}&a_{2,2}&a_{3,2}\\\vdots&\vdots&\vdots\\a_{1,500}&a_{2,500}&a_{3,500}\end{bmatrix} \begin{bmatrix}\theta_{2,11}&\theta_{2,12}&\theta_{2,13}\\\theta_{2,21}&\theta_{2,22}&\theta_{2,23}\\\theta_{3,31}&\theta_{3,32}&\theta_{3,33}\end{bmatrix}'\\ + \begin{bmatrix}a_{1,3}&a_{2,3}&a_{3,3}\\a_{1,4}&a_{2,4}&a_{3,4}\\\vdots&\vdots&\vdots\\a_{1,502}&a_{2,502}&a_{3,502}\end{bmatrix}\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\quad\ \]
- 其中截距 \(\phi_0\) ， AR部分系数和MA部分的系数 \(\phi_p, \theta_q\)都已经知道, 简单起见，假设所有截距为0，其他数据如下
## [1] "phi_1"## [,1] [,2] [,3]
## [1,] 0.8 -0.3 0.9
## [2,] 4.0 -0.6 -0.6
## [3,] 0.5 0.7 -0.8## [1] "theta_1"## [,1] [,2] [,3]
## [1,] -1.04 -5.20 -0.65
## [2,] 0.39 0.78 -0.91
## [3,] -1.17 0.78 1.04## [1] "theta_2"## [,1] [,2] [,3]
## [1,] -0.225 -0.0975 -0.18
## [2,] -0.725 0.2725 -1.02
## [3,] -0.775 0.1775 -0.18根据以上假设通过以下方式模拟一组原始的多元序列:
library(mvtnorm)
library(MTS)
sig = matrix(c(2, .59, 1.3, .59, 1.5, .69, 1.3, .69, 3 ), 3, 3)
nV = nrow(sig)
nT = 500 + 2 # max(p = 1, q = 2) = 2
phi0 = 0
phi1 = matrix(c(.8 ,4, .5, -.3, -.6, .7, .9, -.6, -.8), 3, 3)
theta1 = t(phi1)*1.3*-1 # arbitrary for simplication
theta2 = solve(phi1)*0.6*-1.1
at = rmvnorm(nT, rep(0, nrow(sig)), sigma = sig, set.seed(1000))
head(at)## [,1] [,2] [,3]
## [1,] -0.8175811 -1.5143637 -0.3743104
## [2,] 0.5368340 -0.8908759 -0.5419779
## [3,] -0.5008043 0.7580067 -0.0772550
## [4,] -2.2617668 -1.5590128 -1.7082311
## [5,] -0.4213646 -0.4035595 -2.2022471
## [6,] 0.2684216 0.2235346 0.1457351at1 = at[2:(nrow(at)-1),]
at2 = at[1:(nrow(at)-2),]
at3 = at[3:nrow(at),]
atm = at3 - at1 %*% t(theta1) - at2 %*% t(theta2)
zt = matrix(0, 500, 3)
z1 = atm[1, ]
zt[1,] = z1
i = 2
while(i <= 500){
zt[i,] = z1 %*% t(phi1) + atm[i, ]
z1 = zt[i, ]
i = i + 1
}
head(zt)## [,1] [,2] [,3]
## [1,] -5.326320 0.1884459 1.377179
## [2,] -2.032941 -24.1907492 -5.961443
## [3,] -11.777927 9.4511796 -15.573340
## [4,] -30.941852 -47.6971575 13.659760
## [5,] 1.281234 -105.5183267 -56.810364
## [6,] -16.705131 107.2690428 -32.978388使用MTS包中的VARMAsim函数如下, 便捷得到同样的结果，同时可以顺便看一下\(z_t\)的plots以及cross corelation:
library(MTS)
set.seed(1000)
mtsim = VARMAsim(500, skip = 2, arlags = c(1), malags = c(1,2),
cnst=NULL, phi = phi1, theta = cbind(theta1, theta2),
sigma = sig)
head(mtsim$series)## [,1] [,2] [,3]
## [1,] -5.326320 0.1884459 1.377179
## [2,] -2.032941 -24.1907492 -5.961443
## [3,] -11.777927 9.4511796 -15.573340
## [4,] -30.941852 -47.6971575 13.659760
## [5,] 1.281234 -105.5183267 -56.810364
## [6,] -16.705131 107.2690428 -32.978388MTSplot(mtsim$series)
图中看出，接近结尾处起伏及其大，说明这个序列是几何级增长的，非staionary， 不适用初级的多元分析。用qm函数和ccm函数，会报错，Error: \(\rho_0\) 不可逆。