Comment adapter un modèle ARIMAX avec R?

33

J'ai quatre séries chronologiques différentes de mesures horaires:

  1. La consommation de chaleur à l'intérieur d'une maison
  2. La température à l'extérieur de la maison
  3. Le rayonnement solaire
  4. La vitesse du vent

Je veux pouvoir prédire la consommation de chaleur à l'intérieur de la maison. Il y a une nette tendance saisonnière, à la fois sur une base annuelle et sur une base quotidienne. Puisqu'il existe une corrélation claire entre les différentes séries, je souhaite les ajuster à l'aide d'un modèle ARIMAX. Cela peut être fait dans R, en utilisant la fonction arimax du paquet TSA.

J'ai essayé de lire la documentation sur cette fonction et de lire sur les fonctions de transfert, mais jusqu'à présent, mon code:

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)

Donne moi: entrez la description de l'image ici

où la ligne noire représente les données mesurées réelles et la ligne verte représente mon modèle ajusté en comparaison. Non seulement ce n'est pas un bon modèle, mais il est clair que quelque chose ne va pas.

Je reconnais que ma connaissance des modèles ARIMAX et des fonctions de transfert est limitée. Dans la fonction arimax (), (si j'ai bien compris), xtransf est la série temporelle exogène que je souhaite utiliser (à l'aide de fonctions de transfert) pour prédire ma série temporelle principale. Mais quelle est la différence entre xreg et xtransf?

Plus généralement, qu'est-ce que j'ai mal fait? Je voudrais pouvoir obtenir un meilleur ajustement que celui obtenu à partir de lm (chaleur ~ temp radi vent * temps).

Modifications: sur la base de certains commentaires, j'ai supprimé le transfert et ajouté xreg à la place:

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)

où dayy est le "numéro jour de l'année" et time, l'heure du jour. La température est à nouveau la température extérieure. Cela me donne le résultat suivant:

entrez la description de l'image ici

ce qui est meilleur, mais pas à la hauteur de ce que je pensais voir.

inutile
la source

Réponses:

34

Vous aurez un peu de difficulté à modéliser une série à 2 niveaux de saisonnalité à l'aide d'un modèle ARIMA. Obtenir ce droit est très dépendant de l’organisation correcte des choses. Avez-vous déjà envisagé un modèle linéaire simple? Ils sont beaucoup plus rapides et faciles à ajuster que les modèles ARIMA, et si vous utilisez des variables nominales pour vos différents niveaux de saisonnalité, elles sont souvent assez précises.

  1. Je suppose que vous avez des données horaires, alors assurez-vous que votre objet TS est configuré avec une fréquence de 24.
  2. Vous pouvez modéliser d'autres niveaux de saisonnalité à l'aide de variables nominales. Par exemple, vous voudrez peut-être un ensemble de maquettes 0/1 représentant le mois de l'année.
  3. Incluez les variables nominales dans l' xregargument, avec toutes les covariables (comme la température).
  4. Ajustez le modèle avec la fonction arima en base R. Cette fonction peut gérer les modèles ARMAX grâce à l'utilisation de l' xregargument.
  5. Essayez les fonctions Arima et auto.arima dans le paquet de prévisions. auto.arima est bien car il trouvera automatiquement les bons paramètres pour votre modèle arima. Cependant, il vous faudra FOREVER pour vous adapter à votre jeu de données.
  6. Essayez la fonction tslm dans le package arima en utilisant des variables nominales pour chaque niveau de saisonnalité. Cela ira beaucoup plus vite que le modèle Arima et pourrait même mieux fonctionner dans votre situation.
  7. Si 4/5/6 ne fonctionne pas, ALORS commencez à vous soucier des fonctions de transfert. Vous devez ramper avant de pouvoir marcher.
  8. Si vous envisagez de prévoir dans le futur, vous devez d’abord prévoir vos variables xreg. C'est facile pour les mannequins saisonniers, mais vous devrez réfléchir à la manière de faire de bonnes prévisions météorologiques. Peut-être utiliser la médiane des données historiques?

Voici un exemple de la façon dont j'aborderais ceci:

#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),]
Zach
la source
Quelle est la fonction équipée (). Si je l'utilise, j'obtiens des résultats bien meilleurs que ceux de Predict (model10, newxreg = regParams).
Utdiscant
@utdiscant: predict()est utilisé pour les prévisions, alors que fitted()le modèle est ajusté sur la période historique. Si vous souhaitez une aide plus spécifique, vous devez poster un exemple reproductible avec du code.
Zach
@utdiscant: de plus, si vous utilisez dayy comme xreg, vous courez le risque de sur-adapter car vous n'avez que 24 observations par jour. Vous obtiendrez peut-être de meilleurs résultats si vous utilisez le mois de l'année.
Zach
@utdiscant: De plus, vos xreg basés sur le temps doivent être des variables factices . Selon votre modélisation actuelle, vous vous attendez heatà augmenter linéairement avec l'heure du jour, puis à reculer lorsque l'heure revient à 1. Si vous utilisez des variables nominales, chaque heure de la journée produira son propre effet. Parcourez mon exemple de code et faites très attention à la façon dont je construis mon objet xreg.
Zach
Un inconvénient des fonctions ARIMA dans les packages statset forecastest qu’elles ne conviennent pas aux fonctions de transfert de prober. La documentation de la stats::arimafonction indique ce qui suit: Si un terme xreg est inclus, une régression linéaire (avec un terme constant si include.mean est true et en l'absence de différenciation) est équipée d'un modèle ARMA pour le terme d'erreur. Donc, si vous avez réellement besoin d’adapter des fonctions de transfert, il semble que cette TSA::arimaxfonction soit la solution R.
Christoffer
8

J'utilise R depuis un certain temps pour faire des prévisions de charge et je peux vous suggérer d'utiliser forecastpackage et ses fonctions inestimables (comme auto.arima).

Vous pouvez construire un modèle ARIMA avec la commande suivante:

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

avec yvotre prévision (je suppose dayy), orderl'ordre de votre modèle (en tenant compte de la saisonnalité) et exogenous_datavotre température, le rayonnement solaire, etc. La fonction auto.arimavous aide à trouver l'ordre optimal du modèle. Vous pouvez trouver un bref tutoriel sur le paquet `previsions ' ici .

Matteo De Felice
la source
Ce qui doit être prédit est la chaleur (la consommation de chaleur de la maison).
Utdiscant
3

Personnellement, je ne comprends pas les fonctions de transfert, mais je pense que vous avez compris xtransfet xreginversé. Au moins dans la base de R, arimac'est xregqu'il contient vos variables exogènes. J'ai l'impression qu'une fonction de transfert décrit comment (les données décalées affectent les valeurs futures) plutôt que quoi .

Je voudrais essayer d'utiliser xregpour vos variables exogènes, peut-être en utilisant arimasi arimaxdemande une fonction de transfert. Le problème est que votre modèle est quotidien, mais que vos données ont une saisonnalité journalière et annuelle, et je ne suis pas sûr pour le moment si une première différence (le order=(*, 1, *)) va prendre en charge cela ou pas. (Vous ne obtiendrez certainement pas de prévisions magiques tout au long de l'année avec un modèle qui prend uniquement en compte la saisonnalité quotidienne.)

PS Quelle est la timeque vous utilisez dans votre lm? Heure de l'horloge ou un numéro d'observation 1-up? Je pense que vous pourriez obtenir quelque chose en utilisant un modèle à effets mixtes ( lmerdans le lme4package), bien que je n’aie pas compris si cela rend compte correctement de l’autocorrélation qui se produira dans une série chronologique. Si ce n'est pas pris en compte, ce qui lmn'est pas le cas, vous obtiendrez peut-être un ajustement intéressant, mais votre idée de la précision de vos prédictions sera beaucoup trop optimiste.

Wayne
la source
J'ai à la fois l'heure de la mesure et le "jour de l'année" de la mesure.
Utdiscant