R/causality.R

In cheaton/vars2: VAR Modelling

causality <-
function(x, cause = NULL, vcov.=NULL, boot=FALSE, boot.runs=100){
  if(!(class(x)=="varest")){
    stop("\nPlease provide an object of class 'varest', generated by 'var()'.\n")
  }
  K <- x$K
  p <- x$p
  obs <- x$obs
  type <- x$type
  obj.name <- deparse(substitute(x))
  y <- x$y
  y.names <- colnames(x$y)
  if(is.null(cause)){
    cause <- y.names[1]
    warning("\nArgument 'cause' has not been specified;\nusing first variable in 'x$y' (", cause, ") as cause variable.\n")
  } else {
    if(!all(cause%in%y.names)) stop("Argument cause does not match variables names.\n")
  }
  y1.names <- subset(y.names, subset = y.names %in% cause)
  y2.names <- subset(y.names, subset = !(y.names %in% cause))
  Z <- x$datamat[, -c(1 : K)]
  xMlm<-toMlm(x)
  PI <- coef(xMlm)
  PI.vec <- as.vector(PI)

  ###Restriction matrix R for Granger causality
  #build matrix of same size as coef matrix indicating which to be restricted
  R2<-matrix(0, ncol=ncol(PI), nrow=nrow(PI)) 
  g<-which(gsub("\\.l[[:digit:]]", "", rownames(PI))%in%cause) #select cause regressors
  j<-which(colnames(PI)%in%cause) #select cause regressand
  R2[g,-j]<-1	#select coef to be tested
  w<-which(as.vector(R2)!=0)
  #build corresponding matrix as coef are not vectorized
  N <- length(w)
  R<-matrix(0, ncol=ncol(PI)*nrow(PI), nrow=N) #matrix of restrictions
  for(i in 1:N) R[i,w[i]]<-1

  ##
  ## Granger-causality
  ##
  sigma.pi  <- if (is.null(vcov.)) vcov(xMlm)
        else if (is.function(vcov.)) vcov.(xMlm) else vcov.
  df1 <- p * length(y1.names) * length(y2.names)
  df2 <- obs - dim(PI)[1]  # This line is df2 <- K * obs - length(PI) in the original vars package.
  STATISTIC <- t(R %*% PI.vec) %*% solve(R %*% sigma.pi %*% t(R)) %*% R %*% PI.vec / N

  ###bootstrap procedure
  if(boot){
    ###Restricted model: estimation under null of Granger non-causality
    co.names<-Bcoef(x)
    #needs to rebuild another restriction matrix for restrict(), as disposition of coef is different
    k<-which(gsub("\\.l[[:digit:]]", "", colnames(co.names))%in%cause) #select cause regressors
    l<-which(rownames(co.names)%in%cause) #select cause regressand
    R2inv<-matrix(1, ncol=nrow(PI), nrow=ncol(PI)) #exact inverse steps as R2
    R2inv[-l,k]<-0	#select coef to be tested  
    xres<-restrict(x, method = "man", resmat = R2inv)
    pred<-sapply(xres$varresult,predict)
    res<-residuals(xres)

    #bootstrap function for homo case: use more efficient low-level, as XX-1 already computed
    if(is.null(vcov.)){
      Zmlm<-model.matrix(xMlm)
      cross<-crossprod(Zmlm)
      inside<-solve(R %*% sigma.pi %*% t(R))
      boot.fun<-function(x=1){
	Ynew<-pred+res*rnorm(n=obs, mean=0, sd=x)
	PI.boot<-solve(cross, crossprod(Zmlm,Ynew)) #this could be made more efficent: compute only interest coefs,
	PI.boot.vec<-as.vector(PI.boot)
	t(R %*% PI.boot.vec) %*% inside %*% (R %*% PI.boot.vec) / N
      }
    } else {
    #bootstrap function for hetero case: obliged to run whole lm
    #two next lines as needed as x<-freeny.x; mylm<-lm(freeny.y~x); rm(x);update(mylm) #does not work
      X<-x$datamat 
      if(x$type%in%c("const", "both")) X<-X[, -grep("const", colnames(X))]
      boot.fun<-function(x=1){
	X[,1:K]<-pred+res*rnorm(n=obs, sd=x, mean=1) #workaround as calling it ynew and putting in update() fails
	xMlm.boot<-update(xMlm, .~.)
	sigma.pi.boot  <- if (is.function(vcov.)) vcov.(xMlm.boot) else {vcov.;
	  warning("vcov. should be function, not an object, when used with boot=TRUE")}
	PI.boot.vec <- as.vector(coef(xMlm.boot))
	t(R %*% PI.boot.vec) %*% solve(R %*% sigma.pi.boot %*% t(R)) %*% R %*% PI.boot.vec / N
      }
    }    
    res.rep<-replicate(boot.runs, boot.fun(x=1))
    pval<-mean(res.rep>as.numeric(STATISTIC))
  }
  names(STATISTIC) <- "F-Test"
  if(!boot){
    PARAMETER1 <- df1
    PARAMETER2 <- df2 
    names(PARAMETER1) <- "df1"
    names(PARAMETER2) <- "df2"  
    PVAL <- 1 - pf(STATISTIC, PARAMETER1, PARAMETER2)
    PARAM<-c(PARAMETER1, PARAMETER2)
  } else {
    PARAMETER1 <- boot.runs
    names(PARAMETER1) <- "boot.runs"
    PVAL <- pval
    PARAM<-PARAMETER1
  }
  METHOD <- paste("Granger causality H0:", paste(y1.names, collapse=" "), "do not Granger-cause", paste(y2.names, collapse=" "))
  result1 <- list(statistic = STATISTIC, parameter = PARAM, p.value = PVAL, method = METHOD, data.name = paste("VAR object", obj.name))
  class(result1) <- "htest"
  ##
  ## Instantaneous Causality
  ##
  sigma.u <- crossprod(resid(x)) / (obs - ncol(Z))
  colnames(sigma.u) <- y.names
  rownames(sigma.u) <- y.names
  select <- sigma.u[rownames(sigma.u) %in% y2.names, colnames(sigma.u) %in% y1.names ]
  sig.vech <- sigma.u[lower.tri(sigma.u, diag = TRUE)]
  index <- which(sig.vech %in% select)
  N <- length(index)
  Cmat <- matrix(0, nrow = N, ncol = length(sig.vech))
  for(i in 1 : N){
    Cmat[i, index[i]] <- 1
  }
  Dmat <- .duplicate(K)
  Dinv <- ginv(Dmat)
  lambda.w <- obs %*% t(sig.vech) %*% t(Cmat) %*% solve(2 * Cmat %*% Dinv %*% kronecker(sigma.u, sigma.u) %*% t(Dinv) %*% t(Cmat)) %*% Cmat %*% sig.vech
  STATISTIC <- lambda.w
  names(STATISTIC) <- "Chi-squared"
  PARAMETER <- N
  names(PARAMETER) <- "df"
  PVAL <- 1 - pchisq(STATISTIC, PARAMETER)
  METHOD <- paste("H0: No instantaneous causality between:", paste(y1.names, collapse=" "), "and", paste(y2.names, collapse=" "))
  result2 <- list(statistic = STATISTIC, parameter = PARAMETER, p.value = PVAL, method = METHOD, data.name = paste("VAR object", obj.name))
  class(result2) <- "htest"
  result2
  return(list(Granger = result1, Instant = result2))
}