R/testCCF.R

#' Calulate the CCF second test value
#'
#' Takes in the result of CCF
#' @param ccfResult The result of CCF function
#' @return A number of CCF sencond test value
#' @details CCF second value = (max(cffResult(ro>0))-max(cffResult(ro<0)))/(max(cffResult(ro>0))+max(cffResult(ro<0)))
#' @export
calCCF<-function(ccfResult){
  maxPo<-max(ccfResult$CCF[which(ccfResult$lag>0)])
  maxNe<-max(ccfResult$CCF[which(ccfResult$lag<0)])
  ccfValue<-(maxPo-maxNe)/(maxPo+maxNe)
  return(ccfValue)
}

#' Statistic Test for CCF
#'
#' Takes in five parameters, and test whether there is causality between two variables.
#'
#' @param x A vector represents a time sequence
#' @param y A vector represents another sequence
#' @param lag.max Max time lag
#' @param alpha Significant level
#' @param sampleTimes The times of permutation sample for second order test
#'
#'@return When pho < 0, it will return a sentence for no causality. When first order
#'  test is unsatisfied, it will return a sentence. When the second order test is
#'  unsatisfied, it will return a sentence and list of Q and sampleQ vector. When all
#'  the test is passed, it will return Q
#'@details First order test is to test whether the ccfvalue is different from zero.
#'  Second Test is to test whether the ccfvalue satisfies the 3-sigma test
#'@examples data("BMData")
#' i=8
#' u<-BMData[[i]]$u
#' y<-BMData[[i]]$y
#' resultw<-testCCF(u,y,lag.max=97,alpha = 0.05,sampleTimes=10000)
#' resultq<-testCCF(u,y,lag.max=97,alpha = 0.05,sampleTimes=1000)
#' resultb<-testCCF(u,y,lag.max=97,alpha = 0.05,sampleTimes=100)
#' datab<-cbind(ccf=resultb[[2]],sampleTime=rep(100,100))
#' dataq<-cbind(ccf=resultq[[2]],sampleTime=rep(1000,1000))
#' dataw<-cbind(ccf=resultw[[2]],sampleTime=rep(10000,10000))
#' data<-rbind(datab,dataq,dataw)
#' data<-as.data.frame(data)
#' data$sampleTime<-as.factor(data$sampleTime)
#' library(plyr)
#' cdata <- ddply(data, "sampleTime", summarise, ccf.mean=mean(ccf),
#'                ccf.3sigma=mean(ccf)+3*sd(ccf))
#' library(ggplot2)
#' ggplot(data) + geom_density(aes(x = ccf,
#'                                 colour = sampleTime),
#'                             size=1)+labs(x=expression(paste("C"[pi],""^"CCF")))+
#'   labs(y=expression(paste("p","(C"[pi],""^"CCF",")")))+
#'   labs(title = "2nd order system")+
#'   geom_vline(data=cdata, aes(xintercept=ccf.3sigma,  colour=sampleTime),
#'              linetype="dashed", size=1)+
#'   geom_vline(aes(xintercept=result[[1]]),
#'              linetype="dashed", size=1)
#'@export
testCCF<-function(x,y,lag.max,alpha,sampleTimes,...){
  ccfResult<-CCF(x,y,lag.max)
  location<-which.max(ccfResult$CCF)
  N<-which.max(ccfResult$lag)+1
  pho<-ccfResult$lag[location]
  if(pho>0){
    ro<-ccfResult$CCF[location]
    tscore<-ro*sqrt((N-2)/1-ro^2)
    threhold<-qt(1-alpha/2,N-2)
    if(tscore>threhold){
      ccfValue<-calCCF(ccfResult)
      ccfSample<-unlist(lapply(c(1:sampleTimes),function(i,x,y,lag.max){
        xSample<-sample(x)
        ccfSResult<-CCF(xSample,y,lag.max)
        ccfSValue<-calCCF(ccfSResult)
        return(ccfSValue)
      },x,y,lag.max))
      meanSample<-mean(ccfSample)
      sdSample<-sd(ccfSample)
      threholdJudge<-meanSample+3*sdSample
      if(ccfValue>threholdJudge){
        type<-0
        judge<-"There is causality"
        result<-list("type"=type,"judge"=judge,"ccfValue"=ccfValue,"ccfSample"=ccfSample)
        class(result)<-"testCCF"
        attr(result,"heading")<-c("Cross Correlation Function test:")
        return(result)
      }
      else{
        type<-1
        judge<-"There is no causality becuase second test is unsatisfied"
        result<-list("type"=type,"judge"=judge,"ccfValue"=ccfValue,"ccfSample"=ccfSample)
        class(result)<-"testCCF"
        attr(result,"heading")<-c("Cross Correlation Function test:")
        return(result)
      }
    }
    else{
      type<-2
      judge<-"There is no causality because student test is unsatisfied"
      result<-list("type"=type,"judge"=judge)
      class(result)<-"testCCF"
      attr(result,"heading")<-c("Cross Correlation Function test:")
      return(result)
    }
  }
  else{
    type<-3
    judge<-"There is no causality because pho is negative"
    result<-list("type"=type,"judge"=judge)
    class(result)<-"testCCF"
    attr(result,"heading")<-c("Cross Correlation Function test:")
    return(result)
  }
}

print.testCCF<-function(x,...){
  print(attr(x,"heading"))
  print(x$judge)
}
chengfeifan/Causality documentation built on May 13, 2019, 3:40 p.m.