Automatically resolve perfect multi-collinearity?

[christophsax]

In case of perfect multi-collinearity, td() runs into a ‘singular matrix’ issue (example by @felicitaskemeny):

library(tempdisagg)
y <- ts(rnorm(n = 20), start = 1995)
x <- ts(rnorm(n = 80), start = 1995, freq = 4)
dummy <- ts(rep(0, 80), start = 1995, freq = 4)
disagg <- td(y ~  x + dummy)
#> Error in backsolve(R, c1 - C1 %*% v): singular matrix in 'backsolve'. First zero in diagonal [3]

The solution of course is to drop the dummy. But lm() is smarter and does that autmatically for you.

y2 <- ts(rnorm(n=80), start=1995, freq=4)
lm(y2 ~  x + dummy)
#> 
#> Call:
#> lm(formula = y2 ~ x + dummy)
#> 
#> Coefficients:
#> (Intercept)            x        dummy  
#>     0.01285      0.19844           NA

Should we have something similar in td()?

^{Created on 2020-09-09 by the reprex package (v0.3.0)}

[felicitaskemeny]

Hallo Christoph Das klingt super. Wenn ich es richtig verstehe, würde der Dummy erst dann für das Recursive berücksichtigt, wenn er zum ersten Mal den Wert 1 annimmt und davor ignoriert. Das fände ich gut. Liebe Grüsse Felicitas

[christophsax]

Here is an idea how constants could be droped. The main point is that a RHS variable needs to be constant within the LHS range.

drop_constants <- function(formula, envir = parent.frame()) {
  # which variables have 0 variance?
  ff <- formula

  # relevant time period
  tsp <- tsp(eval(ff[[2]], envir = envir))
  var_within_lhs <- function(x) {
    var(window(x, start = tsp[1], end = tsp[2]))
  }
  ff[[2]] <- NULL
  rhs_vars <- names(model.frame(ff))
  var <- vapply(lapply(setNames(rhs_vars, rhs_vars), get), var_within_lhs, 1)
  var0_vars <- names(var)[var == 0]
  if (length(var0_vars) == 0) return(formula)

  # update formula (e.g. . ~ -dummy)
  upd <- as.formula(paste("~ . ", paste(paste0("- ", var0_vars), collapse = " ")))
  update.formula(formula, upd)
}

It works for ts objects:

library(tempdisagg)
y <- ts(rnorm(n = 20), start = 1995)
x <- ts(rnorm(n = 80), start = 1995, freq = 4)
dummy <- ts(rep(0, 80), start = 1995, freq = 4)
td(drop_constants(y ~  x + dummy))
#> 
#> Call:
#> td(formula = drop_constants(y ~ x + dummy))
#> 
#> Coefficients:
#> (Intercept)            x  
#>    0.072943    -0.003609  
#> 
#> Use summary() for details. 
#> Use predict() to extract the final series.
#>       
#> Use ?td to see the help file.

A proper integration to td() is probably quite complicated with all the different time series classes tempdisagg.

The alternative would be to drop them further down, when we are only dealing with matrices.

^{Created on 2020-09-09 by the reprex package (v0.3.0)}