Hay varias maneras de hacer esto en R. He aquí dos ejemplos utilizando el "Seatbelts" time series dataset in the datasets package que viene con R.
La función arima()
viene en el paquete: las estadísticas que se incluyen con R. La función toma un argumento de el formulario order=c(p, d, q)
donde puede especificar el orden del componente autoregresivo, integrado y promedio móvil. En su pregunta, sugiere que desea crear un modelo AR (1) para corregir la autocorrelación de primer orden en los errores y eso es todo. Podemos hacer eso con el siguiente comando:
arima(Seatbelts[,"drivers"], order=c(1,0,0),
xreg=Seatbelts[,c("kms", "PetrolPrice", "law")])
El valor para la orden especifica que queremos un (1) modelo AR. El componente xreg debe ser una serie de otras X que queremos agregar como parte de una regresión. La salida se ve un poco como la salida de summary.lm()
puesta de lado.
Otro proceso alternativo podría ser más familiar para la forma en que usted ha ajustado los modelos de regresión es usar gls()
en el nlme package. El siguiente código convierte el objeto de series de tiempo del cinturón de seguridad en una trama de datos y luego extrae y añade una nueva columna (t) que es sólo un contador en el objeto de series de tiempo Ordenada:
Seatbelts.df <- data.frame(Seatbelts)
Seatbelts.df$t <- 1:(dim(Seatbelts.df)[1])
Las dos líneas anteriores sólo son poniendo los datos en forma. Como la función arima()
está diseñada para series de tiempo, puede leer objetos de series de tiempo más fácilmente. Para ajustar el modelo con nlme usted entonces ejecute:
library(nlme)
m <- gls(drivers ~ kms + PetrolPrice + law,
data=Seatbelts.df,
correlation=corARMA(p=1, q=0, form=~t))
summary(m)
La línea que comienza con "correlación" es la forma en que se pasa en la estructura de correlación ARMA de GLS. Los resultados no serán exactamente los mismos porque arima()
usa la máxima probabilidad de estimar modelos y gls()
usa la máxima verosimilitud restringida de forma predeterminada. Si agrega method="ML"
a la llamada al gls()
, obtendrá las estimaciones idénticas que obtuvo con la función ARIMA anterior.
El paquete GLMMarp se eliminó del repositorio CRAN. ¿Conoces algún otro paquete que logre esto? – TMS
Todos los paquetes que conozco están listados en http://cran.r-project.org/web/views/TimeSeries.html –