如何向量化实现基于历史数据的累积回归？（R代码场景）

嘿，这个滚动（准确说是累积窗口）回归的需求太常见了！原来的循环写法虽然直观，但数据量一大就会慢到让人头疼。我给你准备了两种更高效的替代方案，尤其是第二种递推OLS的方法，速度提升特别明显：

方案一：用data.table简化循环逻辑 #

data.table的分组操作可以帮我们把循环逻辑包装得更简洁，虽然底层还是迭代，但代码可读性更好，而且比原生for循环略快一点：

library(data.table)
set.seed(123) # 加种子保证结果可复现
dat <- data.frame(t=1:100,y=rnorm(100),x1=rnorm(100),x2=rnorm(100))
setDT(dat) # 转成data.table格式

# 生成累积窗口的系数列表
coefs <- dat[, 
             {
               # 对每个行号idx，用1到idx的数据拟合回归
               if (.I >= 11) list(coef = coefficients(lm(y ~ x1 + x2, data = dat[1:.I, ])))
               else list(coef = NA_real_)
             },
             by = .I]$coef

# 把列表转成矩阵，和原始循环结果对齐
reg_dt <- matrix(NA, nrow = nrow(dat), ncol = 3)
reg_dt[11:nrow(dat), ] <- do.call(rbind, coefs[11:nrow(dat)])
colnames(reg_dt) <- c("(Intercept)", "x1", "x2")

方案二：递推OLS（真正的高效向量化实现） #

如果数据量很大（比如上万行），上面的方法还是不够快。这时候可以用OLS的递推公式——利用Sherman-Morrison公式增量更新回归系数，不需要每次重新拟合整个模型，速度能提升一个数量级：

set.seed(123)
dat <- data.frame(t=1:100,y=rnorm(100),x1=rnorm(100),x2=rnorm(100))

# 构造带截距项的设计矩阵
X <- cbind(1, dat$x1, dat$x2)
y <- dat$y

n <- nrow(X)
k <- ncol(X)
# 初始化X'X的逆矩阵和X'y向量
XX_inv <- matrix(0, nrow = k, ncol = k)
Xy <- numeric(k)
beta <- matrix(NA, nrow = n, ncol = k)

# 递推计算每个累积窗口的系数
for (i in 1:n) {
  xi <- X[i, , drop = FALSE]
  yi <- y[i]
  # 用Sherman-Morrison公式更新X'X的逆矩阵
  denom <- 1 + t(xi) %*% XX_inv %*% xi
  XX_inv <- XX_inv + (xi %*% t(xi)) - (XX_inv %*% xi %*% t(xi) %*% XX_inv) / denom
  # 更新X'y向量
  Xy <- Xy + xi * yi
  # 计算当前窗口的系数
  beta[i, ] <- XX_inv %*% Xy
}

# 前10行设为NA，和原始循环逻辑一致
beta[1:10, ] <- NA
colnames(beta) <- c("(Intercept)", "x1", "x2")

验证结果一致性 #

你可以用下面的代码验证两种方案和原始循环的结果是否一致（浮点误差范围内完全相同）：

# 原始循环的结果
reg <- matrix(rep(NA,3*nrow(dat)),ncol=3)
for(i in 11:nrow(dat)){ 
  reg[i,] <- coefficients(lm(y ~ x1 + x2, data=dat[1:i,])) 
}
colnames(reg) <- c("(Intercept)", "x1", "x2")

all.equal(reg[11:100, ], reg_dt[11:100, ]) # TRUE
all.equal(reg[11:100, ], beta[11:100, ])   # TRUE

内容的提问来源于stack exchange，提问作者bumblebee