【问题标题】：Applying a rolling window regression to an XTS series in R将滚动窗口回归应用于 R 中的 XTS 系列
【发布时间】：2012-03-10 04:47:24
【问题描述】：

对于 5 个货币对，我有 1033 个每日回报点的 xts，我想在这些货币对上运行滚动窗口回归，但 rollapply 不适用于我定义的使用 lm() 的函数。这是我的数据：

> 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

我可以轻松地对整个数据集运行 lm 以针对其他货币对对 USDZAR 进行建模：

> 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

但是我想运行一个滚动的 62 天窗口来了解这些系数随时间的演变，所以我创建了一个函数 dolm 来执行此操作：

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

但是，当我对此运行 rollapply 时，我得到以下信息：

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

即使 dolm(fxr) 本身工作正常：

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

这里发生了什么？如果 dolm 是一个更简单的函数（例如 mean），它似乎可以正常工作：

> 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

非常感谢任何帮助。基本上我想要的是在滚动的 62 天窗口内获得 USDZAR ~ USDEUR + USDGBP + USDCHF + USDCAD 回归的权重。

【问题讨论】：

标签： r regression xts

【解决方案1】：

新答案

G. Grothendieck's answer 是正确的，但您可以使用 rollRegres 包更快地做到这一点，如下例所示（roll_regres.fit 函数快约 118 倍）

# 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)
library(rollRegres)
all.equal(
  rollapply(Z, wdth, FUN = dolm,
            by.column = FALSE,  align = "right", fill = NA_real_),
  roll_regres.fit(X, y, wdth)$coefs,
  check.attributes = FALSE)
#R [1] TRUE

# benchmark
library(compiler)
dolm <- cmpfun(dolm)

microbenchmark::microbenchmark(
  newnew = roll_regres.fit(X, y, wdth),
  prev   = rollapply(Z, wdth, FUN = dolm,
                     by.column = FALSE,  align = "right", fill = NA_real_),
  times = 10)
#R Unit: microseconds
#R expr        min         lq       mean     median         uq        max neval
#R newnew    884.938    950.914   1026.134   1025.581   1057.581   1242.075    10
#R   prev 111057.822 111903.649 118867.761 116857.726 122087.160 141362.229    10

如果您想改用 R 公式，也可以使用包中的 roll_regres 函数。

旧答案

第三种选择是在 QR 分解中更新 R 矩阵，如下面的代码所示。您可以通过在 C++ 中执行此操作来加快速度，但您需要来自 LINPACK 的 dchud 和 dchdd 子例程（或其他更新 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 R 
      R <- .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

上述解决方案要求您首先形成model.matrix 和model.response，但这只是在调用roll_coef 之前的三个调用（model.frame 的一个额外调用）。

【讨论】：

【解决方案2】：

这里有几个问题：

rollapply 传递一个矩阵，但 lm 需要一个 data.frame。
rollapply 将函数分别应用于每一列，除非我们指定by.column=FALSE。
您可能希望也可能不希望结果是与日期右对齐，但如果您使用 rollapplyr ：

1) 结合以上我们有：

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

2) 上面dolm 中lm 的替代方法是使用lm.fit，它直接与矩阵一起工作并且速度更快：

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

【讨论】：

非常感谢。是的，我只是在玩了很多之后才解决了这个问题。傻我。 by.column = FALSE 当然！非常感谢。顺便说一句，正在阅读您的动物园文档。好东西。我猜 rollapply 有点令人困惑的地方是，虽然 lm() 可以在整个 xts 上工作，但它不能在 rollapply() 返回的部分上工作。人们可以合理地期望 rollapply 返回另一个在 lm() 下仍然可以工作的 xts 或者我错过了什么？虽然在 by.column FALSE 上有过失。没有任何借口......
遗漏的是 rollapply 不是 xts 的一部分，但它是 zoo 及其调度 rollapply.zoo 的一部分。
感谢您澄清这一点。然而： > fxr class(fxr) [1] "zoo" > rollapply(fxr, 62, function(x) coef(lm(USDZAR ~ x, data = x)), by.column = FALSE) model.frame.default(formula = USDZAR ~ x, data = x, drop.unused.levels = TRUE) 中的错误：'data' 必须是 data.frame，而不是矩阵或数组所以我们仍然有这个问题。我明白... R 有很多这样的问题，但仍然存在。我们这里的 lm 适用于整个 zoo 对象，但不适用于它的 rollapply 子集。
这只是用户错误。没有理由认为lm 可以与除文档记录之外的任何东西一起使用。 lm 在 data 参数中也不是通用的（也许你觉得它应该是）所以没有理由认为特定的包可以扩展它，尽管确实存在两个包 - dyn 和 dynlm - 每个都允许您可以使用 zoo 对象而不是矩阵进行线性回归（dyn 还允许许多其他类型的回归）。如果您确实想使用矩阵，那么 lm.fit 会这样做（如我的回复中所述）。
谢谢。 lm.fit 的行为似乎完全一致，所以我将使用它。