2012-08-08 6 views
5

Ho un problema trovare il modo più efficiente per calcolare una regressione lineare mobile su un oggetto XTS con più colonne. Ho cercato e letto diverse domande precedenti qui su StackOverflow.rotolamento di regressione su più colonne

Questo question and answer si avvicina ma non abbastanza a mio parere, come voglio calcolare regressioni multiple con la variabile dipendente invariato in tutte le regressioni. Ho cercato di riprodurre un esempio con dati casuali:

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 

La matrice viene creata per memorizzare variabili multiple (residui, coefficienti ecc) nel tempo e per ogni fattore.

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 

Come loop mostra l'idea è quella di eseguire una regressione 30 osservazioni rotazione con data[, 1] come variabile dipendente (fattore) ogni volta contro uno degli altri fattori. Devo memorizzare i 30 residui in un oggetto temporaneo per standardizzarli come fastLm non calcola i residui standardizzati.

Il ciclo è estremamente lento e diventa ingombrante se il numero di colonne (fattori) nell'oggetto xts aumenta a circa 100 - 1.000 colonne richiederebbe un'eternità. Spero che uno abbia un codice più efficiente per creare regressioni rolling su un set di dati di grandi dimensioni.

+0

Si potrebbe rendere più 2x più veloce, non eseguendo la regressione due volte ... che ho curato nella sua domanda. –

+0

Sì, certo! È tardi qui in Europa. Grazie Joshua. La modifica ha aumentato le prestazioni di 2-2,5 volte. Tuttavia, ritieni che questo codice abbia prestazioni adeguate per un set di dati di 2500 osservazioni giornaliere e circa 1.000 fattori? o siete a conoscenza di qualsiasi guadagno in prestazioni utilizzando rollapply rispetto all'approccio di cui sopra? Credo che se il set di dati diventa molto grande si deve applicare ricorsiva almeno filtro piazze o qualcosa legato - qualche idea su questo? –

risposta

8

dovrebbe essere abbastanza veloce se si scende al livello della matematica della regressione lineare. Se X è la variabile indipendente e Y è la variabile dipendente. I coefficienti sono dati da

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

Sono un po 'confuso su quale variabile che si desidera essere quello dipendente e che è indipendente ma si spera di risolvere un problema simile di seguito vi aiuterà pure.

Nell'esempio che segue prendo 1000 variabili al posto dell'originale 5 e non creino NA di.

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 

Ora possiamo calcolare i coefficienti utilizzando il pacchetto TTR di Joshua.

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 

differenza di tempo di 3.934461 secondi

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)) 

Se non ho fatto errori nell'indicizzazione res.sd dovrebbe darvi i residui standardizzati. Non esitate a risolvere questa soluzione per correggere eventuali bug.

+0

+1 per l'approccio diretto. – ricardo