¿Cómo ajustar un modelo ARIMAX con R?

Question

Tengo cuatro series temporales diferentes de mediciones horarias:

1. El consumo de calor dentro de una casa
2. La temperatura fuera de la casa
3. La radiación solar
4. La velocidad del viento

Quiero poder predecir el consumo de calor dentro de la casa. Hay una clara tendencia estacional, tanto a nivel anual como diario. Como hay una clara correlación entre las distintas series, quiero ajustarlas mediante un modelo ARIMAX. Esto se puede hacer en R, utilizando la función arimax del paquete TSA.

He intentado leer la documentación de esta función, y leer sobre las funciones de transferencia, pero hasta ahora, mi código:

regParams = ts.union(ts(dayy))
transferParams = ts.union(ts(temp))
model10 = arimax(heat,order=c(2,1,1),seasonal=list(order=c(0,1,1),period=24),xreg=regParams,xtransf=transferParams,transfer=list(c(1,1))
pred10 = predict(model10, newxreg=regParams)

me da:

donde la línea negra son los datos reales medidos, y la línea verde es mi modelo ajustado en comparación. No sólo no es un buen modelo, sino que claramente algo está mal.

Reconozco que mis conocimientos sobre los modelos ARIMAX y las funciones de transferencia son limitados. En la función arimax(), (por lo que he entendido), xtransf es la serie temporal exógena que quiero utilizar (mediante funciones de transferencia) para predecir mi serie temporal principal. Pero, ¿cuál es la diferencia entre xreg y xtransf realmente?

En general, ¿qué he hecho mal? Me gustaría ser capaz de obtener un mejor ajuste que el logrado a partir de lm(heat ~ temp radio viento*tiempo).

Editado: Basándome en algunos de los comentarios, he eliminado la transferencia y he añadido xreg en su lugar:

regParams = ts.union(ts(dayy), ts(temp), ts(time))
model10 = arimax(heat,order=c(2,1,1),seasonal=list(order=c(0,1,1),period=24),xreg=regParams)

donde díay es el "número de día del año", y hora es la hora del día. Temp es de nuevo la temperatura exterior. Esto me da el siguiente resultado:

que es mejor, pero no es lo que esperaba ver.

Answer 1

Vas a tener un poco de problemas para modelar una serie con 2 niveles de estacionalidad utilizando un modelo ARIMA. Conseguirlo va a depender en gran medida de la configuración correcta de las cosas. ¿Ha considerado ya un modelo lineal simple? Son mucho más rápidos y fáciles de ajustar que los modelos ARIMA, y si se utilizan variables ficticias para los diferentes niveles de estacionalidad, suelen ser bastante precisos.

1. Asumo que tienes datos por hora, así que asegúrate de que tu objeto TS está configurado con una frecuencia de 24.
2. Puede modelar otros niveles de estacionalidad utilizando variables ficticias. Por ejemplo, puede querer un conjunto de variables ficticias 0/1 que representen el mes del año.
3. Incluya las variables ficticias en el xreg junto con cualquier covariable (como la temperatura).
4. Ajuste el modelo con la función arima en R base. Esta función puede manejar modelos ARMAX mediante el uso de la función xreg argumento.
5. Pruebe el Arima y auto.arima auto.arima es bueno porque encontrará automáticamente buenos parámetros para su modelo arima. Sin embargo, tardará una eternidad en ajustarse a su conjunto de datos.
6. Pruebe la función tslm del paquete arima, utilizando variables ficticias para cada nivel de estacionalidad. Esto se ajustará mucho más rápido que el modelo Arima, e incluso puede funcionar mejor en su situación.
7. Si 4/5/6 no funcionan, entonces empieza a preocuparte por las funciones de transferencia. Hay que arrastrarse antes de poder caminar.
8. Si desea realizar una previsión en el futuro, primero tendrá que prever sus variables xreg. Esto es fácil para los dummies estacionales, pero tendrá que pensar en cómo hacer una buena previsión meteorológica. ¿Quizás utilizar la mediana de los datos históricos?

He aquí un ejemplo de cómo lo enfocaría yo:

#Setup a fake time series
set.seed(1)
library(lubridate)
index <- ISOdatetime(2010,1,1,0,0,0)+1:8759*60*60
month <- month(index)
hour <- hour(index)
usage <- 1000+10*rnorm(length(index))-25*(month-6)^2-(hour-12)^2
usage <- ts(usage,frequency=24)

#Create monthly dummies.  Add other xvars to this matrix
xreg <- model.matrix(~as.factor(month))[,2:12]
colnames(xreg) <- c('Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct','Nov','Dec')

#Fit a model
library(forecast)
model <- Arima(usage, order=c(0,0,0), seasonal=list(order=c(1,0,0), period=24), xreg=xreg)
plot(usage)
lines(fitted(model),col=2)

#Benchmark against other models
model2 <- tslm(usage~as.factor(month)+as.factor(hour))
model3 <- tslm(usage~as.factor(month))
model4 <- rep(mean(usage),length(usage))

#Compare the 4 models
library(plyr) #for rbind.fill
ACC <- rbind.fill(  data.frame(t(accuracy(model))),
                    data.frame(t(accuracy(model2))),
                    data.frame(t(accuracy(model3))),
                    data.frame(t(accuracy(model4,usage)))
                )
ACC <- round(ACC,2)
ACC <- cbind(Type=c('Arima','LM1','Monthly Mean','Mean'),ACC)
ACC[order(ACC$MAE),]

Answer 2

He estado usando R para hacer previsiones de carga durante un tiempo y puedo sugerirte que uses forecast y sus valiosas funciones (como auto.arima ).

Puede construir un modelo ARIMA con el siguiente comando:

model = arima(y, order, xreg = exogenous_data)

con y su predicción (supongo que dayy ), order el orden de su modelo (considerando la estacionalidad) y exogenous_data su temperatura, la radiación solar, etc. La función auto.arima le ayuda a encontrar el orden óptimo del modelo. Puede encontrar un breve tutorial sobre el paquete `forecast'. aquí .

Answer 3

Personalmente no entiendo de funciones de transferencia, pero creo que tienes la xtransf y xreg invertido. Al menos en la base de R arima es xreg que contiene sus variables exógenas. Tengo la impresión de que una función de transferencia describe comment (los datos retardados afectan a los valores futuros) en lugar de qué .

Yo intentaría usar xreg para sus variables exógenas, quizás utilizando arima si arimax exige una función de transferencia. El problema es que su modelo es diario, pero sus datos tienen estacionalidad diaria y anual, y no estoy seguro ahora mismo de si una primera diferencia (la order=(*, 1, *) ) se encargará de ello o no. (Desde luego, no se obtendrán previsiones mágicas para todo el año con un modelo que sólo tenga en cuenta la estacionalidad diaria).

P.D. ¿Qué es el time que utiliza en su lm ? ¿Tiempo de reloj literal o un número de observación de 1 arriba? Creo que se podría obtener algo utilizando un modelo de efectos mixtos ( lmer en el lme4 ), aunque no he averiguado si al hacerlo se tiene en cuenta correctamente la autocorrelación que se producirá en una serie temporal. Si no se tiene en cuenta, que un lm no lo hace, es posible que obtenga un ajuste interesante, pero su concepto de la precisión de su predicción será demasiado optimista.