矩阵中不同时间序列的互相关

Cross correlation of different time series in a matrix

在寻找我的问题的解决方案时,我发现了一个旧的 post (Cross correlation of different time series data values in R),它确切地询问了我需要什么,但不幸的是它没有得到任何答案所以我会再次询问希望得到一些指导.

我从大量具有相同大小的时间序列创建了一个大矩阵,每一列都是一个不同的时间序列(类似于以下内容但更大且更多的值不同于零):

      [,1]  [,2]  [,3]  [,4]  [,5]  [,6]  [,7]  [,8]  [,9] [,10] [,11] [,12] [,13] [,14] [,15] [,16] [,17] [,18] [,19]
[1,]    0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0    NA    NA    NA   0.0    NA   0.0   0.0   0.0   0.0
[2,]    0   6.0   0.0   9.0   0.0   0.0   0.0   0.0   0.0   0.0    NA     0    NA   0.0    NA   0.0   0.0   0.0   0.0
[3,]    0   0.0   0.0   5.0   0.0   0.0   0.0   0.0   0.0   0.0    NA     0    NA   0.0    NA   0.0   0.0   0.0   0.0
[4,]    0   0.0   0.0  10.0   0.0   0.0   0.0   0.0   0.0   0.0    NA     0    NA   0.0    NA   0.0   0.0   0.0   0.0
[5,]    0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0    NA     0    NA   0.0    NA   0.0   0.0   0.0   0.0
[6,]    0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0    NA     0    NA   0.0    NA   0.0   0.0   0.0   0.0
[7,]    0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0    NA     0    NA   0.0    NA   0.0   0.0   0.0   0.0
[8,]    0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0    NA     0    NA   0.0    NA   0.0   0.0   0.0   0.0
[9,]    0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0    NA     0    NA  10.0    NA   0.0   0.0   0.0   0.0
.
.
.

我想确定所有时间序列之间的相关性,我将它们放在一个矩阵中,因为我认为这可能是执行互相关程序的最佳方式,但我可能错了。

所以,我也知道函数 "ccf" 和 "diss()":

  1. ccf() #in base packages
  2. diss(meter_daywise,METHOD = "CORT",deltamethod = "DTW")#in TSclust 包

但和旧的 post 一样,我也有同样的问题:

  1. ccf不以全矩阵作为输入
  2. diss() 获取输入矩阵并生成一些矩阵,但在观察值时我发现它不是互相关矩阵,因为值不在 -1 和 1 之间。

所以问题是我们如何在 R 中计算和执行不同时间序列之间的互相关?

一种可能性是 运行 ccf 使用 combn 跨列的所有组合。以下代码针对link中的问题进行了测试:

myResults <- combn(seq_len(nrow(meter_daywise)), 2,
                   FUN=function(x) ccf(meter_daywise[x[1],], meter_daywise[x[2],]),
                   simplify=FALSE)

并生成这样的嵌套列表

str(myResults)
List of 10
 $ :List of 6
  ..$ acf   : num [1:17, 1, 1] 0.0241 0.0895 0.1463 0.0583 -0.0613 ...
  ..$ type  : chr "correlation"
  ..$ n.used: int 15
  ..$ lag   : num [1:17, 1, 1] -8 -7 -6 -5 -4 -3 -2 -1 0 1 ...
  ..$ series: chr "X"
  ..$ snames: chr "meter_daywise[x[1], ] & meter_daywise[x[2], ]"
  ..- attr(*, "class")= chr "acf"
 $ :List of 6
  ..$ acf   : num [1:17, 1, 1] -0.445 -0.493 -0.239 0.465 0.49 ...
  ..$ type  : chr "correlation"
  ..$ n.used: int 15
  ..$ lag   : num [1:17, 1, 1] -8 -7 -6 -5 -4 -3 -2 -1 0 1 ...
  ..$ series: chr "X"
  ..$ snames: chr "meter_daywise[x[1], ] & meter_daywise[x[2], ]"
  ..- attr(*, "class")= chr "acf"

...

列表中的每个外部元素都是 ccf 两对的输出。对于您的应用程序,由于时间序列存储在列中,因此您将其切换为

myResults <- combn(seq_len(ncol(myMat)), 2,
                   FUN=function(x) ccf(myMat[, x[1]], myMat[, x[2]]), simplify=FALSE)

其中 myMat 是矩阵的名称。您可以通过更简单地调用 combn 来查看这些对,例如

myPairs <- combn(seq_len(ncol(myMat)), 2)

ccf returns 每个偏移量的成对相关性(即滞后),但我认为你想要的是最大(绝对值(相关性)来自那个。因为你有 NA,你需要设置na.action 参数。

mat <- matrix(rnorm(100000), ncol=100)
mat[sample(1:length(mat), 100)] <- NA 

res <- sapply(1:ncol(mat), function(x) {
  sapply(1:ncol(mat), function(z){
    resTmp <- ccf(x = mat[, x], y = mat[, z], plot=F, na.action = na.pass)
    resTmp$acf[which.max(abs(resTmp$acf))]
  })
})

来自 ccf 帮助:

By default, no missing values are allowed. If the na.action function passes through missing values (as na.pass does), the covariances are computed from the complete cases. This means that the estimate computed may well not be a valid autocorrelation sequence, and may contain missing values.