Wie rüste ich ein ARIMAX-Modell mit R aus?

33

Ich habe vier verschiedene Zeitreihen von Stundenmessungen:

  1. Der Wärmeverbrauch in einem Haus
  2. Die Temperatur außerhalb des Hauses
  3. Die Sonnenstrahlung
  4. Die Windgeschwindigkeit

Ich möchte den Wärmeverbrauch im Haus vorhersagen können. Es gibt einen klaren saisonalen Trend, sowohl auf jährlicher Basis als auch auf täglicher Basis. Da es eine eindeutige Korrelation zwischen den verschiedenen Serien gibt, möchte ich sie mit einem ARIMAX-Modell anpassen. Dies kann in R mit der Funktion arimax aus dem Paket TSA erfolgen.

Ich habe versucht, die Dokumentation zu dieser Funktion und die Übertragungsfunktionen zu lesen, aber bisher meinen 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)

gibt mir: Bildbeschreibung hier eingeben

wo die schwarze Linie die tatsächlichen gemessenen Daten sind und die grüne Linie mein angepasstes Modell im Vergleich ist. Es ist nicht nur kein gutes Modell, sondern offensichtlich stimmt etwas nicht.

Ich gebe zu, dass meine Kenntnisse über ARIMAX-Modelle und Übertragungsfunktionen begrenzt sind. In der Funktion arimax () ist xtransf (soweit ich verstanden habe) die exogene Zeitreihe, mit der ich (unter Verwendung von Übertragungsfunktionen) meine Hauptzeitreihe vorhersagen möchte. Aber was ist der Unterschied zwischen xreg und xtransf wirklich?

Was habe ich allgemein falsch gemacht? Ich möchte in der Lage sein, eine bessere Passform zu erreichen als die, die ich mit lm (wärmetemp. Radi wind * time) erreicht habe.

Änderungen: Aufgrund einiger Kommentare habe ich die Übertragung entfernt und stattdessen xreg hinzugefügt:

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)

Wobei dayy der "Zahlentag des Jahres" und time die Stunde des Tages ist. Temp ist wieder die Außentemperatur. Dies gibt mir das folgende Ergebnis:

Bildbeschreibung hier eingeben

Das ist besser, aber nicht annähernd das, was ich erwartet hatte.

utdiscant
quelle

Antworten:

34

Mit einem ARIMA-Modell können Sie eine Reihe mit zwei Saisonalitätsstufen leicht modellieren. Um dies zu erreichen, müssen die Dinge richtig eingerichtet werden. Haben Sie schon über ein einfaches lineares Modell nachgedacht? Sie lassen sich viel schneller und einfacher anpassen als ARIMA-Modelle. Wenn Sie Dummy-Variablen für Ihre unterschiedlichen Saisonalitätsstufen verwenden, sind sie häufig recht genau.

  1. Ich gehe davon aus, dass Sie stündliche Daten haben, stellen Sie also sicher, dass Ihr TS-Objekt mit einer Häufigkeit von 24 eingerichtet ist.
  2. Sie können andere Saisonalitätsstufen mithilfe von Dummy-Variablen modellieren. Sie möchten beispielsweise einen Satz von 0/1-Dummys, die den Monat des Jahres darstellen.
  3. Fügen Sie die Dummy-Variablen xregzusammen mit etwaigen Kovariaten (z. B. der Temperatur) in das Argument ein.
  4. Passen Sie das Modell mit der arima-Funktion in die Basis R ein. Diese Funktion kann ARMAX-Modelle mithilfe des xregArguments verarbeiten.
  5. Probieren Sie die Funktionen Arima und auto.arima im Vorhersagepaket aus. auto.arima ist nett, weil es automatisch gute Parameter für Ihr Arima-Modell findet. Es wird jedoch FÜR IMMER dauern, bis es auf Ihren Datensatz passt.
  6. Probieren Sie die Funktion tslm im Paket arima aus und verwenden Sie für jede Saisonalitätsstufe Dummy-Variablen. Dies passt viel schneller als das Arima-Modell und funktioniert möglicherweise sogar in Ihrer Situation besser.
  7. Wenn 4/5/6 nicht funktioniert, machen Sie sich DANN Gedanken über die Übertragungsfunktionen. Sie müssen kriechen, bevor Sie gehen können.
  8. Wenn Sie in die Zukunft prognostizieren möchten, müssen Sie zuerst Ihre xreg-Variablen prognostizieren. Dies ist leicht für Saison-Dummies, aber Sie müssen sich überlegen, wie Sie eine gute Wettervorhersage erstellen können. Vielleicht den Median der historischen Daten verwenden?

Hier ist ein Beispiel, wie ich das angehen würde:

#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
quelle
Was ist die gepaßte () Funktion. Wenn ich das benutze, erhalte ich weitaus bessere Ergebnisse als mit predict (model10, newxreg = regParams).
utdiscant
@utdiscant: predict()wird für Prognosen verwendet, während fitted()das Modell über den historischen Zeitraum angepasst wird . Wenn Sie spezifischere Hilfe benötigen, sollten Sie ein reproduzierbares Beispiel mit Code veröffentlichen.
Zach
@utdiscant: Wenn Sie dayy als xreg verwenden, besteht die Gefahr einer Überanpassung, da Sie nur 24 Beobachtungen pro Tag haben. Sie erhalten möglicherweise bessere Prognoseergebnisse, wenn Sie den Monat des Jahres verwenden.
Zach
@utdiscant: Außerdem müssen Ihre zeitbasierten Xregs Dummy-Variablen sein . Die Art und Weise, wie Sie es jetzt modellieren, ist, dass Sie erwarten heat, mit der Stunde des Tages linear zuzunehmen und dann wieder nach unten zu springen, wenn die Stunde auf 1 zurückkehrt. Wenn Sie Dummy-Variablen verwenden, wird jede Stunde des Tages ihren eigenen Effekt bekommen. Führen Sie meinen Beispielcode durch und achten Sie sorgfältig darauf, wie ich mein xreg-Objekt konstruiere.
Zach
Ein Nachteil der ARIMA-Funktionen in den Paketen statsund forecastist, dass sie nicht für Prober-Transfer-Funktionen geeignet sind. Die Dokumentation der stats::arimaFunktion gibt Folgendes an: Wenn ein xreg-Term enthalten ist, wird eine lineare Regression (mit einem konstanten Term, wenn include.mean wahr ist und keine Differenzierung vorliegt) mit einem ARMA-Modell für den Fehlerterm ausgestattet. Wenn Sie also tatsächlich Übertragungsfunktionen anpassen müssen, ist die Funktion anscheinend der TSA::arimaxrichtige Weg R.
Christoffer
8

Ich verwende R schon seit einiger Zeit, um Lastprognosen zu erstellen, und ich kann Ihnen vorschlagen, das forecastPaket und seine unschätzbaren Funktionen (wie auto.arima) zu verwenden.

Sie können ein ARIMA-Modell mit dem folgenden Befehl erstellen:

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

mit yIhrem Vorhersagewert (nehme ich an dayy), orderder Reihenfolge Ihres Modells (unter Berücksichtigung der Saisonalität) und exogenous_dataIhrer Temperatur, Sonneneinstrahlung usw. Die Funktion auto.arimahilft Ihnen, die optimale Modellreihenfolge zu finden. Ein kurzes Tutorial zum "Forecast" -Paket finden Sie hier .

Matteo De Felice
quelle
Was vorhergesagt werden muss, ist Wärme (der Wärmeverbrauch des Hauses).
utdiscant
3

Ich persönlich verstehe keine Übertragungsfunktionen, aber ich denke, Sie haben die xtransfund xregumgekehrt. Zumindest in R-Base arimaist es , xregdass Ihre exogenen Variablen enthält. Ich habe den Eindruck, dass eine Übertragungsfunktion eher beschreibt, wie (verzögerte Daten zukünftige Werte beeinflussen), als was .

Ich würde versuchen, xregfür Ihre exogenen Variablen zu verwenden, vielleicht mit, arimawenn arimaxeine Übertragungsfunktion verlangt. Das Problem ist, dass Ihr Modell täglich ist, Ihre Daten jedoch sowohl die tägliche als auch die jährliche Saison haben, und ich bin mir derzeit nicht sicher, ob ein erster Unterschied (die order=(*, 1, *)) dies erledigen wird oder nicht. (Mit einem Modell, das nur die tägliche Saisonabhängigkeit berücksichtigt, werden Sie auf keinen Fall magische Prognosen für das ganze Jahr erhalten.)

PS Was ist das time, was du in deinem verwendest lm? Wörtliche Uhrzeit oder eine 1-fache Beobachtungszahl? Ich denke, Sie könnten etwas mit einem Mischeffektmodell ( lmerim lme4Paket) erreichen, obwohl ich nicht herausgefunden habe, ob dies die Autokorrelation, die in einer Zeitreihe auftreten wird, korrekt berücksichtigt. Wenn dies nicht berücksichtigt wird lm, erhalten Sie möglicherweise eine interessante Übereinstimmung, aber Ihr Konzept, wie genau Ihre Vorhersage ist, ist viel zu optimistisch.

Wayne
quelle
Ich habe sowohl die Stunde der Messung als auch den "Tag des Jahres" der Messung.
16.