如何用R代码提取使分位数回归rho值最大化的权重w及优化建议

我正在针对tau序列（0.01到0.99，步长0.01）进行分位数回归，处理1000个响应变量y（公式为y_it = w^T_i * r_t，i=1,...,1000）。每个tau值对应1000个y的分位数回归结果，最终得到1000×99的rho值矩阵（每列对应一个tau）。现有R代码如下：

现有代码 #

生成响应变量y #

# 计算y观测值
all_yts1 <- matrix(nrow = nrow(return),
                   ncol = 1000 )
for( i in 1:1000){
  all_yts1[,i] <- rowSums(all_wts1[i,] * return) 
}

计算分位数回归rho值矩阵 #

tau <- seq(0.01, .99, by = .01)
quant.matrix.reg3 <- matrix(nrow = np1, ncol = length(tau) )
for( i in 1:np1){
  quant.matrix.reg3[i,] <- rq(all_yts1[,i][-1] ~ x.reg[,3], tau)$rho
  colnames(quant.matrix.reg3) <- tau
}

找出每列rho最大值对应的行号 #

opt.w3.alltaus <- matrix(ncol = length(tau) )
for(j in 1: length(tau)){
  opt.w3.alltaus[,j] <- which(quant.matrix.reg3[,j] == max(quant.matrix.reg3[,j]), arr.ind = T) 
}

需求与问题 #

1. 自动提取每个tau对应的最大rho值行号对应的w值（all_wts1中对应行），避免手动提取的繁琐。
2. 优化现有代码，提升运行效率。

解决方案与优化建议 #

1. 自动提取对应w值 #

直接通过行号索引从all_wts1中提取，无需循环：

# 提取每个tau对应的最优w值，结果为99行（对应每个tau），列数与all_wts1一致
opt_w_matrix <- all_wts1[as.vector(opt.w3.alltaus), ]
# 为行命名对应tau值，便于查看
rownames(opt_w_matrix) <- tau

2. 代码优化建议 #

（1）向量化生成all_yts1

替换循环，用矩阵乘法直接计算，大幅提升速度：

# 原循环等价于矩阵内积运算，直接用矩阵乘法替代
all_yts1 <- return %*% t(all_wts1)

解释：return为T×N矩阵（T为时间点，N为因子数），all_wts1为1000×N矩阵，两者运算后直接得到T×1000的all_yts1矩阵，与原循环结果完全一致。

（2）并行化分位数回归计算

针对1000个y变量的分位数回归，用并行计算加速：

library(quantreg)
library(doParallel)

# 根据机器核心数设置并行集群
cl <- makeCluster(detectCores() - 1)
registerDoParallel(cl)

tau <- seq(0.01, .99, by = .01)
np1 <- ncol(all_yts1) # 对应1000个y变量

# 并行计算每个y的rho值，合并为矩阵
quant.matrix.reg3 <- foreach(i = 1:np1, .combine = rbind) %dopar% {
  rq(all_yts1[,i][-1] ~ x.reg[,3], tau)$rho
}
colnames(quant.matrix.reg3) <- tau

# 关闭并行集群
stopCluster(cl)

多核心并行可显著减少计算时间，尤其适合大样本场景。

（3）简化行号查找逻辑

用apply函数替代循环，代码更简洁：

# 按列查找rho最大值对应的行号，若有多个最大值取第一个
opt.w3.alltaus <- apply(quant.matrix.reg3, 2, function(col) which(col == max(col))[1])
# 若需保留所有最大值行号，去掉末尾的[1]即可

优化后完整代码示例 #

library(quantreg)
library(doParallel)

# 1. 向量化生成响应变量y
all_yts1 <- return %*% t(all_wts1)

# 2. 并行计算分位数回归rho矩阵
tau <- seq(0.01, .99, by = .01)
np1 <- ncol(all_yts1)

cl <- makeCluster(detectCores() - 1)
registerDoParallel(cl)

quant.matrix.reg3 <- foreach(i = 1:np1, .combine = rbind) %dopar% {
  rq(all_yts1[,i][-1] ~ x.reg[,3], tau)$rho
}
colnames(quant.matrix.reg3) <- tau

stopCluster(cl)

# 3. 查找每列最大值行号
opt.w3.alltaus <- apply(quant.matrix.reg3, 2, function(col) which(col == max(col))[1])

# 4. 自动提取对应w值
opt_w_matrix <- all_wts1[opt.w3.alltaus, ]
rownames(opt_w_matrix) <- tau

内容的提问来源于stack exchange，提问作者A.F.R.S2022