Hi,
A not unusual task is performing a multiple regression in a rolling window
on a time-series. A standard piece of advice for doing in R is something
like the code that follows at the end of the email. I am currently using
an "embed" variant of that code and that piece of advice is out there
too.
But, it occurs to me that for such an easily specified matrix operation
standard R code is really slow. rollapply constantly returns to R
interpreter at each window step for a new lm. All lm is at its heart is
(X^t X)^(-1) * Xy, and if you think about doing that with Rcpp in rolling
window you are just incrementing a counter and peeling off rows (or columns
of X and y) of a particular window size, and following that up with some
matrix multiplication in a loop. The psuedo-code for that Rcpp
practically writes itself and you might want a wrapper of something like:
rolling_lm (y=y, x=x, width=4).
My question is this: has any of the thousands of R packages out there
published anything like that. Rolling window multiple regressions that
stay in C/C++ until the rolling window completes? No sense and writing it
if it exist.
Thanks,
Jeremiah
Standard (slow) advice for "rolling window regression" follows:
set.seed(1)
z <- zoo(matrix(rnorm(10), ncol = 2))
colnames(z) <- c("y", "x")
## rolling regression of width 4
rollapply(z, width = 4,
function(x) coef(lm(y ~ x, data = as.data.frame(x))),
by.column = FALSE, align = "right")
## result is identical to
coef(lm(y ~ x, data = z[1:4,]))
coef(lm(y ~ x, data = z[2:5,]))
[[alternative HTML version deleted]]
Just replacing lm with a faster version would speed it up. Try lm.fit or even faster is fastLm in the RcppArmadillo package. On Thu, Jul 21, 2016 at 2:02 PM, jeremiah rounds <roundsjeremiah at gmail.com> wrote:> Hi, > > A not unusual task is performing a multiple regression in a rolling window > on a time-series. A standard piece of advice for doing in R is something > like the code that follows at the end of the email. I am currently using > an "embed" variant of that code and that piece of advice is out there too. > > But, it occurs to me that for such an easily specified matrix operation > standard R code is really slow. rollapply constantly returns to R > interpreter at each window step for a new lm. All lm is at its heart is > (X^t X)^(-1) * Xy, and if you think about doing that with Rcpp in rolling > window you are just incrementing a counter and peeling off rows (or columns > of X and y) of a particular window size, and following that up with some > matrix multiplication in a loop. The psuedo-code for that Rcpp > practically writes itself and you might want a wrapper of something like: > rolling_lm (y=y, x=x, width=4). > > My question is this: has any of the thousands of R packages out there > published anything like that. Rolling window multiple regressions that > stay in C/C++ until the rolling window completes? No sense and writing it > if it exist. > > > Thanks, > Jeremiah > > Standard (slow) advice for "rolling window regression" follows: > > > set.seed(1) > z <- zoo(matrix(rnorm(10), ncol = 2)) > colnames(z) <- c("y", "x") > > ## rolling regression of width 4 > rollapply(z, width = 4, > function(x) coef(lm(y ~ x, data = as.data.frame(x))), > by.column = FALSE, align = "right") > > ## result is identical to > coef(lm(y ~ x, data = z[1:4,])) > coef(lm(y ~ x, data = z[2:5,])) > > [[alternative HTML version deleted]] > > ______________________________________________ > R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code.-- Statistics & Software Consulting GKX Group, GKX Associates Inc. tel: 1-877-GKX-GROUP email: ggrothendieck at gmail.com
Jeremiah, for this purpose there are the "roll" and "RcppRoll" packages. Both use Rcpp and the former also provides rolling lm models. The latter has a generic interface that let's you define your own function. One thing to pay attention to, though, is the numerical reliability. Especially on large time series with relatively short windows there is a good chance of encountering numerically challenging situations. The QR decomposition used by lm is fairly robust while other more straightforward matrix multiplications may not be. This should be kept in mind when writing your own Rcpp code for plugging it into RcppRoll. But I haven't check what the roll package does and how reliable that is... hth, Z On Thu, 21 Jul 2016, jeremiah rounds wrote:> Hi, > > A not unusual task is performing a multiple regression in a rolling window > on a time-series. A standard piece of advice for doing in R is something > like the code that follows at the end of the email. I am currently using > an "embed" variant of that code and that piece of advice is out there too. > > But, it occurs to me that for such an easily specified matrix operation > standard R code is really slow. rollapply constantly returns to R > interpreter at each window step for a new lm. All lm is at its heart is > (X^t X)^(-1) * Xy, and if you think about doing that with Rcpp in rolling > window you are just incrementing a counter and peeling off rows (or columns > of X and y) of a particular window size, and following that up with some > matrix multiplication in a loop. The psuedo-code for that Rcpp > practically writes itself and you might want a wrapper of something like: > rolling_lm (y=y, x=x, width=4). > > My question is this: has any of the thousands of R packages out there > published anything like that. Rolling window multiple regressions that > stay in C/C++ until the rolling window completes? No sense and writing it > if it exist. > > > Thanks, > Jeremiah > > Standard (slow) advice for "rolling window regression" follows: > > > set.seed(1) > z <- zoo(matrix(rnorm(10), ncol = 2)) > colnames(z) <- c("y", "x") > > ## rolling regression of width 4 > rollapply(z, width = 4, > function(x) coef(lm(y ~ x, data = as.data.frame(x))), > by.column = FALSE, align = "right") > > ## result is identical to > coef(lm(y ~ x, data = z[1:4,])) > coef(lm(y ~ x, data = z[2:5,])) > > [[alternative HTML version deleted]] > > ______________________________________________ > R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. >
Thanks all. roll::roll_lm was essentially what I wanted. I think maybe I would prefer it to have options to return a few more things, but it is the coefficients, and the remaining statistics you might want can be calculated fast enough from there. On Thu, Jul 21, 2016 at 12:36 PM, Achim Zeileis <Achim.Zeileis at uibk.ac.at> wrote:> Jeremiah, > > for this purpose there are the "roll" and "RcppRoll" packages. Both use > Rcpp and the former also provides rolling lm models. The latter has a > generic interface that let's you define your own function. > > One thing to pay attention to, though, is the numerical reliability. > Especially on large time series with relatively short windows there is a > good chance of encountering numerically challenging situations. The QR > decomposition used by lm is fairly robust while other more straightforward > matrix multiplications may not be. This should be kept in mind when writing > your own Rcpp code for plugging it into RcppRoll. > > But I haven't check what the roll package does and how reliable that is... > > hth, > Z > > > On Thu, 21 Jul 2016, jeremiah rounds wrote: > > Hi, >> >> A not unusual task is performing a multiple regression in a rolling window >> on a time-series. A standard piece of advice for doing in R is >> something >> like the code that follows at the end of the email. I am currently using >> an "embed" variant of that code and that piece of advice is out there too. >> >> But, it occurs to me that for such an easily specified matrix operation >> standard R code is really slow. rollapply constantly returns to R >> interpreter at each window step for a new lm. All lm is at its heart is >> (X^t X)^(-1) * Xy, and if you think about doing that with Rcpp in rolling >> window you are just incrementing a counter and peeling off rows (or >> columns >> of X and y) of a particular window size, and following that up with some >> matrix multiplication in a loop. The psuedo-code for that Rcpp >> practically writes itself and you might want a wrapper of something like: >> rolling_lm (y=y, x=x, width=4). >> >> My question is this: has any of the thousands of R packages out there >> published anything like that. Rolling window multiple regressions that >> stay in C/C++ until the rolling window completes? No sense and writing it >> if it exist. >> >> >> Thanks, >> Jeremiah >> >> Standard (slow) advice for "rolling window regression" follows: >> >> >> set.seed(1) >> z <- zoo(matrix(rnorm(10), ncol = 2)) >> colnames(z) <- c("y", "x") >> >> ## rolling regression of width 4 >> rollapply(z, width = 4, >> function(x) coef(lm(y ~ x, data = as.data.frame(x))), >> by.column = FALSE, align = "right") >> >> ## result is identical to >> coef(lm(y ~ x, data = z[1:4,])) >> coef(lm(y ~ x, data = z[2:5,])) >> >> [[alternative HTML version deleted]] >> >> ______________________________________________ >> R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see >> https://stat.ethz.ch/mailman/listinfo/r-help >> PLEASE do read the posting guide >> http://www.R-project.org/posting-guide.html >> and provide commented, minimal, self-contained, reproducible code. >> >>[[alternative HTML version deleted]]
This may be useful: Sven Hammarling and Craig Lucas "Updating the QR factorization and the least squares problem" http://eprints.ma.man.ac.uk/1192/01/covered/MIMS_ep2008_111.pdf http://www.maths.manchester.ac.uk/~clucas/updating/ 2016-07-21 20:02 GMT+02:00 jeremiah rounds <roundsjeremiah at gmail.com>:> Hi, > > A not unusual task is performing a multiple regression in a rolling window > on a time-series. A standard piece of advice for doing in R is something > like the code that follows at the end of the email. I am currently using > an "embed" variant of that code and that piece of advice is out there too. > > But, it occurs to me that for such an easily specified matrix operation > standard R code is really slow. rollapply constantly returns to R > interpreter at each window step for a new lm. All lm is at its heart is > (X^t X)^(-1) * Xy, and if you think about doing that with Rcpp in rolling > window you are just incrementing a counter and peeling off rows (or columns > of X and y) of a particular window size, and following that up with some > matrix multiplication in a loop. The psuedo-code for that Rcpp > practically writes itself and you might want a wrapper of something like: > rolling_lm (y=y, x=x, width=4). > > My question is this: has any of the thousands of R packages out there > published anything like that. Rolling window multiple regressions that > stay in C/C++ until the rolling window completes? No sense and writing it > if it exist. > > > Thanks, > Jeremiah > > Standard (slow) advice for "rolling window regression" follows: > > > set.seed(1) > z <- zoo(matrix(rnorm(10), ncol = 2)) > colnames(z) <- c("y", "x") > > ## rolling regression of width 4 > rollapply(z, width = 4, > function(x) coef(lm(y ~ x, data = as.data.frame(x))), > by.column = FALSE, align = "right") > > ## result is identical to > coef(lm(y ~ x, data = z[1:4,])) > coef(lm(y ~ x, data = z[2:5,])) > > [[alternative HTML version deleted]] > > ______________________________________________ > R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code.