多元时间序列分析基础二

来源：互联网发布：网络管理规章制度编辑：程序博客网时间：2024/06/01 09:28

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

上篇主要内容为:

根据均值，协方差矩阵模拟white noise 属性的多元时间序列，
根据给定的多元序列计算cross covariance/correlation matrix，以及
多元时间序列的Ljung Box Test的计算。

本篇首先模拟一个初始的多元时间序列，同样3个观测对象， 500个观测点。首先在假定:

White Noise属性的残差序列 at 的协方差给定

  matrix(c(2, .5, .3, .5, 1.5, .6, .3, .6, 3 ), 3, 3)

原始时间zt序列可以通过ARMA(1, 2)模型进行模拟和预测

⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ z 1, 3 z 1, 4 ⋮ z 1, 502 z 2, 3 z 2, 4 ⋮ z 2, 502 z 3, 3 z 3, 4 ⋮ z 3, 502 ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ = ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ϕ 1, 0 ϕ 1, 0 ⋮ ϕ 1, 0 ϕ 2, 0 ϕ 2, 0 ⋮ ϕ 2, 0 ϕ 3, 0 ϕ 3, 0 ⋮ ϕ 3, 0 ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ + ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ z 1, 2 z 1, 3 ⋮ z 1, 501 z 2, 2 z 2, 3 ⋮ z 2, 501 z 3, 2 z 3, 3 ⋮ z 3, 501 ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎡ ⎣ ⎢ ⎢ ϕ 1, 1 ϕ 2, 1 ϕ 3, 1 ϕ 1, 2 ϕ 2, 2 ϕ 3, 2 ϕ 1, 3 ϕ 2, 3 ϕ 3, 3 ⎤ ⎦ ⎥ ⎥' - ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ a 1, 2 a 1, 3 ⋮ a 1, 501 a 2, 2 a 2, 3 ⋮ a 2, 501 a 3, 2 a 3, 3 ⋮ a 3, 501 ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎡ ⎣ ⎢ θ 1, 11 θ 1, 21 θ 1, 31 θ 1, 12 θ 1, 22 θ 1, 32 θ 1, 13 θ 1, 23 θ 1, 33 ⎤ ⎦ ⎥' - ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ a 1, 1 a 1, 2 ⋮ a 1, 500 a 2, 1 a 2, 2 ⋮ a 2, 500 a 3, 1 a 3, 2 ⋮ a 3, 500 ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎡ ⎣ ⎢ θ 2, 11 θ 2, 21 θ 3, 31 θ 2, 12 θ 2, 22 θ 3, 32 θ 2, 13 θ 2, 23 θ 3, 33 ⎤ ⎦ ⎥' + ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ a 1, 3 a 1, 4 ⋮ a 1, 502 a 2, 3 a 2, 4 ⋮ a 2, 502 a 3, 3 a 3, 4 ⋮ a 3, 502 ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥

其中截距 ϕ0 ， AR部分系数和MA部分的系数 ϕp,θq都已经知道, 简单起见，假设所有截距为0，其他数据如下

library(MTS)print("phi_1")temp = matrix(c(.8 ,4, .5, -.3, -.6, .7, .9, -.6, -.8), 3, 3)tempprint("theta_1")t(temp)*1.3*-1print("theta_2")solve(temp)*0.6*-1.1

根据以上假设通过以下方式模拟一组原始的多元序列:

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) = 2phi0 = 0phi1 = matrix(c(.8 ,4, .5, -.3, -.6, .7, .9, -.6, -.8), 3, 3)theta1 = t(phi1)*1.3*-1 # arbitrary for simplicationtheta2 = solve(phi1)*0.6*-1.1at = rmvnorm(nT, rep(0, nrow(sig)), sigma = sig, set.seed(1000))head(at)at1 = 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,] =  z1i = 2while(i <= 500){  zt[i,] = z1 %*% t(phi1) + atm[i, ]  z1 = zt[i, ]  i = i + 1}head(zt)

使用MTS包中的VARMAsim函数如下, 便捷得到同样的结果，同时可以顺便看一下zt的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)MTSplot(mtsim$series)

图中看出，接近结尾处起伏及其大，说明这个序列是几何级增长的，非staionary，不适用初级的多元分析。用qm函数和ccm函数，会报错，Error: ρ0 不可逆。

备注

Rmarkdwon
书目：Tsay, Ruey S(2014): Multivariate Time Series Analysis With R and Financial Applications, Chapter 1.

阅读全文

0 0