[Vars-commits] r107 - pkg/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sun Sep 12 14:44:15 CEST 2021
Author: bpfaff
Date: 2021-09-12 14:44:15 +0200 (Sun, 12 Sep 2021)
New Revision: 107
Modified:
pkg/R/causality.R
Log:
Fixed regular expression (github patch)
Modified: pkg/R/causality.R
===================================================================
--- pkg/R/causality.R 2021-09-12 12:41:30 UTC (rev 106)
+++ pkg/R/causality.R 2021-09-12 12:44:15 UTC (rev 107)
@@ -25,10 +25,22 @@
###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
+ 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)
@@ -38,8 +50,13 @@
##
## Granger-causality
##
- sigma.pi <- if (is.null(vcov.)) vcov(xMlm)
- else if (is.function(vcov.)) vcov.(xMlm) else vcov.
+ 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
@@ -49,10 +66,18 @@
###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
+ 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
+ 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)
@@ -59,6 +84,7 @@
#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))
@@ -68,20 +94,35 @@
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 {
- #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))]
+ xtmp <- x
+ # 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")}
+ 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))
}
@@ -88,9 +129,9 @@
names(STATISTIC) <- "F-Test"
if(!boot){
PARAMETER1 <- df1
- PARAMETER2 <- df2
+ PARAMETER2 <- df2
names(PARAMETER1) <- "df1"
- names(PARAMETER2) <- "df2"
+ names(PARAMETER2) <- "df2"
PVAL <- 1 - pf(STATISTIC, PARAMETER1, PARAMETER2)
PARAM<-c(PARAMETER1, PARAMETER2)
} else {
@@ -117,7 +158,7 @@
Cmat[i, index[i]] <- 1
}
Dmat <- .duplicate(K)
- Dinv <- ginv(Dmat)
+ 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"
@@ -130,4 +171,3 @@
result2
return(list(Granger = result1, Instant = result2))
}
-
More information about the Vars-commits
mailing list