2012-08-08 69 views
5

我有一个问题找到最有效的方法来计算多列xts对象的滚动线性回归。我已经在这里搜索和阅读了几个以前的问题在stackoverflow。多列滚动回归

这个question and answer接近但在我看来是不够的,因为我想计算多个回归,因变量在所有回归中保持不变。我试图再现与随机数据的示例:为了随着时间的推移,每个因子多个变量(残差系数等)存储创建

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个观察滚动回归。由于fastLm不计算标准化残差,因此我必须将30个残差存储在临时对象中以便将它们标准化。

如果xts对象中的列(因子)数量增加到大约100 - 1000列,则循环非常缓慢并且变得非常麻烦,这将需要一个永恒。我希望有一个更高效的代码来创建大型数据集上的滚动回归。

+0

通过不运行回归两次,我可以将它编辑成您的问题。 – 2012-08-08 21:25:02

+0

当然可以!这在欧洲很晚。谢谢约书亚。 该变化使性能提高了2-2.5倍。但是,您是否认为该代码对2500个每日观察数据集和1000个因子具有足够的性能? 或者您是否知道使用rollapply与上述方法相比有哪些性能提升? 我想如果数据集变得非常大,你必须应用递归最小二乘滤波器或相关的东西 - 对此的任何想法? – 2012-08-08 21:42:09

回答

8

如果你下降到线性回归的数学水平,它应该是非常快的。如果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的任何错误应该给你的标准化残差。请随时修复此解决方案以纠正任何错误。

+0

+1为直接的方法。 – ricardo 2012-08-27 20:01:26