# Rolling regressions in R

I want to run rolling regressions over each group and store the coefficient. The following works, but it's slow, since I have too many series and I want to run too many regressions for each group.

Using foreach(), speeds things up (and also getting the coefficient with (X'X)^{-1}X'Y but is there a way to vectorize this operation? Also open for any and all improvements.

library(data.table)

run.rolling.regressions <- function(x) {
DT <- data.table(   Y = rnorm(10000),
X = rnorm(10000),
key.group = rep(LETTERS[1:10], each = 1000))

window.length <- 12
names.of.groups <- unique(DT$key.group) number.of.groups <- length(names.of.groups) X.coefficients <- list() for(j in 1:length(names.of.groups)) { regressed.DT <- DT[key.group == names.of.groups[j]] nrows.of.group <- nrow(regressed.DT) print(paste0(j, ', key.group: ', names.of.groups[j])) for (i in window.length:nrows.of.group) { if(i == window.length) { X.coefficients[[names.of.groups[j]]] <- c(rep(NA, nrows.of.group)) } X.coefficients[[names.of.groups[j]]][i] <- lm(Y ~ 1 + X, data = regressed.DT[(i - 11):i])$coefficients['X']

}
}
return(X.coefficients)
}
system.time(X.coef <- run.rolling.regressions())

• Welcome to Code Review. Are you looking for improvements on all aspects of your code, or for an answer to your specific programming question? – Phrancis Apr 12 '16 at 21:06
• Everything you can think of is welcome. I 've read "R inferno" so this is why I preallocate with X.coefficients[[names.of.groups[j]]] <- c(rep(NA, nrows.of.group)) }. I do not consider myself a programmer, though, since I do not do anything creative other than that. :) – Konstantinos Apr 12 '16 at 21:10
• Okay I have made a small edit to your post to clarify this, sometimes questions will get closed here if they are too specific about one thing only. I hope you get some good answers! – Phrancis Apr 12 '16 at 21:26

For each group in your data table, your code computes the coefficient b1 from a linear regression y = b0 + b1*x + epsilon, and you want to run this regression and obtain b1 for observations 1-12, 2-13, 3-14, ..., 989-1000. Right now you are separately calling lm for each data subset, which is a non-vectorized approach.

Vectorization of prediction models across datasets is in general not straightforward, but for the special case you have here (simple linear regression) is it possible because there is a simple closed-form expression for b1, the coefficient of interest. In particular, for given vectors x and y we have b1 = (mean(x*y) - mean(x)*mean(y)) / (mean(x^2) - mean(x)^2). The rolling coefficient value can therefore be computed using the rolling means of x*y, x, y, and x^2 with the appropriate window width.

The end result is a fully vectorized version of the code (I use the RcppRoll package to obtain rolling means):

library(RcppRoll)
rolling2 <- function(DT, window.length) {
setNames(lapply(unique(DT$key.group), function(g) { regressed.DT <- DT[key.group == g] xyBar = roll_mean(regressed.DT$X*regressed.DT$Y, window.length) xBar = roll_mean(regressed.DT$X, window.length)
yBar = roll_mean(regressed.DT$Y, window.length) x2Bar = roll_mean(regressed.DT$X^2, window.length)
c(rep(NA, window.length-1), (xyBar - xBar*yBar) / (x2Bar - xBar^2))
}), unique(DT$key.group)) }  We can confirm that this yields identical results to the code from the original post about 3 orders of magnitude more quickly: set.seed(144) DT <- data.table( Y = rnorm(10000), X = rnorm(10000), key.group = rep(LETTERS[1:10], each = 1000)) system.time(X.coef <- run.rolling.regressions(DT, 12)) # user system elapsed # 13.321 0.098 13.504 system.time(X.coef2 <- rolling2(DT, 12)) # user system elapsed # 0.010 0.000 0.011 all.equal(X.coef, X.coef2) # [1] TRUE  Note that I slightly modified the provided run.rolling.regressions function to take DT and window.length as input and to not print progress updates; I think it makes sense to separate the generation of the dataset from the function that computes the rolling means, and down the road it might be useful to have the window length as an adjustable argument instead of a fixed value. • Welcome to Code Review! Great job on your first answer. – SirPython Apr 13 '16 at 0:21 • Today I am sure I was one of the happiest persons on earth. Thank you very much! :) – Konstantinos Apr 13 '16 at 22:45 • The computations you do in the denominator can lead to catastrophic cancellation as I show in very similar example here. An easy solution is to use the rollRegres package I have written instead or another more stable method. – Benjamin Christoffersen Jan 24 at 23:38 Here is another solution which uses the rollRegres package # simulate data library(data.table) set.seed(144) DT <- data.table( Y = rnorm(10000), X = rnorm(10000), key.group = rep(LETTERS[1:10], each = 1000)) # assign suggested function library(rollRegres) f <- function(SD, width = 12L) roll_regres.fit(x = cbind(1, SD$X), y = SD$Y, width = 12L)$coefs[, 2]
o <- DT[, .(beta = f(.SD)), by = key.group]

# assign function to compare with
library(RcppRoll)
rolling2 <- function(DT, window.length) {
setNames(lapply(unique(DT$key.group), function(g) { regressed.DT <- DT[key.group == g] xyBar = roll_mean(regressed.DT$X*regressed.DT$Y, window.length) xBar = roll_mean(regressed.DT$X, window.length)
yBar = roll_mean(regressed.DT$Y, window.length) x2Bar = roll_mean(regressed.DT$X^2, window.length)
c(rep(NA, window.length-1), (xyBar - xBar*yBar) / (x2Bar - xBar^2))
}), unique(DT$key.group)) } X.coef2 <- rolling2(DT, 12) # gets the same all.equal(unlist(X.coef2), o$beta, check.attributes = FALSE)
#R TRUE

microbenchmark::microbenchmark(
X.coef2 = rolling2(DT, 12),
roll_regres.fit = DT[, .(beta = f(.SD)), by = key.group])
#R Unit: milliseconds
#R             expr   min    lq   mean median     uq   max neval
#R          X.coef2 9.453 9.726 20.729 10.529 13.705 440.0   100
#R  roll_regres.fit 4.235 4.312  7.173  4.374  4.517 136.9   100


I make a comparision with the answer from josliber. The implementation is faster as shown above and more numerically stable.

# Update

The package uses the LINPACK routines dchdd and dchud to update the Cholesky decomposition in the triangular matrix of the QR decomposition of the design matrix. See the LINPACK user guide at chapter 8 and 9. It is also be numerically stable unlike the aforementioned answer which can be subject have catastrophic cancellation as shown below

set.seed(60679655)
X <- rnorm(100) + 1e6
Y <- rnorm(100, sd = .01) + X + 1e6 # large intercept!

library(data.table)
DT <- data.table(X, Y, key.group = "a")
o1 <- DT[, .(beta = f(.SD))]
o2 <- rolling2(DT, 12)
library(zoo)
lm_res <- rollapply(data.frame(1, X, Y), 12L, function(x)
lm.fit(x[, c("X1", "X")], x[, "Y"])coefficients[2L], by.column = FALSE, fill = NA_real_, align = "right") # compare the error c_func <- function(x, y) mean(abs(x[[1L]] - y), na.rm = TRUE) c_func(o1, lm_res) #R [1] 2.067e-09 c_func(o2, lm_res) #R [1] 0.0003046 # redo above but with larger intercept set.seed(60679655) X <- rnorm(100) + 1e7 Y <- rnorm(100, sd = .01) + X + 1e7 # slope is 1 DT <- data.table(X, Y, key.group = "a") o1 <- DT[, .(beta = f(.SD))] o2 <- rolling2(DT, 12) lm_res <- rollapply(data.frame(1, X, Y), 12L, function(x) lm.fit(x[, c("X1", "X")], x[, "Y"])coefficients[2L],
by.column = FALSE, fill = NA_real_, align = "right")
c_func(o1, lm_res)
#R [1] 1.879e-08
c_func(o2, lm_res)
#R [1] 0.01634


The above is run with this version of the package which is not yet on CRAN. It does not matter at all for moderate means relative to the noise

set.seed(60679655)
X <- rnorm(100) + 1e3
Y <- rnorm(100, sd = 2) + X + 1e3 # slope is 1
DT <- data.table(X, Y, key.group = "a")
o1 <- DT[, .(beta = f(.SD))]
o2 <- rolling2(DT, 12)
lm_res <- rollapply(data.frame(1, X, Y), 12L, function(x)
lm.fit(x[, c("X1", "X")], x[, "Y"])\$coefficients[2L],
by.column = FALSE, fill = NA_real_, align = "right")
c_func(o1, lm_res)
#R [1] 1.982e-12
c_func(o2, lm_res)
#R [1] 3.047e-10