rendertext - ¿Agregando variables retrasadas a un modelo lm?

tags$style shiny

Question

Estoy usando lm en una serie de tiempo, que funciona bastante bien en realidad, y es super súper rápido.

Digamos que mi modelo es:

> formula <- y ~ x

Entreno esto en un set de entrenamiento:

> train <- data.frame( x = seq(1,3), y = c(2,1,4) )
> model <- lm( formula, train )

... y puedo hacer predicciones para nuevos datos:

> test <- data.frame( x = seq(4,6) )
> test$y <- predict( model, newdata = test )
> test
x y
1 4 4.333333
2 5 5.333333
3 6 6.333333

Esto funciona muy bien, y es muy rápido.

Quiero agregar variables retrasadas al modelo. Ahora, puedo hacer esto aumentando mi conjunto de entrenamiento original:

> train$y_1 <- c(0,train$y[1:nrow(train)-1])
> train
x y y_1
1 1 2 0
2 2 1 2
3 3 4 1

actualizar la fórmula:

formula <- y ~ x * y_1

... y el entrenamiento funcionará bien:

> model <- lm( formula, train )
> # no errors here

Sin embargo, el problema es que no hay manera de usar ''predecir'', porque no hay manera de llenar y_1 en un conjunto de prueba de manera discontinua.

Ahora, para muchas otras cosas de regresión, hay formas muy convenientes de expresarlas en la fórmula, como poly(x,2) etc., y estas funcionan directamente utilizando los datos de pruebas y entrenamiento no modificados.

Entonces, me pregunto si hay alguna forma de expresar las variables retrasadas en la fórmula, ¿para que se pueda usar la predict ? Idealmente:

formula <- y ~ x * lag(y,-1)
model <- lm( formula, train )
test$y <- predict( model, newdata = test )

... ¿sin tener que aumentar (no estoy seguro de si esa es la palabra correcta) los conjuntos de datos de entrenamiento y prueba, y solo poder usar la predict directamente?

Answer 1

Aquí hay un pensamiento:

¿Por qué no creas un nuevo marco de datos? Rellene un marco de datos con los regresores que necesita. Podría tener columnas como L1, L2, ..., Lp para todos los retrasos de cualquier variable que desee y, luego, podrá usar sus funciones exactamente como lo haría para un tipo de regresión de sección transversal.

Debido a que no tendrá que operar con sus datos cada vez que llame a las funciones de ajuste y predicción, sino que habrá transformado los datos una vez, será considerablemente más rápido. Sé que Eviews y Stata proporcionan operadores rezagados. Es cierto que hay algo de conveniencia para ello. Pero también es ineficiente si no necesita que todo funcione como el cálculo ''lm''. Si tiene que realizar cientos de miles de iteraciones y solo necesita el pronóstico, o el pronóstico y el valor de los criterios de información como BIC o AIC, puede mejorar la velocidad de ''lm'' evitando realizar cálculos que no hará. uso: solo escribe un estimador OLS en una función y listo.

Answer 2

Eche un vistazo a, por ejemplo, el paquete dynlm que le da operadores de retardo. Más generalmente, las Vistas de tareas en econometría y series temporales tendrán mucho más que ver.

Aquí está el comienzo de sus ejemplos: un retraso de uno y doce meses:

R> data("UKDriverDeaths", package = "datasets")
R> uk <- log10(UKDriverDeaths)
R> dfm <- dynlm(uk ~ L(uk, 1) + L(uk, 12))
R> dfm
Time series regression with "ts" data:
Start = 1970(1), End = 1984(12)
Call:
dynlm(formula = uk ~ L(uk, 1) + L(uk, 12))
Coefficients:
(Intercept) L(uk, 1) L(uk, 12)
0.183 0.431 0.511
R> 

Answer 3

Prueba la función ARIMA. El parámetro AR es para auto-regresivo, lo que significa y retrasado. xreg = le permite agregar otras variables X. Puedes obtener predicciones con predict.ARIMA.

Answer 4

Siguiendo la sugerencia de Dirk en dynlm , no pude averiguar cómo predecir, pero la búsqueda de eso me llevó al paquete dyn través de https://stats.stackexchange.com/questions/6758/1-step-ahead-predictions-with-dynlm-r-package

Luego, después de varias horas de experimentación, se me ocurrió la siguiente función para manejar la predicción. Hubo un buen número de ''gotcha'' en camino, por ejemplo, parece que no puedes encontrar rbind serie de tiempo, y el resultado de predicción se compensa con el start y un montón de cosas así, por lo que creo que esta respuesta aumenta significativamente en comparación con solo nombrando un paquete, aunque he subido la respuesta de Dirk.

Entonces, una solución que funciona es:

• usa el paquete dyn
• Utilice el siguiente método para la predicción

método predictDyn:

# pass in training data, test data,
# it will step through one by one
# need to give dependent var name, so that it can make this into a timeseries
predictDyn <- function( model, train, test, dependentvarname ) {
Ntrain <- nrow(train)
Ntest <- nrow(test)
# can''t rbind ts''s apparently, so convert to numeric first
train[,dependentvarname] <- as.numeric(train[,dependentvarname])
test[,dependentvarname] <- as.numeric(test[,dependentvarname])
testtraindata <- rbind( train, test )
testtraindata[,dependentvarname] <- ts( as.numeric( testtraindata[,dependentvarname] ) )
for( i in 1:Ntest ) {
result <- predict(model,newdata=testtraindata,subset=1:(Ntrain+i-1))
testtraindata[Ntrain+i,dependentvarname] <- result[Ntrain + i + 1 - start(result)][1]
}
return( testtraindata[(Ntrain+1):(Ntrain + Ntest),] )
}

Ejemplo de uso:

library("dyn")
# size of training and test data
N <- 6
predictN <- 10
# create training data, which we can get exact fit on, so we can check the results easily
traindata <- c(1,2)
for( i in 3:N ) { traindata[i] <- 0.5 + 1.3 * traindata[i-2] + 1.7 * traindata[i-1] }
train <- data.frame( y = ts( traindata ), foo = 1)
# create testing data, bunch of NAs
test <- data.frame( y = ts( rep(NA,predictN) ), foo = 1)
# fit a model
model <- dyn$lm( y ~ lag(y,-1) + lag(y,-2), train )
# look at the model, it''s a perfect fit. Nice!
print(model)
test <- predictDyn( model, train, test, "y" )
print(test)
# nice plot
plot(test$y, type=''l'')

Salida:

> model
Call:
lm(formula = dyn(y ~ lag(y, -1) + lag(y, -2)), data = train)
Coefficients:
(Intercept) lag(y, -1) lag(y, -2)
0.5 1.7 1.3
> test
y foo
7 143.2054 1
8 325.6810 1
9 740.3247 1
10 1682.4373 1
11 3823.0656 1
12 8686.8801 1
13 19738.1816 1
14 44848.3528 1
15 101902.3358 1
16 231537.3296 1

Edit: hmmm, esto es super lento aunque. Incluso si limito los datos en el subset a unas pocas filas constantes del conjunto de datos, toma aproximadamente 24 milisegundos por predicción o, para mi tarea, 0.024*7*24*8*20*10/60/60 1.792 hours = 1.792 hours : -O