Nothing
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))
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.