R/causality.R

In vars: VAR Modelling

causality <-
function(x, cause = NULL, vcov.=NULL, boot=FALSE, boot.runs=100){
  if(!is(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\\d+", "", 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
  #If the model already has restriction, overlay with the new ones
  if (!is.null(x$restrictions)) {
      xr <- t(x$restrictions)
      xr <- abs(xr - 1)
      # match positions of variables
      rownames(xr)[rownames(xr) == "const"] <- "(Intercept)"
      xr <- xr[rownames(PI), colnames(PI)]
      # overlay
      xr <- xr + R2
      xr[xr == 2] <- 1
      R2 <- xr
  }
  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
  ##
  if (is.null(vcov.)) {
      sigma.pi  <- vcov(xMlm)
  } else if (is.function(vcov.)) {
      sigma.pi  <- vcov.(xMlm)
  } else {
      sigma.pi  <- vcov.
  }
  df1 <- p * length(y1.names) * length(y2.names)
  df2 <- K * obs - length(PI)#K^2 * p - detcoeff
  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\\d+", "", 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
    #If the model already has restriction, overlay with the new ones
    if (!is.null(x$restrictions)) {
        xr <- x$restrictions
        # match positions of variables
        xr <- xr[rownames(co.names), colnames(co.names)]
        # overlay
        R2inv <- xr * R2inv
    }
    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.)){
        if (is.null(x$restrictions)) { #haven't figured out how to adjust these lines to account for existing restrictions
      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 { #if restrictions already exist; reestimate models (slower), use vcov by default
            xtmp <- x
            boot.fun <- function(x = 1) {
                xtmp$datamat[,1:K] <- pred + res * rnorm(n = obs, mean = 0, sd = x)
                xMlm.boot <- toMlm(xtmp)
                sigma.pi.boot  <- vcov(xMlm.boot)
                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
            }
        }
    } else {
    #two next lines as needed as x<-freeny.x; mylm<-lm(freeny.y~x); rm(x);update(mylm) #does not work
        xtmp <- x
      # X<-x$datamat
      # if(x$type%in%c("const", "both")) X<-X[, -grep("const", colnames(X))]
      boot.fun<-function(x=1){
          xtmp$datamat[,1:K]<-pred+res*rnorm(n=obs, sd=x, mean=0) #workaround as calling it ynew and putting in update() fails
	# xMlm.boot<-update(xMlm, .~.) #replace with the row below to account for possible restrictions
          xMlm.boot <- toMlm(xtmp)
          if (is.function(vcov.)) {
              sigma.pi.boot <- vcov.(xMlm.boot)
          } else {
              sigma.pi.boot <- 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 <- MASS::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))
}