Application of rolling regression to the XTS series in R

Question

I have xts of 1033 daily return points for 5 currency pairs on which I want to run rolling regression, but rollapply does not work for my specific function that uses lm (). Here are my details:

> head(fxr) USDZAR USDEUR USDGBP USDCHF USDCAD 2007-10-18 -0.005028709 -0.0064079963 -0.003878743 -0.0099537170 -0.0006153215 2007-10-19 -0.001544470 0.0014275520 -0.001842564 0.0023058211 -0.0111410271 2007-10-22 0.010878027 0.0086642116 0.010599365 0.0051899551 0.0173792230 2007-10-23 -0.022783987 -0.0075236355 -0.010804304 -0.0041668499 -0.0144788687 2007-10-24 -0.006561223 0.0008545792 0.001024275 -0.0004261666 0.0049525483 2007-10-25 -0.014788901 -0.0048523001 -0.001434280 -0.0050425302 -0.0046422944 > tail(fxr) USDZAR USDEUR USDGBP USDCHF USDCAD 2012-02-10 0.018619309 0.007548205 0.005526184 0.006348533 0.0067151342 2012-02-13 -0.006449463 -0.001055966 -0.002206810 -0.001638002 -0.0016995755 2012-02-14 0.006320364 0.006843933 0.006605875 0.005992935 0.0007001751 2012-02-15 -0.001666872 0.004319096 -0.001568874 0.003686840 -0.0015009759 2012-02-16 0.006419616 -0.003401364 -0.005194817 -0.002709588 -0.0019044761 2012-02-17 -0.004339687 -0.003675992 -0.003319899 -0.003043481 0.0000000000 

I can easily run lm on it so that the entire dataset simulates USDZAR against other pairs:

 > lm(USDZAR ~ ., data = fxr)$coefficients (Intercept) USDEUR USDGBP USDCHF USDCAD -1.309268e-05 5.575627e-01 1.664283e-01 -1.657206e-01 6.350490e-01 

However, I want to run a maximized 62-day window to get the evolution of these coefficients over time, so I create a dolm function that does this:

 > dolm function(x) { return(lm(USDZAR ~ ., data = x)$coefficients) } 

However, when I run rollapply, I get the following:

 > rollapply(fxr, 62, FUN = dolm) Error in terms.formula(formula, data = data) : '.' in formula and no 'data' argument 

even if dolm (fxr) itself works fine:

 > dolm(fxr) (Intercept) USDEUR USDGBP USDCHF USDCAD -1.309268e-05 5.575627e-01 1.664283e-01 -1.657206e-01 6.350490e-01 

What's going on here? This seems to be normal if dolm is a simpler function, for example, means:

 > dolm <- edit(dolm) > dolm function(x) { return(mean(x)) } > rollapply(fxr, 62, FUN = dolm) USDZAR USDEUR USDGBP USDCHF USDCAD 2007-11-29 -1.766901e-04 -6.899297e-04 6.252596e-04 -1.155952e-03 7.021468e-04 2007-11-30 -1.266130e-04 -6.512204e-04 7.067767e-04 -1.098413e-03 7.247315e-04 2007-12-03 8.949942e-05 -6.406932e-04 6.637066e-04 -1.154806e-03 8.727564e-04 2007-12-04 2.042046e-04 -5.758493e-04 5.497422e-04 -1.116308e-03 7.124593e-04 2007-12-05 7.343586e-04 -4.899982e-04 6.161819e-04 -1.057904e-03 9.915495e-04 

Any help is greatly appreciated. In fact, I want to get the weighting coefficients for the regression of USDZAR ~ USDEUR + USDGBP + USDCHF + USDCAD during the rolling 62-day window.

Answer 1

There are several issues here:

• rollapply passes the matrix, but lm requires data.frame .
• rollapply applies the function to each column separately if we specify by.column=FALSE .
• you may or may not want the result to be right aligned with dates, but if you use rollapplyr :

1) Including the above, we have:

 dolm <- function(x) coef(lm(USDZAR ~ ., data = as.data.frame(x)))) rollapplyr(fxr, 62, dolm, by.column = FALSE) 

2) An alternative to lm in dolm above is to use lm.fit , which works directly with matrices and is also faster:

 dolm <- function(x) coef(lm.fit(cbind(Intercept = 1, x[,-1]), x[,1])) 

Answer 2

The third option is to update the R matrix in the QR decomposition, as is done in the code below. You can speed it up by doing this in C ++, but you will need the dchud and dchdd from LINPACK (or another function to update R)

 library(SamplerCompare) # for LINPACK `chdd` and `chud` roll_coef <- function(X, y, width){ n <- nrow(X) p <- ncol(X) out <- matrix(NA_real_, n, p) is_first <- TRUE i <- width while(i <= n){ if(is_first){ is_first <- FALSE qr. <- qr(X[1:width, ]) R <- qr.R(qr.) # Use X^T for the rest X <- t(X) XtY <- drop(tcrossprod(y[1:width], X[, 1:width])) } else { x_new <- X[, i] x_old <- X[, i - width] # update RR <- .Fortran( "dchud", R, p, p, x_new, 0., 0L, 0L, 0., 0., numeric(p), numeric(p), PACKAGE = "SamplerCompare")[[1]] # downdate R R <- .Fortran( "dchdd", R, p, p, x_old, 0., 0L, 0L, 0., 0., numeric(p), numeric(p), integer(1), PACKAGE = "SamplerCompare")[[1]] # update XtY XtY <- XtY + y[i] * x_new - y[i - width] * x_old } coef. <- .Internal(backsolve(R, XtY, p, TRUE, TRUE)) out[i, ] <- .Internal(backsolve(R, coef., p, TRUE, FALSE)) i <- i + 1 } out } # simulate data set.seed(101) n <- 1000 wdth = 100 X <- matrix(rnorm(10 * n), n, 10) y <- drop(X %*% runif(10)) + rnorm(n) Z <- cbind(y, X) # assign other function dolm <- function(x) coef(lm.fit(x[, -1], x[, 1])) # show that they yield the same library(zoo) all.equal( rollapply(Z, wdth, FUN = dolm, by.column = FALSE, align = "right", fill = NA_real_), roll_coef(X, y, wdth), check.attributes = FALSE) #R> [1] TRUE # benchmark library(compiler) roll_coef <- cmpfun(roll_coef) dolm <- cmpfun(dolm) microbenchmark::microbenchmark( new = roll_coef(X, y, wdth), prev = rollapply(Z, wdth, FUN = dolm, by.column = FALSE, align = "right", fill = NA_real_), times = 10) #R> Unit: milliseconds #R> expr min lq mean median uq max neval cld #R> new 8.631319 9.010579 9.808525 9.659665 9.973741 11.87083 10 a #R> prev 118.257128 121.734860 124.489826 122.882318 127.195410 135.21280 10 b 

To solve the above, you first need to form model.matrix and model.response , but these are just three calls (one additional to model.frame ) before calling roll_coef .