2012-08-08 20 views
5

Tengo un problema para encontrar la forma más eficiente de calcular una regresión lineal continua sobre un objeto xts con múltiples columnas. He buscado y leído varias preguntas anteriores aquí en stackoverflow.Regresión continua sobre múltiples columnas

Este question and answer se acerca, pero no lo suficiente, en mi opinión, ya que quiero calcular regresiones múltiples con la variable dependiente sin cambios en todas las regresiones. He intentado reproducir un ejemplo con datos aleatorios:

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 matriz se creó con el fin de almacenar múltiples variables (residuos, coeficientes etc.) con el tiempo y por factor.

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 

Como el bucle de muestra la idea es ejecutar una regresión de rodadura 30 observaciones con data[, 1] como la variable dependiente (factor) cada vez en contra de uno de los otros factores. Tengo que almacenar los 30 residuales en un objeto temporal para estandarizarlos ya que fastLm no calcula los residuos estandarizados.

El ciclo es extremadamente lento y se vuelve engorroso si el número de columnas (factores) en el objeto xts aumenta a alrededor de 100 - 1.000 columnas tardarían una eternidad. Espero que uno tenga un código más eficiente para crear regresiones continuas sobre un gran conjunto de datos.

+0

Puede hacer que sea 2 veces más rápido al no ejecutar la regresión dos veces ... que he editado en su pregunta. –

+0

Sí, por supuesto! Es tarde aquí en Europa. Gracias Joshua. El cambio ha aumentado el rendimiento en 2-2.5x. Sin embargo, ¿considera que este código tiene un rendimiento adecuado para un conjunto de datos de 2500 observaciones diarias y alrededor de 1.000 factores? ¿O conoce alguna ganancia en el rendimiento al usar rollapply en comparación con el enfoque anterior? Supongo que si el conjunto de datos se vuelve muy grande, tiene que aplicar filtro recursivo de mínimos cuadrados o algo relacionado, ¿tiene alguna idea al respecto? –

Respuesta

8

Debería ser bastante rápido si baja al nivel de las matemáticas de la regresión lineal. Si X es la variable independiente e Y es la variable dependiente. Los coeficientes están dados por

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

Estoy un poco confundido acerca de qué variable que desea ser el dependiente y que es el independiente pero esperemos que la solución por debajo de un problema similar le ayudará también.

En el siguiente ejemplo tomo 1000 variables en vez de las 5 originales y no introduzco ninguna 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 

Ahora podemos calcular los coeficientes utilizando el paquete TTR de 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 

diferencia de tiempo de 3.934461 segundos

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

Si no he cometido errores en la indexación res.sd debe darle los residuos estandarizados. Por favor, siéntase libre de arreglar esta solución para corregir cualquier error.

+0

+1 para el enfoque directo. – ricardo

Cuestiones relacionadas