[Vars-commits] r67 - pkg/R

[noreply at r-forge.r-project.org]

Author: matthieu
Date: 2010-02-04 02:15:47 +0100 (Thu, 04 Feb 2010)
New Revision: 67

Modified:
   pkg/R/causality.R
Log:
Allow for bootstrap granger causality

Modified: pkg/R/causality.R
===================================================================
--- pkg/R/causality.R	2010-02-04 01:15:45 UTC (rev 66)
+++ pkg/R/causality.R	2010-02-04 01:15:47 UTC (rev 67)
@@ -1,5 +1,5 @@
-"causality" <-
-function(x, cause = NULL, vcov.=NULL){
+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")
   }
@@ -44,14 +44,51 @@
   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
+  if(boot){
+    pred<-predict(xMlm)
+    res<-residuals(xMlm)
+    Zmlm<-model.matrix(xMlm)
+    cross<-crossprod(Zmlm)
+    inside<-solve(R %*% sigma.pi %*% t(R))
+    #bootstrap function for homo case: use more efficient low-level, as XX-1 already computed
+    bootit.Homo<-function(x){
+      Ynew<-pred+res*rnorm(n=obs, sd=0)
+      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
+    }
+    #bootstrap function for hetero case: obliged to run whole lm
+    bootit.Hetero<-function(x){
+      Ynew<-pred+res*rnorm(n=obs)
+      #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))]
+      xMlm.boot<-update(xMlm, Ynew~.)
+      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
+    }
+    boot.fun<-if(is.null(vcov.)) bootit.Homo else bootit.Hetero
+    res.rep<-replicate(boot.runs, boot.fun(Xmlm))
+    pval<-mean(res.rep>as.numeric(STATISTIC))
+  }
   names(STATISTIC) <- "F-Test"
-  PARAMETER1 <- df1
-  PARAMETER2 <- df2 
-  names(PARAMETER1) <- "df1"
-  names(PARAMETER2) <- "df2"  
-  PVAL <- 1 - pf(STATISTIC, PARAMETER1, PARAMETER2)
+  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 = c(PARAMETER1, PARAMETER2), p.value = PVAL, method = METHOD, data.name = paste("VAR object", obj.name))
+  result1 <- list(statistic = STATISTIC, parameter = PARAM, p.value = PVAL, method = METHOD, data.name = paste("VAR object", obj.name))
   class(result1) <- "htest"
   ##
   ## Instantaneous Causality