R中的部分互相关答案

【问题标题】：Partial Cross-correlation in RR中的部分互相关
【发布时间】：2016-08-24 20:09:16
【问题描述】：

我认为标题是不言自明的。我想计算为其他滞后值控制的两个时间序列之间的互相关。我找不到任何现有的 R 代码来执行此操作，而且我对自己的统计知识（或 R）没有足够的信心尝试自己编写一些东西。它类似于偏自相关函数，只是用于互相关而不是自相关。

如果它有帮助的话，我更大的目标是寻找物理系统的不同测量值之间的滞后相关性（首先，来自耀变体伽马射线测量的通量和光子指数），目标是建立一个一般线性模型来尝试预测燃烧事件。

【问题讨论】：

我在这里有点不够深入，但您可能可以查看函数 acf 的文档。尤其是使用 Union acf(ts.union(mdeaths, fdeaths)) 的例子
编辑：抱歉，在这种情况下，工会不会有任何好处。
除非您有特定的、可实施的方法来计算您想要的值，否则这可能不是问您问题的最佳地点。如果您在分析数据方面需要帮助，最好在Cross Validated 询问统计方法。
我会把它贴在那里，谢谢。

【解决方案1】：

查看我的answer 对我自己的问题（与您发布的问题相同）。

您可以使用 R 中的 pacf 函数，将其扩展到具有 2 个或更多时间序列的矩阵。我检查了多变量 acf 和 ccf 函数之间的结果，它们产生了相同的结果，因此对于多变量 pacf 和不存在的 pccf 可以得出相同的结论。

【讨论】：

【解决方案2】：

我相信这项工作，

  pccf <- function(x,y,nlags=7,partial=TRUE){
  
  # x (numeric): variable that leads y
  # y (numeric): variable of interest
  # nlags (integer): number of lags (uncluding zero)
  # partial (boolean): partial or absolute correlation
    
  # trim y
  y <- y[-(1:(nlags-1))]
  
  # lagged matrix of x
  x_lagged <- embed(x,nlags)
  
  # process for each lag
  rho <- lag <- NULL
  for(i in 1:(nlags)){
    
    if(partial){
      # residuals of x at lag of interest regressed on all other lags of x
      ex <- lm(x_lagged[,i] ~ x_lagged[,-i])$residuals
      
      # residuals of y regressed on all lags of x but the one of interest
      ey <- lm(y ~ x_lagged[,-i])$residuals
    }else{
      ex <- x_lagged[,i]
      ey <- y
    }
    
    # calculate correlation
    rho[i] = cor(ex,ey, use="pairwise.complete.obs")
    lag[i] = i-1
  }
  
  return(
    tibble(lag=lag, rho=rho) %>%
           arrange(lag)
    )
  
}

# test
n <- 200 # count
nlag <- 6 # number of lags
x <- as.numeric(arima.sim(n=n,list(ar=c(phi=0.9)),sd=1)) # simulate times series x
y <- lag(x,nlag) + rnorm(n,0,0.5) # simulate y to lag x
y <- y[(nlag+1):n] # remove NAs from lag
x <- x[(nlag+1):n] # align with y

pccf(x,y,nlags=10,partial=FALSE) %>%
 mutate(type='Cross correlation') %>%
 bind_rows(
  pccf(x,y,nlags=10,partial=TRUE) %>%
  mutate(type='Partial cross correlation') 
) %>%
ggplot() +
geom_col(aes(-lag,rho),width=0.1) +
facet_wrap(~type,scales='free_y', ncol=1) +
scale_x_continuous(breaks=-10:0) +
theme_bw(base_size=20)

【讨论】：