¿Cómo adaptar un modelo ARIMAX con R?

Tengo cuatro series temporales diferentes de mediciones por hora:

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. Existe una clara tendencia estacional, tanto anualmente como a diario. Dado que existe una clara correlación entre las diferentes series, quiero ajustarlas usando un modelo ARIMAX. Esto se puede hacer en R, utilizando la función arimax del paquete TSA.

Intenté leer la documentación sobre 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 medidos reales, y la línea verde es mi modelo ajustado en comparación. No solo no es un buen modelo, sino que claramente algo está mal.

Admitiré que mi conocimiento de los modelos ARIMAX y las funciones de transferencia es limitado. En la función arimax (), (por lo que he entendido), xtransf es la serie de tiempo exógena que quiero usar (usando funciones de transferencia) para predecir mi serie de tiempo principal. Pero, ¿cuál es la diferencia entre xreg y xtransf realmente?

En general, ¿qué he hecho mal? Me gustaría poder obtener un mejor ajuste que el que se obtiene de lm (heat ~ temp radi wind * time).

Ediciones: Basado en algunos de los comentarios, eliminé la transferencia y agregué 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 dayy es el "día número del año" y time es la hora del día. La temperatura vuelve a ser la temperatura exterior. Esto me da el siguiente resultado:

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

Tendrás algunos problemas para modelar una serie con 2 niveles de estacionalidad utilizando un modelo ARIMA. Hacer esto bien depende en gran medida de configurar las cosas correctamente. ¿Ya ha considerado un modelo lineal simple? Son mucho más rápidos y fáciles de ajustar que los modelos ARIMA, y si utiliza variables ficticias para sus diferentes niveles de estacionalidad, a menudo son bastante precisas.

1. Supongo que tiene datos por hora, así que asegúrese de que su objeto TS esté configurado con una frecuencia de 24.
2. Puede modelar otros niveles de estacionalidad utilizando variables ficticias. Por ejemplo, es posible que desee un conjunto de 0/1 dummies que representan el mes del año.
3. Incluya las variables ficticias en el xregargumento, junto con cualquier covariable (como la temperatura).
4. Ajuste el modelo con la función arima en la base R. Esta función puede manejar modelos ARMAX mediante el uso de xreg argumento.
5. Pruebe las funciones Arima y auto.arima en el paquete de pronóstico. auto.arima es bueno porque encontrará automáticamente buenos parámetros para su modelo arima. Sin embargo, llevará SIEMPRE encajar en su conjunto de datos.
6. Pruebe la función tslm en el 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 el 4/5/6 no funciona, ENTONCES comience a preocuparse por las funciones de transferencia. Tienes que gatear antes de poder caminar.
8. Si planea pronosticar en el futuro, primero deberá pronosticar sus variables xreg. Esto es fácil para los maniquíes estacionales, pero tendrá que pensar en cómo hacer un buen pronóstico del tiempo. ¿Quizás usar la mediana de los datos históricos?

Aquí hay un ejemplo de cómo abordaría esto:

#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),]

He estado usando R para hacer pronósticos de carga durante un tiempo y puedo sugerirle que use el forecastpaquete y sus funciones invaluables (como auto.arima).

Puede construir un modelo ARIMA con el siguiente comando:

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

con ysu predicción (supongo dayy), orderel orden de su modelo (considerando la estacionalidad) y exogenous_datasu temperatura, radiación solar, etc. La función lo auto.arimaayuda a encontrar el orden óptimo del modelo. Puede encontrar un breve tutorial sobre el paquete de "pronóstico" aquí .

Personalmente, no entiendo las funciones de transferencia, pero creo que obtuviste xtransfy xregrevertiste. Al menos en la base de R arimaes xregque contiene sus variables exógenas. Tengo la impresión de que una función de transferencia describe cómo (los datos rezagados afectan los valores futuros) en lugar de qué .

Intentaría usar xregpara sus variables exógenas, tal vez usando arimaif arimaxexige una función de transferencia. El problema es que su modelo es diario, pero sus datos tienen una estacionalidad diaria y anual, y no estoy seguro en este momento si una primera diferencia (la order=(*, 1, *)) se encargará de eso o no. (Ciertamente no obtendrá pronósticos mágicos durante todo el año de un modelo que solo considera la estacionalidad diaria).

PD: ¿Cuál es el timeque usas en tu lm? ¿Hora de reloj literal o un número de observación de 1 arriba? Creo que podría obtener algo utilizando un modelo de efectos mixtos ( lmeren el lme4paquete), aunque no he descubierto si hacerlo representa correctamente la autocorrelación que ocurrirá en una serie de tiempo. Si no se tiene en cuenta, lo que lmno es así, es posible que obtenga un ajuste interesante, pero su concepto de cuán precisa es su predicción será demasiado optimista.