在R中对XTS序列应用滚动窗口回归

Question

在R中对XTS序列应用滚动窗口回归

10

我有一个由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是一个更简单的函数，例如平均值，它似乎可以正常工作：

> 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的回归权重。

- Thomas Browne

2个回答

3

新答案

G. Grothendieck的回答是正确的，但您可以使用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（或其他用于更新 R 的函数）获取 dchud 和 dchdd 子例程。

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）。

- Benjamin Christoffersen

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- G. Grothendieck · Accepted Answer

10

这里存在几个问题：

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]))

- G. Grothendieck

太棒了，谢谢。是的，我也经过很多尝试才解决了它。我真傻。by.column = FALSE 当然没问题！非常感谢。顺便说一下，我刚刚在阅读你的zoo文档。非常好的东西。我想，在rollapply方面有点令人困惑的是，虽然lm()适用于整个xts，但它不适用于由rollapply()返回的部分。人们可以合理地期望rollapply返回另一个仍然可以在lm()下工作的xts，或者我错过了什么？对于by.column FALSE，我承认错误。没有任何借口... - Thomas Browne

错过的是 rollapply 不是 xts 的一部分，而是 zoo 的一部分，并且它的调度是 rollapply.zoo。 - G. Grothendieck

谢谢您澄清这个问题。然而：

fxr <- zoo(fxr) class(fxr) [1] "zoo" rollapply(fxr, 62, function(x) coef(lm(USDZAR ~ x, data = x)), by.column = FALSE) Error in model.frame.default(formula = USDZAR ~ x, data = x, drop.unused.levels = TRUE) : 'data'必须是数据框，而不是矩阵或数组。

所以我们仍然有这个问题。我明白...R在这方面确实存在很多问题，但是我们现在的问题是lm可以作用于整个zoo对象，但无法作用于rollapply的子集。 - Thomas Browne

这只是用户错误。没有理由认为lm可以与其文档中所述的以外的任何东西一起使用。此外，lm在数据参数上不是通用的（也许您觉得它应该是），因此没有理由认为特定的软件包可以扩展它，尽管确实存在两个软件包--dyn和dynlm--它们将允许您使用zoo对象进行线性回归（dyn还允许进行多种其他类型的回归），但不能使用矩阵。如果您确实想使用矩阵，则lm.fit可以实现（如我所提到的）。 - G. Grothendieck

非常感谢。lm.fit 在其行为上似乎完全一致，因此我将使用它。 - Thomas Browne

显示剩余2条评论