多列滚动回归

Question

我遇到了一个问题，需要寻找一种计算多列xts对象中滚动线性回归最有效的方法。我在StackOverflow上搜索并阅读了多个相关问题。

这个问题和答案接近我的需求，但是不够完美，因为我想要计算多个回归，其中自变量在所有回归中都保持不变。我尝试使用随机数据复制了一个示例：

require(xts)
require(RcppArmadillo)  # Load libraries

data <- matrix(sample(1:10000, 1500), 1500, 5, byrow = TRUE)  # Random data
data[1000:1500, 2] <- NA  # insert NAs to make it more similar to true data
data <- xts(data, order.by = as.Date(1:1500, origin = "2000-01-01"))

NR <- nrow(data)  # number of observations
NC <- ncol(data)  # number of factors
obs <- 30  # required number of observations for rolling regression analysis
info.names <- c("res", "coef")

info <- array(NA, dim = c(NR, length(info.names), NC))
colnames(info) <- info.names

该数组被创建用于按时间和因素存储多个变量（残差、系数等）。

loop.begin.time <- Sys.time()

for (j in 2:NC) {
  cat(paste("Processing residuals for factor:", j), "\n")
  for (i in obs:NR) {
    regression.temp <- fastLm(data[i:(i-(obs-1)), j] ~ data[i:(i-(obs-1)), 1])
    residuals.temp <- regression.temp$residuals
    info[i, "res", j] <- round(residuals.temp[1] / sd(residuals.temp), 4)
    info[i, "coef", j] <- regression.temp$coefficients[2]
  } 
}

loop.end.time <- Sys.time()
print(loop.end.time - loop.begin.time)  # prints the loop runtime

这个循环的想法是每次将data [,1]作为因变量（因子）与其他因子之一进行30次观测滚动回归。我必须将这30个残差存储在临时对象中，以便将它们标准化，因为fastLm不会计算标准化残差。

如果xts对象中的列数（因子）增加到约100，则循环非常缓慢，并且当列数增加到1000时，处理速度就会变得非常缓慢。希望有更高效的代码来创建大型数据集的滚动回归。

- P. Garnry

---

你可以通过不运行两次回归来使它快2倍...我已经将其编辑到你的问题中了。 - Joshua Ulrich

当然可以！现在欧洲已经很晚了。谢谢你，乔舒亚。这个改变提高了2-2.5倍的性能。但是，您认为这段代码对于每天2500个观测值和约1000个因素的数据集来说，性能是否足够？或者，您是否知道使用rollapply与上述方法相比，是否会提高性能？我猜如果数据集变得非常大，您必须应用递归最小二乘滤波器或类似的东西 - 对此有什么想法吗？ - P. Garnry

Answer 1

如果您降到线性回归的数学层面，那么翻译应该会很快完成。如果X是自变量，Y是因变量，那么系数为： 我有点困惑哪个变量是自变量，哪个是因变量，但希望下面的类似问题的解决方案也可以帮助您。

在下面的示例中，我使用了1000个变量而不是原始的5个变量，并且没有引入任何NA值。

require(xts)

data <- matrix(sample(1:10000, 1500000, replace=T), 1500, 1000, byrow = TRUE)  # Random data
data <- xts(data, order.by = as.Date(1:1500, origin = "2000-01-01"))

NR <- nrow(data)  # number of observations
NC <- ncol(data)  # number of factors
obs <- 30  # required number of observations for rolling regression analysis

现在，我们可以使用Joshua的TTR软件包来计算系数。

library(TTR)

loop.begin.time <- Sys.time()

in.dep.var <- data[,1]
xx <- TTR::runSum(in.dep.var*in.dep.var, obs)
coeffs <- do.call(cbind, lapply(data, function(z) {
    xy <- TTR::runSum(z * in.dep.var, obs)
    xy/xx
}))

loop.end.time <- Sys.time()

print(loop.end.time - loop.begin.time)  # prints the loop runtime

3.934461 秒的时间差

res.array = array(NA, dim=c(NC, NR, obs))
for(z in seq(obs)) {
  res.array[,,z] = coredata(data - lag.xts(coeffs, z-1) * as.numeric(in.dep.var))
}
res.sd <- apply(res.array, c(1,2), function(z) z / sd(z))

如果我没有在索引中犯任何错误，res.sd 应该会给你标准化残差。请随意修正此解决方案以纠正任何错误。

- Sameer

Answer 2

这里有一个更快的方法，使用rollRegres包来完成它

library(xts)
library(RcppArmadillo)

#####
# simulate data
set.seed(50554709)
data <- matrix(sample(1:10000, 1500), 1500, 5, byrow = TRUE)  # Random data
# data[1000:1500, 2] <- NA # only focus on the parts that are computed
data <- xts(data, order.by = as.Date(1:1500, origin = "2000-01-01"))

#####
# setup for solution in OP
NR <- nrow(data)
NC <- ncol(data)
obs <- 30L
info.names <- c("res", "coef")

info <- array(NA, dim = c(NR, length(info.names), NC))
colnames(info) <- info.names

#####
# solve with rollRegres
library(rollRegres)

loop.begin.time <- Sys.time()

X <- cbind(1, drop(data[, 1]))
out <- lapply(2:NC, function(j){
  fit <- roll_regres.fit(
    y = data[, j], x = X, width = obs, do_compute = c("sigmas"))

  # are you sure you want the residual of the first and not the last
  # observation in each window?
  idx <- 1:(nrow(data) - obs + 1L)
  idx_tail <- idx + obs - 1L
  resids <- c(rep(NA_real_, obs - 1L),
                  data[idx, j] - rowSums(fit$coefs[idx_tail, ] * X[idx, ]))

  # the package uses the unbaised estimator so we have to time by this factor
  # to get the same
  sds <-  fit$sigmas * sqrt((obs - 2L) / (obs - 1L))

  unclass(cbind(coef = fit$coefs[, 2L], res = drop(round(resids / sds, 4))))
})

loop.end.time <- Sys.time()
print(loop.end.time - loop.begin.time)
#R Time difference of 0.03123808 secs

#####
# solve with original method
loop.begin.time <- Sys.time()

for (j in 2:NC) {
  cat(paste("Processing residuals for factor:", j), "\n")
  for (i in obs:NR) {
    regression.temp <- fastLm(data[i:(i-(obs-1)), j] ~ data[i:(i-(obs-1)), 1])
    residuals.temp <- regression.temp$residuals
    info[i, "res", j] <- round(residuals.temp[1] / sd(residuals.temp), 4)
    info[i, "coef", j] <- regression.temp$coefficients[2]
  }
}
#R Processing residuals for factor: 2
#R Processing residuals for factor: 3
#R Processing residuals for factor: 4
#R Processing residuals for factor: 5

loop.end.time <- Sys.time()
print(loop.end.time - loop.begin.time)  # prints the loop runtime
#R Time difference of 7.554767 secs

#####
# check that results are the same
all.equal(info[, "coef", 2L], out[[1]][, "coef"])
#R [1] TRUE
all.equal(info[, "res" , 2L], out[[1]][, "res"])
#R [1] TRUE

all.equal(info[, "coef", 3L], out[[2]][, "coef"])
#R [1] TRUE
all.equal(info[, "res" , 3L], out[[2]][, "res"])
#R [1] TRUE

all.equal(info[, "coef", 4L], out[[3]][, "coef"])
#R [1] TRUE
all.equal(info[, "res" , 4L], out[[3]][, "res"])
#R [1] TRUE

all.equal(info[, "coef", 5L], out[[4]][, "coef"])
#R [1] TRUE
all.equal(info[, "res" , 5L], out[[4]][, "res"])
#R [1] TRUE

请注意上述解决方案中的此注释。

# are you sure you want the residual of the first and not the last
# observation in each window?

---

这里是与Sameer的回答的比较

library(rollRegres)
require(xts)

data <- matrix(sample(1:10000, 1500000, replace=T), 1500, 1000, byrow = TRUE)  # Random data
data <- xts(data, order.by = as.Date(1:1500, origin = "2000-01-01"))

NR <- nrow(data)  # number of observations
NC <- ncol(data)  # number of factors
obs <- 30  # required number of observations for rolling regression analysis

loop.begin.time <- Sys.time()

X <- cbind(1, drop(data[, 1]))
out <- lapply(2:NC, function(j){
  fit <- roll_regres.fit(
    y = data[, j], x = X, width = obs, do_compute = c("sigmas"))

  # are you sure you want the residual of the first and not the last
  # observation in each window?
  idx <- 1:(nrow(data) - obs + 1L)
  idx_tail <- idx + obs - 1L
  resids <- c(rep(NA_real_, obs - 1L),
              data[idx, j] - rowSums(fit$coefs[idx_tail, ] * X[idx, ]))

  # the package uses the unbaised estimator so we have to time by this factor
  # to get the same
  sds <-  fit$sigmas * sqrt((obs - 2L) / (obs - 1L))

  unclass(cbind(coef = fit$coefs[, 2L], res = drop(round(resids / sds, 4))))
})

loop.end.time <- Sys.time()
print(loop.end.time - loop.begin.time)
#R Time difference of 0.9019711 secs

时间包括用于计算标准化残差的时间。

- Benjamin Christoffersen