compute the cross-correlation coefficient and check the p.value

fun_ccfmax.R

## Compute the cross-correlation Pearson product-moment correlation coefficient
## and check the p.value 
## This function is an alternative of the native ccf() function which does not 
## return p.value by default

fun.ccfmax <- function(x) {
 a <- x[1:(length(x) / 2)]
 b <- x[-(1:(length(x) / 2))]
 if (all(is.na(a)) | all(is.na(b))) {
   return(c(NA, NA, NA))
 }
 else
   p.value <- vector('numeric', 9)
   r <- vector('numeric', 9)

   ## from lag-0 to lag-8
   cort1 <- cor.test(a, b)
   p.value[1]<- as.numeric(cort1$p.value)
   r[1]<- as.numeric(cort1$estimate)

   cort2 <- cor.test(a[-(length(a))], b[-1])
   p.value[2]<- as.numeric(cort2$p.value)
   r[2]<- as.numeric(cort2$estimate)

   cort3 <- cor.test(a[-((length(a) - 1):(length(a)))], b[-(1:(1+1))])
   p.value[3]<- as.numeric(cort3$p.value)
   r[3]<- as.numeric(cort3$estimate)

   cort4 <- cor.test(a[-((length(a) - 2):(length(a)))], b[-(1:(1+2))])
   p.value[4]<- as.numeric(cort4$p.value)
   r[4]<- as.numeric(cort4$estimate)

   cort5 <- cor.test(a[-((length(a) - 3):(length(a)))], b[-(1:(1+3))])
   p.value[5]<- as.numeric(cort5$p.value)
   r[5]<- as.numeric(cort5$estimate)

   cort6 <- cor.test(a[-((length(a) - 4):(length(a)))], b[-(1:(1+4))])
   p.value[6]<- as.numeric(cort6$p.value)
   r[6]<- as.numeric(cort6$estimate)

   cort7 <- cor.test(a[-((length(a) - 5):(length(a)))], b[-(1:(1+5))])
   p.value[7]<- as.numeric(cort7$p.value)
   r[7]<- as.numeric(cort7$estimate)

   cort8 <- cor.test(a[-((length(a) - 6):(length(a)))], b[-(1:(1+6))])
   p.value[8]<- as.numeric(cort8$p.value)
   r[8]<- as.numeric(cort8$estimate)

   cort9 <- cor.test(a[-((length(a) - 7):(length(a)))], b[-(1:(1+7))])
   p.value[9]<- as.numeric(cort9$p.value)
   r[9]<- as.numeric(cort9$estimate)

   r[which(p.value > 0.05)] <- NA
   r.max <- max(r, na.rm = T)
   lag <- 1 - which.max(r)
   p <- p.value[which.max(r)]
   return(c(r.max, lag, p))
}