Como ajustar um modelo ARIMAX ao R?

Eu tenho quatro séries temporais diferentes de medições horárias:

1. O consumo de calor dentro de uma casa
2. A temperatura fora de casa
3. A radiação solar
4. A velocidade do vento

Quero poder prever o consumo de calor dentro de casa. Existe uma clara tendência sazonal, tanto anualmente como diariamente. Como existe uma correlação clara entre as diferentes séries, desejo ajustá-las usando um modelo ARIMAX. Isso pode ser feito em R, usando a função arimax do pacote TSA.

Tentei ler a documentação desta função e ler as funções de transferência, mas até agora meu 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 dá:

onde a linha preta é os dados reais medidos e a linha verde é meu modelo ajustado em comparação. Não é apenas um bom modelo, mas claramente algo está errado.

Admito que meu conhecimento dos modelos ARIMAX e das funções de transferência é limitado. Na função arimax (), (tanto quanto eu entendi), xtransf é a série temporal exógena que eu quero usar (usando funções de transferência) para prever minha principal série temporal. Mas qual é realmente a diferença entre xreg e xtransf?

De maneira mais geral, o que fiz de errado? Gostaria de conseguir um ajuste melhor do que o obtido com lm ​​(heat ~ temp radi wind * time).

Edições: com base em alguns comentários, removi a transferência e adicionei o xreg:

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)

onde dayy é o "número do dia do ano" e time é a hora do dia. Temp é novamente a temperatura externa. Isso me dá o seguinte resultado:

o que é melhor, mas não é o que eu esperava ver.

Você terá problemas para modelar uma série com 2 níveis de sazonalidade usando um modelo ARIMA. Conseguir isso certo depende muito da configuração correta das coisas. Você já considerou um modelo linear simples? Eles são muito mais rápidos e fáceis de ajustar do que os modelos ARIMA, e se você usar variáveis ​​fictícias para seus diferentes níveis de sazonalidade, elas geralmente são bastante precisas.

1. Suponho que você tenha dados por hora, portanto, verifique se o seu objeto TS está configurado com uma frequência de 24.
2. Você pode modelar outros níveis de sazonalidade usando variáveis ​​fictícias. Por exemplo, você pode querer um conjunto de 0/1 manequins representando o mês do ano.
3. Inclua as variáveis ​​fictícias no xregargumento, juntamente com quaisquer covariáveis ​​(como temperatura).
4. Ajuste o modelo com a função arima na base R. Essa função pode manipular modelos ARMAX através do uso do xregargumento.
5. Experimente as funções Arima e auto.arima no pacote de previsão. O auto.arima é bom porque encontrará automaticamente bons parâmetros para o seu modelo arima. No entanto, o FOREVER levará para caber no seu conjunto de dados.
6. Experimente a função tslm no pacote arima, usando variáveis ​​fictícias para cada nível de sazonalidade. Isso se ajustará muito mais rápido que o modelo Arima e pode até funcionar melhor em sua situação.
7. Se 4/5/6 não funcionar, ENTÃO comece a se preocupar com as funções de transferência. Você precisa rastejar antes de poder andar.
8. Se você planeja prever no futuro, primeiro precisará prever suas variáveis ​​xreg. Isso é fácil para manequins sazonais, mas você terá que pensar em como fazer uma boa previsão do tempo. Talvez use a mediana dos dados históricos?

Aqui está um exemplo de como eu abordaria isso:

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

Estou usando o R para fazer previsão de carga há um tempo e posso sugerir que você use o forecastpacote e suas funções inestimáveis ​​(como auto.arima).

Você pode construir um modelo ARIMA com o seguinte comando:

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

com ysua previsão (suponho dayy), ordera ordem do seu modelo (considerando a sazonalidade) e exogenous_datasua temperatura, radiação solar etc. A função auto.arimaajuda a encontrar a ordem ideal do modelo. Você pode encontrar um breve tutorial sobre o pacote `forecast ' aqui .

Pessoalmente, não entendo as funções de transferência, mas acho que você entendeu xtransfe xregreverteu. Pelo menos na base de R arimaé xregque contém suas variáveis ​​exógenas. Tenho a impressão de que uma função de transferência descreve como (dados atrasados ​​afetam valores futuros) e não o que .

Eu tentaria usar xregpara suas variáveis ​​exógenas, talvez usando arimase arimaxexige uma função de transferência. O problema é que seu modelo é diário, mas seus dados têm sazonalidade diária e anual, e não tenho certeza no momento se uma primeira diferença (o order=(*, 1, *)) resolverá isso ou não. (Você certamente não obterá previsões mágicas durante todo o ano de um modelo que considera apenas a sazonalidade diária.)

PS O timeque você usa no seu lm? Hora literal do relógio ou número de observação de 1 em 1? Eu acho que você pode obter algo usando um modelo de efeito misto ( lmerno lme4pacote), embora eu não tenha descoberto se fazer isso corretamente é responsável pela autocorrelação que ocorrerá em uma série temporal. Se não for considerado, o que lmnão é verdade, você pode ter um ajuste interessante, mas seu conceito de quão precisa é sua previsão será otimista demais.