gpuGranger: every pair in double by the Frisch-Waugh route; NA where unresolvable

DESCRIPTION

@@ -1,5 +1,5 @@
 Package: gputools
-Version: 2.1.0
+Version: 2.1.1
 Title: GPU Drop-Ins for Statistics on Wide Data
 Authors@R: c(person("Josh", "Buckner", email = "nullsatz@gmail.com",
                     role = c("aut", "cre")),

NEWS.md

@@ -1,3 +1,20 @@
+# gputools 2.1.1
+
+- `gpuGranger()` reported F = 0, p = 1, "no causality", for a pair it
+  could not resolve (collinear lags, a constant series, a cause identical
+  to the effect), computed in single precision by default on uncentred
+  data (F statistics of long random walks off by more than 10 percent),
+  and built one regression per pair with memory growing as the square of
+  the number of series. It now computes every pair in double by the
+  Frisch-Waugh route, with each effect's own past projected out once, and
+  the cause's test taken from a few matrix products, agreeing with `lm()`
+  to about 1e-8; a pair it cannot resolve in double is `NA` with a warning
+  (series integrated twice or more should be differenced first), a series
+  against itself or a constant series `NA` without one. `lag` is checked
+  (1 to 16). The `precision` argument is accepted but the tests always
+  compute in double. The method and its limits follow the Python package
+  widestats.
+
 # gputools 2.1.0
 
 - New `gpuCorPerm()`: permutation p-values for every correlation between

R/gpuGranger.R

@@ -1,21 +1,31 @@
 gpuGranger <- function(x, y=NULL, lag,
     precision = getOption("gputools.precision", "single"))
 {
+    # the argument is accepted for the common interface; the Granger tests
+    # always compute in double (see the help page)
     precision <- .precision(precision)
     x <- as.matrix(x)
     rows <- nrow(x)
     colsx <- ncol(x)
 
-    if(rows - lag <= 2 * lag + 1) {
+    if (ncol(x) < 1) stop("x needs at least one column")
+    if (length(lag) != 1 || !is.numeric(lag) || is.na(lag) ||
+        lag != round(lag) || lag < 1) {
+        stop("lag must be a positive whole number")
+    }
+    if (lag > 16) {
+        stop("lag must be at most 16")
+    }
+    lag <- as.integer(lag)
+    if (rows - lag <= 2 * lag + 1) {
         stop("time sequence too short for lag: use longer sequences or smaller lag")
     }
 
     if(!is.null(y)) {
         y <- as.matrix(y)
+        if (nrow(y) != rows) stop("x and y have different numbers of rows")
     }
 
-    lag <- as.integer(lag)
-
     if(is.null(y)) {
         colsy <- colsx
         cRetVal <- .Call(C_gpuGranger, x, as.integer(rows), as.integer(colsx),
@@ -26,6 +36,14 @@ gpuGranger <- function(x, y=NULL, lag,
                          as.integer(colsx), y, as.integer(colsy), lag,
                          precision)
     }
+    if (cRetVal$unresolved > 0) {
+        warning(sprintf(paste(
+            "%d pair(s) are too ill-conditioned to resolve in double",
+            "precision: a cause almost exactly a linear function of the",
+            "effect's own past, or nearly collinear lags (a series integrated",
+            "twice or more: difference it; or an exactly periodic one); their",
+            "tests are NA"), cRetVal$unresolved), call. = FALSE)
+    }
     # [i, j] tests whether x[, i] Granger-causes x[, j] (or y[, j]); the
     # dimnames say so when printed
     names <- function(m, prefix)
@@ -35,10 +53,9 @@ gpuGranger <- function(x, y=NULL, lag,
     dn <- list(cause = causes, effect = effects)
     fStats <- matrix(cRetVal$fStats, colsx, colsy, dimnames = dn)
     pValues <- matrix(cRetVal$pValues, colsx, colsy, dimnames = dn)
-    if (is.null(y)) {
-        # a series does not Granger-cause itself: no test, so no value
-        diag(fStats) <- NA
-        diag(pValues) <- NA
-    }
+    # no test (a series against itself, a constant series) or no answer
+    # (unresolvable) is NA, not NaN, as R reports a missing statistic
+    fStats[is.nan(fStats)] <- NA
+    pValues[is.nan(pValues)] <- NA
     return(list(fStatistics = fStats, pValues = pValues))
 }