Zeitreihenvorhersage mit Random Forest

Benjamin Franklin sagte, dass nur zwei Dinge im Leben sicher sind: der Tod und die Steuern. Das erklärt, warum meine Kolleg:innen bei statworx nicht gerade begeistert waren, als sie mir vor ein paar Wochen von ihren Wochenendplänen erzählten: die Einkommenssteuererklärung machen. Mann, dachte ich, das ist echt mies, ich würde diese Zeit viel lieber draußen verbringen. Und dann war die Idee geboren.

Was könnten Steuern und die Natur wohl gemeinsam haben? Nun ja, ich fragte mich: können wir Steuereinnahmen mit Random Forest vorhersagen? (wild kreativ, ich weiß). Wenn wir uns mit Steuereinnahmen beschäftigen, betreten wir das Reich der Zeitreihen, beherrscht von fantastischen Wesen wie ARIMA, VAR, STLM und anderen. Das sind bewährte Methoden – warum also Random Forest?

Nun, du und ich sind uns wahrscheinlich einig, dass Random Forest einer der großartigsten Algorithmen ist: einfach, flexibel und leistungsstark. So sehr, dass Wyner et al. (2015) ihn als das „Off-the-shelf“-Tool für die meisten Data-Science-Anwendungen bezeichnen. Kurz gesagt: es ist einer dieser Algorithmen, die einfach funktionieren (wenn du genau wissen willst, wie – dann lies diesen großartigen Beitrag von meinem Kollegen Andre).

Random Forest ist ein Hammer, aber ist Zeitreihendaten ein Nagel?

Du hast Random Forest wahrscheinlich schon für Regression und Klassifikation verwendet – aber für Zeitreihenvorhersage? Halt, wirst du sagen, Zeitreihendaten sind etwas Besonderes! Und du hast recht. Wenn es um Daten mit einer zeitlichen Dimension geht, wird der Einsatz von Machine-Learning-(ML)-Methoden etwas knifflig.

Warum das? Nun, Random Forests – wie die meisten ML-Methoden – haben kein Zeitverständnis. Im Gegenteil: sie betrachten Beobachtungen als unabhängig und identisch verteilt. Diese Annahme ist bei Zeitreihendaten offensichtlich verletzt, denn diese zeichnen sich durch serielle Abhängigkeit aus.

Außerdem sind Random Forests bzw. baumbasierte Methoden nicht in der Lage, einen Trend vorherzusagen, d.h. sie extrapolieren nicht. Um zu verstehen, warum, erinnere dich daran, dass Bäume durch If-Then-Regeln funktionieren, die den Eingaberaum rekursiv aufteilen. Dadurch können sie keine Werte vorhersagen, die außerhalb des Wertebereichs des Zielwerts im Trainingsdatensatz liegen.

Also doch zurück zu ARIMA? Noch nicht! Mit ein paar Tricks können wir mit Random Forest Zeitreihen vorhersagen. Es braucht nur etwas Vor- und (Nach-)Bearbeitung. Dieser Blogpost zeigt dir, wie du Random Forests fürs Forecasting nutzen kannst!

Dabei sei gesagt, dass es verschiedene Wege gibt, das zu tun. So machen wir es: Wir bedienen uns aus dem Werkzeugkasten der Zeitreihenökonometrie an alten, aber bewährten Techniken – Differenzierung und statistischen Transformationen. Diese sind die Eckpfeiler der ARIMA-Modellierung – aber wer sagt, dass wir sie nicht auch für Random Forest verwenden dürfen?

Passend zum Thema nutzen wir eine Zeitreihe des Statistischen Bundesamtes zu den Lohn- und Einkommensteuereinnahmen in Deutschland von 1999–2018 (nach Steuerumverteilung). Du kannst die Daten hier herunterladen. Los geht’s!

Bereit machen für Machine Learning oder was steckt überhaupt in einer Zeitreihe?

Im Wesentlichen ist eine (univariate) Zeitreihe ein Vektor von Werten, die durch die Zeit indiziert sind. Damit sie „lernbar“ wird, müssen wir etwas Vorverarbeitung betreiben. Dazu gehören einige oder alle der folgenden Schritte:

• Statistische Transformationen (Box-Cox-Transformation, Log-Transformation etc.)
• Detrending (Differenzierung, STL, SEATS etc.)
• Time Delay Embedding (mehr dazu weiter unten)
• Feature Engineering (Lags, gleitende Statistiken, Fourier-Terme, Zeit-Dummies etc.)
• Der Kürze und Klarheit halber konzentrieren wir uns in diesem Beitrag auf die Schritte eins bis drei.

Okay, strukturieren wir das ein wenig: um Random Forest für Zeitreihendaten zu verwenden, machen wir TDE: transformieren, differenzieren und einbetten. Starten wir R und laden die benötigten Pakete sowie unsere Daten.

# load the packages
suppressPackageStartupMessages(require(tidyverse))
suppressPackageStartupMessages(require(tsibble))
suppressPackageStartupMessages(require(randomForest))
suppressPackageStartupMessages(require(forecast))

# specify the csv file (your path here)
file <- ".../tax.csv"

# read in the csv file
tax_tbl <- readr::read_delim(
  file = file,
  delim = ";",
  col_names = c("Year", "Type", month.abb),
  skip = 1,
  col_types = "iciiiiiiiiiiii",
  na = c("...")
) %>% 
  select(-Type) %>% 
  gather(Date, Value, -Year) %>% 
  unite("Date", c(Date, Year), sep = " ") %>% 
  mutate(
    Date = Date %>% 
      lubridate::parse_date_time("m y") %>% 
      yearmonth()
  ) %>% 
  drop_na() %>% 
  as_tsibble(index = "Date") %>% 
  filter(Date <= "2018-12-01")

# convert to ts format
tax_ts <- as.ts(tax_tbl)

‍

Bevor wir uns in die Analyse stürzen, sollten wir schnell nach impliziten und expliziten Fehlstellen in den Daten suchen. Das Paket tsibble verfügt über einige praktische Funktionen, um genau das zu tun:

# implicit missings
has_gaps(tax_tbl)

# explicit missings
colSums(is.na(tax_tbl[, "Value"]))

‍

Nein, sieht gut aus! Mit welcher Art von Zeitreihen haben wir es also zu tun?

# visualize
plot_org <- tax_tbl %>% 
  ggplot(aes(Date, Value / 1000)) + # to get the axis on a more manageable scale
  geom_line() +
  theme_minimal() +
  labs(title = "German Wage and Income Taxes 1999 - 2018", x = "Year", y = "Euros")

Differenzierung kann den Unterschied machen

Wenn du schon einmal mit klassischen Zeitreihenmodellen gearbeitet hast, bist du wahrscheinlich auf das Konzept der Differenzierung gestoßen. Der Grund dafür ist, dass klassische Zeitreihenmodelle stationäre Daten erfordern.

Stationarität bedeutet, dass der Mittelwert und die Varianz der Zeitreihe endlich sind und sich über die Zeit nicht verändern. Sie impliziert also eine gewisse Stabilität in den statistischen Eigenschaften der Zeitreihe. Wie wir im Plot sehen, ist unsere Zeitreihe weit davon entfernt! Es gibt einen Aufwärtstrend sowie ein ausgeprägtes saisonales Muster in der Reihe.

Wie hängen diese beiden Konzepte – Differenzierung und Stationarität – zusammen? Du weißt es wahrscheinlich schon oder hast es erraten: Differenzierung ist ein Weg, um nicht-stationäre Zeitreihendaten stationär zu machen. Das ist schön zu wissen, aber im Moment ist für uns vor allem wichtig, dass Differenzierung Veränderungen im Niveau einer Reihe und damit den Trend entfernt. Genau das, was wir für unseren Random Forest brauchen!

Wie macht man das? Hier nehmen wir einfach die ersten Differenzen der Daten, d.h. die Differenz zwischen aufeinanderfolgenden Beobachtungen $y_t' = y_t - y_{t-1}$. Das funktioniert auch mit einem saisonalen Lag $y_t' = y_t - y_{t-s}$, was bedeutet, dass die Differenz zwischen einer Beobachtung und einer früheren Beobachtung aus derselben Saison gebildet wird, z.B. November dieses Jahres und November des Vorjahres.

Während die Differenzierung den Mittelwert einer Zeitreihe stabilisieren kann, kann eine Box-Cox- oder Log-Transformation die Varianz stabilisieren. Die Familie der Box-Cox-Transformationen dreht sich um den Parameter Lambda:

$$\tilde{y}_t = \begin{cases} \log(y_t) & \lambda = 0\\ (y_t^\lambda - 1)/ \lambda & \lambda \neq 0 \end{cases}$$

Wenn Lambda gleich null ist, entspricht die Box-Cox-Transformation der Logarithmierung. Wir wählen diesen Wert, um die Rücktransformation unserer Prognosen einfach zu gestalten. Aber zögere nicht, mit verschiedenen Lambda-Werten zu experimentieren oder den „besten“ Wert mit Hilfe des forecast-Pakets zu schätzen.

# pretend we're in December 2017 and have to forecast the next twelve months
tax_ts_org <- window(tax_ts, end = c(2017, 12))

# estimate the required order of differencing
n_diffs <- nsdiffs(tax_ts_org)

# log transform and difference the data
tax_ts_trf <- tax_ts_org %>% 
  log() %>% 
  diff(n_diffs)

# check out the difference! (pun)
plot_trf <- tax_ts_trf %>% 
  autoplot() +
  xlab("Year") +
  ylab("Euros") +
  ggtitle("German Wage and Income Taxes 1999 - 2018") +
  theme_minimal()

gridExtra::grid.arrange(plot_org, plot_trf)

‍

Fassen wir zusammen, was wir bisher gemacht haben: Zuerst haben wir den Logarithmus der Daten genommen, um die Varianz zu stabilisieren. Dann haben wir die Daten einmal differenziert, um sie im Mittel stationär zu machen. Zusammen haben uns diese eher einfachen Transformationen von nicht-stationär zu stationär geführt.

Was kommt als Nächstes? Nun, jetzt nutzen wir die so transformierten Daten, um unseren Random Forest zu trainieren und Prognosen zu erstellen. Sobald wir die Vorhersagen erhalten, machen wir die Transformationen rückgängig, um sie wieder auf die ursprüngliche Skala der Daten zu bringen.

Nur noch ein Schritt, bevor wir zum Modellierungsteil kommen: Wir haben immer noch nur einen Vektor. Wie bringen wir diese Daten in eine Form, mit der ein ML-Algorithmus umgehen kann?

Enter the matrix: Time Delay Embedding

Um unserem Random Forest die transformierten Daten zuzuführen, müssen wir das, was im Wesentlichen ein Vektor ist, in eine Matrix verwandeln, also eine Struktur, mit der ein ML-Algorithmus arbeiten kann. Dafür nutzen wir ein Konzept namens Time Delay Embedding.

Time Delay Embedding stellt eine Zeitreihe in einem euklidischen Raum mit der Einbettungsdimension $K$ dar. Um dies in R zu tun, verwendet man die Basisfunktion embed(). Alles, was du tun musst, ist das Zeitreihenobjekt einzufügen und die Einbettungsdimension auf einen Wert größer als die gewünschte Anzahl an Lags zu setzen.

lag_order <- 6 # the desired number of lags (six months)
horizon <- 12 # the forecast horizon (twelve months)

tax_ts_mbd <- embed(tax_ts_trf, lag_order + 1) # embedding magic!

‍

Wenn du dir das Objekt tax_ts_mbd ansiehst, wirst du feststellen, dass du eine Matrix erhältst, in der die abhängige Variable in der ersten Spalte auf ihre Lags in den übrigen Spalten regressiert wird:

$$Y^K = \left[ \begin{array}{lllll} y_K \quad & y_{K-1} \quad & \dots \quad & y_2 \quad & y_1 \\ \vdots \quad & \vdots \quad & \vdots \quad & \vdots \quad & \vdots \\ y_i \quad & y_{i-1} \quad & \dots \quad & y_{i-K+2} \quad & y_{i-K+1} \\ \vdots \quad & \vdots \quad & \vdots \quad & \vdots \quad & \vdots \\ y_N \quad & y_{N-1} \quad & \dots \quad & y_{N-K+2} \quad & y_{N-K+1} \end{array} \right]$$

Time Delay Embedding ermöglicht es uns, jede lineare oder nichtlineare Regressionsmethode auf Zeitreihendaten anzuwenden, sei es Random Forest, Gradient Boosting, Support Vector Machines usw. Ich habe mich für ein Lag von sechs Monaten entschieden, aber du kannst auch mit anderen Lags experimentieren. Außerdem beträgt der Prognosehorizont zwölf, da wir die Steuereinnahmen für das Jahr 2018 vorhersagen.

Wenn es ums Forecasting geht, bin ich ziemlich direkt

In diesem Beitrag verwenden wir die sogenannte Direct Forecasting Strategy. Das bedeutet, dass wir $H$ separate Modelle $f_H$ schätzen – eines für jeden Prognosehorizont. Mit anderen Worten: Wir trainieren für jede zeitliche Distanz in den Daten ein separates Modell. Für ein großartiges Tutorial, wie das genau funktioniert, schau dir diesen Beitrag an.

Die Direct Forecasting Strategy ist weniger effizient als die Recursive Forecasting Strategy, bei der nur ein einziges Modell $f$ geschätzt und – wie der Name schon sagt – $H$-mal wiederverwendet wird. Rekursiv bedeutet in diesem Fall, dass jede Vorhersage als Input zurück ins Modell gegeben wird, um die nächste Vorhersage zu erzeugen.

Trotz dieses Nachteils hat die direkte Strategie zwei entscheidende Vorteile: Erstens leidet sie nicht unter einer Akkumulation von Prognosefehlern, und zweitens lässt sie sich leicht mit exogenen Prädiktoren kombinieren.

Wie man die Direct Forecasting Strategy implementiert, wird sehr schön im oben genannten Beitrag demonstriert, daher möchte ich das hier nicht wiederholen. Falls du wenig Zeit hast, hier die Kurzfassung: Wir verwenden die direkte Strategie, um Multi-Step-Ahead-Forecasts zu erzeugen. Das bedeutet, dass wir für jeden Prognosehorizont ein Modell trainieren, indem wir die Trainingsdaten schrittweise so umformen, dass sie die zeitliche Distanz zwischen den Beobachtungen abbilden.

y_train <- tax_ts_mbd[, 1] # the target
X_train <- tax_ts_mbd[, -1] # everything but the target

y_test <- window(tax_ts, start = c(2018, 1), end = c(2018, 12)) # the year 2018
X_test <- tax_ts_mbd[nrow(tax_ts_mbd), c(1:lag_order)] # the test set consisting
# of the six most recent values (we have six lags) of the training set. It's the
# same for all models.

‍

Wenn du bis hierhin mitgekommen bist – Respekt, wir sind fast am Ziel! Jetzt kommen wir zum spaßigen Teil: Wir lassen unseren Random Forest auf die Daten los. Wir trainieren das Modell in einer Schleife, wobei jede Iteration ein Modell für einen Prognosehorizont $h = 1, \dots, H$ erstellt.

Die Random-Forest-Prognose: Es sieht gut aus

Unten verwende ich den Random Forest direkt „out of the box“, ohne mir die Mühe zu machen, ihn zu tunen (ein Thema, dem ich gerne einen eigenen Beitrag widmen würde). Das mag faul wirken (und ist es vielleicht auch), aber ich habe den Prozess bewusst auf das Wesentliche reduziert, in der Hoffnung, möglichst klar zu zeigen, was hier passiert.

forecasts_rf <- numeric(horizon)

for (i in 1:horizon){
  # set seed
  set.seed(2019)

  # fit the model
  fit_rf <- randomForest(X_train, y_train)

  # predict using the test set
  forecasts_rf[i] <- predict(fit_rf, X_test)

  # here is where we repeatedly reshape the training data to reflect the time distance
  # corresponding to the current forecast horizon.
  y_train <- y_train[-1] 

  X_train <- X_train[-nrow(X_train), ] 
}

Alles klar, die Schleife ist durchgelaufen. Wir haben gerade zwölf Modelle trainiert und zwölf Vorhersagen erhalten. Da wir unsere Zeitreihe vor dem Training transformiert haben, müssen wir nun auch die Vorhersagen zurücktransformieren.

Zurück zum Ausgangspunkt oder wie wir Prognosen auf der ursprünglichen Skala erhalten

Da wir zuvor die Log-Transformation angewendet haben, ist die Rücktransformation recht unkompliziert. Wir kehren den Prozess von innen nach außen um, d.h., wir machen zuerst die Differenzierung rückgängig und dann die Log-Transformation. Das tun wir, indem wir die kumulierte Summe unserer transformierten Vorhersagen exponentieren und das Ergebnis mit der letzten Beobachtung unserer Zeitreihe multiplizieren. Mit anderen Worten, wir berechnen:

$$y_{t+H} = y_{t}\exp\left(\sum_{h = 1}^H \tilde{y}_{t+h}\right)$$

# calculate the exp term
exp_term <- exp(cumsum(forecasts_rf))

# extract the last observation from the time series (y_t)
last_observation <- as.vector(tail(tax_ts_org, 1))

# calculate the final predictions
backtransformed_forecasts <- last_observation * exp_term

# convert to ts format
y_pred <- ts(
  backtransformed_forecasts,
  start = c(2018, 1),
  frequency = 12
)

# add the forecasts to the original tibble
tax_tbl <- tax_tbl %>% 
  mutate(Forecast = c(rep(NA, length(tax_ts_org)), y_pred))

# visualize the forecasts
plot_fc <- tax_tbl %>% 
  ggplot(aes(x = Date)) +
  geom_line(aes(y = Value / 1000)) +
  geom_line(aes(y = Forecast / 1000), color = "blue") +
  theme_minimal() +
  labs(
    title = "Forecast of the German Wage and Income Tax for the Year 2018",
    x = "Year",
    y = "Euros"
  )

accuracy(y_pred, y_test)

‍

‍

Sieht so aus, als wäre unsere Prognose ziemlich gut! Wir haben einen MAPE von 2,6 Prozent erreicht. Aber da man sich nie allein auf Genauigkeitsmetriken verlassen sollte, berechnen wir schnell noch einen einfachen Benchmark wie das saisonale naive Modell. Das ist eine niedrige Hürde, also wenn unser Modell sie nicht überspringt, kommt es in die Tonne.

benchmark <- forecast(snaive(tax_ts_org), h = horizon)

tax_ts %>% 
  autoplot() +
  autolayer(benchmark, PI = FALSE)

accuracy(benchmark, y_test)

Die Fehlermetriken sind hier deutlich höher, also kann man mit Sicherheit sagen, dass unser Random Forest gute Arbeit geleistet hat.

Wie geht es jetzt weiter? Nun, wir haben es noch nicht ausprobiert, aber wir könnten unsere Vorhersagen mit etwas Hyperparameter-Tuning weiter verbessern. Und, unter uns gesagt: Vielleicht ist Random Forest ziemlich gut, aber nicht das beste Modell für diese Aufgabe. Wir könnten andere ausprobieren oder ein Ensemble von Modellen bilden.

Wenn du aus diesem Beitrag heute nur eine Sache mitnimmst, dann diese: Wir können effektive Zeitreihenvorhersagen mit Machine Learning machen, ohne gleich die großen Geschütze wie rekurrente neuronale Netze aufzufahren. Es braucht nur ein wenig Vor- und Nachbearbeitung. Warum also nicht beim nächsten Forecasting (oder beim Prokrastinieren deiner Steuererklärung) Random Forest in dein Repertoire aufnehmen?

UPDATE: Hi, einige Leute haben mich nach den Daten gefragt – ich weiß, die sind nicht ganz einfach zu finden (destatis ist nicht gerade benutzerfreundlich). Hier ist daher ein Link zum statworx GitHub, wo du die CSV-Datei herunterladen kannst. Nach Rückmeldungen haben wir dort auch eine ausführliche Erklärung des Codes für den Splitting-Prozess ergänzt.

Literatur

• Wyner, Abraham J., et al. “Explaining the success of adaboost and random forests as interpolating classifiers.” The Journal of Machine Learning Research 18.1 (2017): 1558-1590.