滚动回归多列

2 回答

• 9

  如果你深入到线性回归的数学水平，它应该很快 . 如果X是自变量而Y是因变量 . 系数由下式给出

  Beta = inv(t(X) %*% X) %*% (t(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 应该给你标准化的残差 . 请随时修复此解决方案以纠正任何错误 .

  回复于 2024-04-19T22:43:42+08:00
• 0

  使用 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's answer的比较

  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
  

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

  回复于 2024-04-19T22:43:42+08:00