多元时间序列分析基础一

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

多元时间序列分析的目的可以概括为分析不同时间序列的动态相关性， 从而得到单元序列分析更准确的预测。首先，我们生成k = 3个时间序列，y1,y2,y3, 每个序列长度T为500。从单元时间序列分析得知，分析的重点在残差，即矫正过的时间序列，比较宽松但是很实际的假设是时间序列是一个white noise， 所以这里模拟的序列从简单出发也是white noise。 假定我们已经知道三个序列的Covariance Matrix cm 为

temp = matrix(c(1.2, -0.3, 0.8, -.3, .85, 1.45, .8, 1.45, .68), 3, 3)colnames(temp) = c("z1","z2", "z3")rownames(temp) = colnames(temp)temp 

三个时间序列的平均数 mn=(0，0，0)′, 那么我们可以模拟得到这个三时间序列

require(mvtnorm) cm = matrix(c(1.2, -0.3, 0.1, -.3, 1.5, 0.2, .1, .2, .68),3,3)mn = c(0, 0, 0)set.seed(1000) # for  replication zr = rmvnorm(500, mean = mn, sigma = cm)colnames(zr) = c("z1","z2", "z3") dim(zr) head(zr, n = 5) var(zr) # check covariance matrixapply(zr, 2, mean) # check mean

我们可以从下图中看到，很明显，他们都是stationary的。 因为先验的设定了他们的covariance和mean，他们的走势没有特别相似或不相似。

library(MTS)MTSplot(zr)

很自然，我们会考虑到这些数列和自己的lag value 以及其他序列的lag vlaue 是否确有联系，这个问题可以用cross covariance/correlation matrices 中找到提示， 这两个matrices，相对于lag ℓ， 依序为如下定义：

Γℓ^=1T−1∑t=ℓ+1T(zt−μ^z)(zt−ℓ−μ^z)′

ρℓ^=D^−1Γℓ^D^−1,D^=diag(Γ0^0.5)=diag(σ^ii）

用上面的多元序列zr为基础，下面举例计算 ℓ 为0 和 2 的情况, 首先 ℓ 是0的计算过程

nT = nrow(zr)nTnV = ncol(zr)nVtem_zr = scale(zr, center = T, scale = F) # zr - mean(zr) by colhead(tem_zr)g0 = (t(tem_zr) %*% tem_zr) / (nT - 1)g0 # Gamma_0se = (diag(diag(g0)))^0.5seseinv = solve(se)rho0 = seinv %*% g0 %*% seinvrho0 

从 ρ0 可知，相关性很弱，接着着计算 ℓ 为2的例子

lag = 2zr1 = tem_zr[(lag + 1):nT,]zr2 = tem_zr[1:(nT - lag), ]g2 = (t(zr1) %*% zr2)/(nT - 1)g2rho2 = seinv %*% g2 %*% seinvrho2 

显然，几乎没有相关性存在，以上计算可以很容易的使用MTS包中的ccm函数得出

ccm(zr,2)$ccm

这个包中使用了判别时代替具体数值，具体地，根据某个cross correlation matrices ρℓ的中的元素ρℓ,ij生成新的matrix S：

Sℓ,ij=⎧⎩⎨⎪⎪+−.ρℓ,ij≥2/T−−√ρℓ,ij≤−2/T−−√ρℓ,ij<2/T−−√

类似单元时间序列，在多元情况下，我们也可以用Ljung-Box Test来检验 H0:ρ1=ρ2=...ρm=0, 同样检验的基准为 χ2分布：

Qk(m)=T2∑ℓ=1m1T−ℓtr(Γ^′ℓΓ^−10ΓℓΓ^−10)

对于序列zr，下面举例计算Q3(2)

g1 = t(tem_zr[(1+1): nT, ]) %*% tem_zr[1:(nT-1),]g1 = g1/(nT -1)g1trpd1 = t(g1) %*% solve(g0) %*% g1 %*% solve(g0)trs = sum(diag(trpd1)) / (nT - 1)Q31 = trs * (nT^2)Q31trpd2 = t(g2) %*% solve(g0) %*% g2 %*% solve(g0)trs2 = sum(diag(trpd2)) / (nT - 2) + trsQ32 = trs2 * (nT^2)Q32pv = 1 - pchisq(c(Q31, Q32), df = c(nV^2, 2*nV^2))pv

显然,在百分之五的confidence-level的状况下，可以拒绝说有在ρ1,ρ2中的元素都不相关， LJung Box Test 同样可以使用MTS 中的mq函数轻易实现：

mq(zr, lag = 2)  

####备注
Rmarkdown版本
书目：Tsay, Ruey S(2014): *Multivariate Time Series Analysis
With R and Financial Applications*, Chapter 1.