windows 计算的两个时间序列之间的相关系数在时间上向前移动了 n 个时间单位

correlation coefficients between two time series calculated over windows moved forward in time by n time unit

是否有生成图的包或简单代码 (1) 在 windows 上计算的两个时间序列之间的相关系数在时间上向前移动了 n 个时间单位 (2) 以及它们各自为每次移动计算的 p 值 ?

library(zoo)

x = ts(rnorm(1:121), start = 1900, end = 2021)
y = ts(rnorm(1:121), start = 1900, end = 2021)
data = data.frame(x, y)

# 40-year moving window lagged forward by 15 years per example

rollapply(data, width=40, by = 15, 
          function(x) cor(x[,1],x[,2], method =  "pearson"),
          by.column=FALSE)

[1]  0.92514750  0.5545223 -0.207100231 -0.119647462 -0.125114237  0.041334073

** Hmisc::rcorr 会更好,它也计算 p 值,但我没有设法将它集成到 rollapply.

这里的结果中,第一个系数(0.9251...)对1900:1940有效,第二个对1915:1955有效等

所以问题是:有没有一种快速的方法可以将此结果整合到具有时间、r 和 p 值的阶梯图中?

输出如下:

时间 r P
1900 0.92 0.000001
1901 0.92 0.000001
... ... ...
1915 0.55 0.00045
1916 0.55 0.00045

几点:

  • 从1900年到2021年有2021-1900+1 = 122年,不是121年
  • 40/15 参数不能均匀地使用 122 点所以从 1907 开始

rcorr returns 一个包含 3 个组件的列表,我们需要每个组件的 1,2 个元素。我们可以使用 na.locf 填充 rollapplyr 中的缺失值。输入和输出都是mts/ts系列

library(zoo)
library(Hmisc)

set.seed(123)
tt <- ts(cbind(x = rnorm(115), y = rnorm(115)), start = 1907)

na.locf(rollapplyr(tt, width=40, by = 15, 
          function(x) sapply(rcorr(x), `[`, 1, 2),
          by.column = FALSE, fill = NA), fromLast = TRUE)

上面的 returns 一个系列,其行数与输入 tt 相同,但基于以下年份范围的计算 rcorr:

rollapplyr(1907:2021, 40, by = 15, range)
##      [,1] [,2]
## [1,] 1907 1946
## [2,] 1922 1961
## [3,] 1937 1976
## [4,] 1952 1991
## [5,] 1967 2006
## [6,] 1982 2021