2012-08-08 19 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〜1,000列に増えて永遠になる場合は面倒です。大規模なデータセットに対してローリング回帰を作成するためのより効率的なコードがあることを願っています。

+0

私はあなたの質問に編集した回帰を2回実行しないことで2倍速くすることができます。 –

+0

はい、もちろんです!ここヨーロッパで遅いです。ジョシュアありがとう。 変更のパフォーマンスが2-2.5倍向上しました。しかし、このコードは、1日に2500件の観測データと約1,000件のデータセットに対して適切なパフォーマンスを持っていると考えていますか? 上記の方法と比較してrollapplyを使用してパフォーマンスの向上を確認していますか? データセットが非常に大きくなるかどうかは、再帰最小二乗フィルタなどを適用する必要があると思います。 –

答えて

8

線形回帰の数学のレベルに下がるのはかなり速いはずです。 Xが独立変数であり、Yが従属変数である場合。係数は、私はあなたが独立しているが、うまくいけば、以下の同様の問題を解決するだけでなく、あなたを助ける依存し、どちらになりたいどの変数について少し混乱してい

Beta = inv(t(X) %*% X) %*% (t(X) %*% Y)

によって与えられています。

以下の例では、元の5の代わりに1000個の変数を取りますが、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 

ここで、ジョシュアのTTRパッケージを使用して係数を計算することができます。 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)) 

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 

時間差Iはインデックスres.sdにすべてのエラーを行っていない場合は、あなたの標準化残差を与える必要があります。バグ修正のためにこのソリューションを修正してください。

+0

+1直接アプローチ。 – ricardo

関連する問題