In leopopeo/Time-Series: Analysen von Zeitreihen

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

library(TimeSeries)
library(tibble)
library(ggplot2)
library(itsmr)
library(dplyr)

Im Folgenden geben wir einige Beispiele und Anwendungen der Funktionen des Pakets Time Series an. Im ersten Teil werden wir dazu eine Simulationsstudie durchführen und im zweiten Teil reale Daten analysieren:

1. Simulationsstudie

Autokovarianz-Funktion
Durbin-Levinson Algorithmus
Innovation Prediction

2. Analyse realer Daten

Temperaturdaten aus Mannheim

Simulationsstudie

In dieser Simulationsstudio werden die Fähigkeiten des Paketes durch kleine Veränderungen der Funktions-Parameter veranschaulicht.

Autokovarianz-Funktion

In diesem Teil soll die Funktion acf getestet werden. Dazu wird ein AR(1)-Prozess mithilfe der Funktion arma_sim simuliert. Für diesen Prozess ist es leicht möglich die "wahre", theoretische Autokovarianz-Funktion zu berechnen.

Theoretische ACF eines AR(1)-Prozesses

Wir nehmen an, dass die Zeitreihe der folgenden Vorraussetzung genügt, wobei $\epsilon$ "white noise" ist:

$$ X_t=\phi X_{t-1}+\epsilon_t $$ Für die Varianz $\gamma_0$ können wir folgende Berechnungen anstellen:

$$ \gamma(0)=Cov(X_t,X_t)=Cov(\phi X_{t-1}+\epsilon,\phi X_{t-1} + \epsilon_t)=\phi^2 \gamma(0) +\sigma^2= \frac{\sigma^2}{1-\phi^2} $$ Nun sind auch die weiteren Autokovarianzen leicht bestimmbar:

$$ \gamma(h)=Cov(X_{t},X_{t-h})=Cov(\phi X_{t-1},X_{t-h})+Cov(\epsilon_t,X_{t-h})=\phi\gamma(h-1)+0=... = \phi^h \gamma(0) $$ Es gilt also:

$$ \gamma(1)=Cov(X_{t},X_{t-1})=...=\phi \gamma_0 $$ $$ \gamma(2)=Cov(X_{t},X_{t-2})=...=\phi^2 \gamma_0 $$ $$ .\ .\ . $$ (Brockwell et.al. (2002), S. 17f). Wir berechnen nun die theoretischen, ersten 50 Autokovarianzen für den AR(1)-Prozess mit $\phi=0.9$. Wir verwenden eine Standardnormalverteilung für das Rauschen, daher gilt $\sigma_{\epsilon}^2=1$.

gamma_0 <- 1/(1-0.9^2)
true_acf <- numeric(50)
for(i in 0:49){
        gamma_i <- gamma_0*0.9^i
        true_acf[i+1] <- gamma_i
}
true_tbl <- tibble(ACF=true_acf, Lag=seq_along(true_acf))
true_acf_plot <- ggplot(true_tbl, aes(x=Lag, y=ACF)) +
                geom_bar(stat="identity") + 
                ggtitle("Theoretische ACF: AR(1)")
true_acf_plot

Es ist die typische langsam abflachende Autokovarianzfunktion eines AR(1) zu erkennen.

Vergleich von theoretischer mit geschätzter ACF

Wir simulieren nun Daten mit den oben beschriebenen Eigenschaften und schätzen die ACF-Werte ab. Dafür testen wir, wie sich die Zeitreihe mit unterschiedlichen Längen verhält. Dabei wollen wir untersuchen, ab welcher Stichprobengröße sich der geschätzte ACF-Wert dem theoretischen ACF annähert und betrachten dabei jeweils die ersten 50 Werte. Wir erstellen dafür eine Funktion, welche die simulierten ACFs zusammen mit den theoretischen plottet.

set.seed(13)
ar_50 <- arma_sim(phi=0.9,I=50, sd =1)
ar_100 <- arma_sim(phi=0.9,I=100, sd = 1)
ar_500 <- arma_sim(phi=0.9,I=500, sd =1)
ar_2000 <- arma_sim(phi=0.9,I=2000, sd =1)
ar_5000 <- arma_sim(phi=0.9,I=5000, sd =1)

acf_plot <- function(objekt, ...){
        sim_tbl <- tibble(ACF=ACF(objekt)[1:50],Lag=1:50)
        sim_tbl$Group <- "Simulation"

        true_tbl <- tibble(ACF=true_acf,Lag=1:50)
        true_tbl$Group <- "Theoretisch"

        tbl <- rbind(sim_tbl, true_tbl)

        acf_plot <- ggplot(tbl, aes(x=factor(Lag), y=ACF, fill=Group)) +
                geom_bar(stat='identity', position='dodge') + 
                ggtitle(...) +
                scale_x_discrete(name="Lag", breaks=seq(0,50,5))
acf_plot
}
acf_plot(ar_50, "Länge 50") 
acf_plot(ar_100, "Länge 100")
acf_plot(ar_500, "Länge 500")
acf_plot(ar_2000, "Länge 2000")
acf_plot(ar_5000, "Länge 5000")

Als zusätzliche Veranschaulichung zeigt der folgende Plot den Zusammenhang zwischen der Länge der Zeitreihe und der Differenz zwischen den theoretischen und simulierten Werten an.

set.seed(13)
acf_diff_plot <- function(n){
        res <- tibble(Laenge=n,Differenz=numeric(n))
        zaehler <- 1
        for(i in 1:n){
          a <- ACF(arma_sim(phi=0.9,I=i*50, sd =1))
          diff <- sum(abs(a[1:50]-true_acf))
          res$Laenge[zaehler] <- i*50
          res$Differenz[zaehler] <- diff
          zaehler <- zaehler+1
        }

        acf_plot <- ggplot(res,aes(x=Laenge,y=Differenz)) +
          geom_bar(stat='identity', position='dodge') +
          ggtitle("Differenz der simulierten und theoretischen Daten")

        acf_plot
}
######################unebdingt reinmachen nur in Klammer wegen langer berechnungszeit!!!!!!!
#acf_diff_plot(100)

Der obige Vergleich zeigt deutlich, dass sich die geschätzte ACF der theoretischen ACF des AR(1)-Prozesses bei größerer Länge annähert. Die Übereinstimmung nimmt bis zu einer Länge der AR(1)-Reihe von 1000 sichtbar zu. Ab diesem Zeitpunkt oszilliert die Differenz erstaunlicherweise auf relativ hohem Niveau. Wir nehmen an, dass diese Oszillation durch den zugrundeliegenden Zufallsprozess erzeugt wird.

Durbin-Levinson Algorithmus

In diesem Abschnitt befassen wir uns mit dem Durbin-Levinson Algorithmus und wollen uns seiner Vorhersagen bezüglich der ARMA-Zeitreihen anschauen. Außerdem analyisieren wir, welchen Einfluss die Parameter $\phi$ und $\theta$ haben und für welche Zeitreihen der Durbin-Levinson Algorithmus sich besonders gut eignet. Hierfür erstellen wir eine Prediciton-Funktion, welche uns für die vorherigen $n$ Werte den $n+1$ Wert vorhersagt:

DL_Helper <- function(X){

  prediction <- rep(NA, length(X))

  for(i in 10:(length(X))){
    x_help <- X[1:i]
    dl_save <- DLA(x_help)
    prediction[i] <- sum(dl_save*rev(x_help))
  }
  prediction
}

Nun simulieren wir Daten mit verschiedenen Eigenschaften und beobachten, welche Auswirkungen sich bei der Änderung der $\phi$ und $\theta$ bei gleichbleibender Stichprobengröße ergeben und für welche Zeitreihen der DLA akkurat ist. Hierfür betrachten wir AR(1), AR(2), MA(1), MA(2) und ARMA(1)-Zeitreihen:

set.seed(42)

AR_1 <- arma_sim(phi = 0.6, sd = 0.1, I = 100)
AR_2 <- arma_sim(phi = c(0.5, 0.2), sd = 0.1, I = 100)
MA_1 <- arma_sim(theta = 0.4, sd = 0.1, I = 100)
MA_2 <- arma_sim(theta = c(0.65, 0.25), sd = 0.1, I = 100)
ARMA_1 <- arma_sim(phi = 0.6, theta = 0.4, sd = 0.1, I = 100)

Im folgenden berechnen wir die Vorhersagen der gegebenen ARMA-Zeitreihen anhand des DL Algorithmus:

DL1 <- DL_Helper(AR_1)
DL2 <- DL_Helper(AR_2)
DL3 <- DL_Helper(MA_1)
DL4 <- DL_Helper(MA_2)
DL5 <- DL_Helper(ARMA_1)

Wir betrachten im folgenden immer Ein-Schritt-Prognosen. Zur Beurteilung der Vorhersagegüte ziehen wir den $MSE$ heran, welcher definiert ist als

$$ MSE(X,\hat X) := \frac{1}{n}\sum_{i = 1}^{n} (X_i - \hat X_i)^2. $$

plot_mult <- function(series, forecast1, forecast2 = NULL,  title = "Zeitreihe"){
  tbl <- tibble(Series=series, forecast1 = forecast1, forecast2 = forecast2, Lag=1:length(series))

    a <- ggplot(tbl, aes(x=Lag)) 
    b <- geom_line(aes(y = Series, color = "Series"), size = 0.5, na.rm = T) 
    c <- geom_line(aes(y = forecast1, color = "forecast1"), size = 0.5, na.rm = T)
    if(is.null(forecast2)){d <- NULL }
    else{d <- geom_line(aes(y = forecast2, color = "forecast2"), size = 0.5, na.rm = T)} 
    e <- ggtitle(title)
    plt <- a+b+c+d+e
    plt
}

MSE <- function(X_true, X_pred){
  n = length(X_true)
  return(1/n*sum((X_true-X_pred)^2))
}

MSE(AR_1[-(1:9)], DL1[-(1:9)])
plot_mult(series = AR_1, forecast1 =  DL1,  title = "AR(1)")
MSE(AR_2[-(1:9)], DL2[-(1:9)])
plot_mult(series = AR_2, forecast1 = DL2, title = "AR(2)")
MSE(MA_1[-(1:9)], DL3[-(1:9)])
plot_mult(series = MA_1, forecast1 = DL3, title = "MA(1)")
MSE(MA_2[-(1:9)], DL4[-(1:9)])
plot_mult(series = MA_2, forecast1 = DL4, title = "MA(2)")
MSE(ARMA_1[-(1:9)], DL5[-(1:9)])
plot_mult(series = ARMA_1, forecast1 = DL5, title = "ARMA(1)")

Die Plots verdeutlichen, dass der Durbin-Levinson Algorithmus bereits mit wenigen Inputs sehr gute Vorhersagen trifft. Der $MSE = 4.88e-3$ ist beim $AR(2)$ Prozess am größten, bei den anderen liegt er sehr nahe beieinander. Die beste Vorhersage nach $MSE$ liefert der Durbin-Levinson-Algorithmus beim $MA(2)$-Prozess mit $MSE = 1.78e-3$. Auffallend ist, dass der Durbin-Levinson Algorithmus entgegen unseren Erwartungen für die $MA$-Prozesse besser ist als für die $AR$-Prozesse. Außerdem verringert sich die Ungenauigkeit des Algorithmus, je mehr Input dieser bekommt. Die ist schön zu erkennen, da bereits bei 60+ die Vorhersage fast exakt mit der Zeitreihe übereinstimmt.

Jetzt wollen wir noch überprüfen, welche Ergebnisse der Algorithmus liefert, wenn er mit seinen Vorhersagen weiterrechnet.

DL_Helper_Flaw <- function(X){

  prediction <- rep(NA, length(X))

  for(i in 10:(length(X))){
    x_help <- X[1:i]
    dl_save <- DLA(x_help)
    prediction[i] <- sum(dl_save*rev(x_help))
    X[i] <- prediction[i]
  }
  prediction
}

DL_Flaw1 <- DL_Helper_Flaw(AR_1)
DL_Flaw2 <- DL_Helper_Flaw(MA_1)
DL_Flaw3 <- DL_Helper_Flaw(ARMA_1)

MSE(AR_1[-(1:9)], DL_Flaw1[-(1:9)])
plot_mult(AR_1, forecast1 = DL_Flaw1, title = "AR(1)")
MSE(MA_1[-(1:9)], DL_Flaw2[-(1:9)])
plot_mult(MA_1, forecast1 =DL_Flaw2, title ="MA(1)")
MSE(ARMA_1[-(1:9)], DL_Flaw3[-(1:9)])
plot_mult(ARMA_1, forecast1 = DL_Flaw3,title ="ARMA(1,1)")

Auch hier ist der $MSE$ beim $MA(1)$-Prozess mit $2.00e-3$ am geringsten. Obwohl er optisch die Peaks nicht so gut erfasst, wie beim $AR(1)$-Prozess. Die Vorhersage ist insbesondere immer flacher und nicht so spitz wie in unserer tatsächlichen Zeitreihe. Das hängt unter anderem damit zusammen, dass die ACF Funktion selbst starke Peaks erzeugt und der DLA nicht mit dem Rauschen klarkommt. Folglich schlägt der DLA nicht so stark aus, wie der Innovation Algorithmus. Der deutlich größere $MSE$ beim $ARMA(1,1)$-Prozess wird hauptsächlich durch den großen Vorhersagefehler bei ca. Lag = 40 verursacht.

Innovation Algorithmus

Wir betrachten zunächst die letzten Zeilen und ersten Spalten der Innovationmatrix für den AR(1) und den MA(2) - Prozess.

M = innovation(AR_1, lag = 99)[90:99,1:6]
rownames(M) = 90:99; M

M = innovation(MA_2, lag = 99)[90:99,1:6]
rownames(M) = 90:99; M

Beim $MA(2)$-Prozess entsprechen die ersten beiden Spalten ungefähr den spezifizierten $MA$-Parametern der Zeitreihe. Der Algorithmus sollte also besonders gut für die Vorhersage von MA-Prozessen geeignet sein.

Innovation Prediction

Dieser Abschnitt veranschaulicht die Vorhersagegenauigkeit des Innovation Algorithmus. Um einen Vergleich mit dem Durbin-Levinson Algortihmus zu gewährleisten, nutzen wir die analoge Methodik, sowie die selben Zeitreihen, wie oben:

Hierfür erstellen wir eine Prediciton-Funktion, welche uns für die vorherigen $n$ Werte den $n+1.$ Wert vorhersagt:

innovation_prediction <- function(X){
   inno_prediction <- rep(NA, length(X))

   for(i in 10:(length(X)-1)){
     inno_help <- X[1:i]
     inno_prediction[i+1] <- ts_predict(inno_help,1)
   }
   inno_prediction
}

innovation_sim_AR_1 <- innovation_prediction(AR_1)
innovation_sim_AR_2 <- innovation_prediction(AR_2)
innovation_sim_MA_1 <- innovation_prediction(MA_1)
innovation_sim_MA_2 <- innovation_prediction(MA_2)
innovation_sim_ARMA_1 <- innovation_prediction(ARMA_1)

MSE(AR_1[-(1:10)], innovation_sim_AR_1[-(1:10)])
plot_mult(AR_1, innovation_sim_AR_1, title = "AR(1)")

MSE(AR_2[-(1:10)], innovation_sim_AR_2[-(1:10)])
plot_mult(AR_2, innovation_sim_AR_2, title = "AR(2)")

MSE(MA_1[-(1:10)], innovation_sim_MA_1[-(1:10)])
plot_mult(MA_1, innovation_sim_MA_1, title ="MA(1)")

MSE(MA_2[-(1:10)], innovation_sim_MA_2[-(1:10)])
plot_mult(MA_2, innovation_sim_MA_2, title = "MA(2)")

MSE(ARMA_1[-(1:10)], innovation_sim_ARMA_1[-(1:10)])
plot_mult(ARMA_1, innovation_sim_ARMA_1, title ="ARMA(1,1)")

Zunächst können wir festhalten, dass die Ein-Schritt-Prognose des Innovation-Algorithmus für alle Prozesse einen größeren $MSE$ liefert als mit Durbin-Levinson. Bei den $MA$-Prozessen ist der $MSE$ geringer als bei den Prozessen mit $AR$-Anteil. Dies entspricht auch nach Brockwell unseren Erwartungen. Die Ausschläge der Extremwerte sind bei der Vorhersage ebenfalls nicht so stark, wie bei der tatsächlichen Zeitreihe. Dies liegt daran, dass zur Prognose ein gewichtetes Mittel der vergangenen Differenzen $\mathbf{X_n} - \mathbf{\hat X_n}$ herangezogen wird, was dazu führt, dass die Prognose tendenziell immer zum Mittelwert strebt. Dieses Verhalten der Vorhersage des Innovation Algorithmus, Bewegungen zu antizipieren, aber zu dämpfen, sieht man besonders gut an folgendem Beispiel eines vollständig deterministischen, periodischen Prozesses:

X = sin(1:65)
forc = ts_predict(X, steps = 30)

plot_timeseries(X,forc,"Sinus-Funktion mit Vorhersage")

Zur besseren Übersichtlichkeit stellen wir hier für die$AR(1)$-, $MA(1)$- und $ARMA(1,1)$-Prozesse die Vorhersage des Durbin-Levinson- und des Innovation-Algorithmus sowie die $MSE$s gegenüber.

MSE(AR_1[-(1:10)], innovation_sim_AR_1[-(1:10)])
MSE(AR_1[-(1:10)], DL1[-(1:10)])
plot_mult(AR_1, innovation_sim_AR_1,DL1, title = "AR(1)")

MSE(MA_1[-(1:10)], innovation_sim_MA_1[-(1:10)])
MSE(MA_1[-(1:10)], DL3[-(1:10)])
plot_mult(MA_1, innovation_sim_MA_1,DL3, title = "MA(1)")

MSE(ARMA_1[-(1:10)], innovation_sim_ARMA_1[-(1:10)])
MSE(ARMA_1[-(1:10)], DL5[-(1:10)])
plot_mult(ARMA_1, innovation_sim_ARMA_1,DL5, title = "ARMA(1)")

Nun wollen wir unseren Innovation-Vorhersage mit der bestehenden Funktion "forecast" aus dem Paket "itsmr" (https://mparison <- georgeweigt.github.io/itsmr-refman.pdf, Seite 14) vergleichen: Zunächst betrachten wir das Ergebnis einer Ein-Schritt-Prognose der forecast-Funktion mit unserer ts_predict-Funktion eines $MA(1)$-Prozesses. Der kleine $MSE$ von $9.73e-3$ ist bei der forecast-Prognose nur geringfügig kleiner als bei unserer Funktion ts_predict ($MSE = 10.8e-3$).

#Funktion zur Ein-Schritt-Prognose mit forecast
forecast_prediction <- function(X, p = 0, q = 0){
   inno_prediction <- rep(NA, length(X))

   for(i in 10:(length(X)-1)){
     inno_help <- X[1:i]
     inno_prediction[i+1] <- forecast(inno_help,
                                      M = NULL,
                                      a = arma(inno_help, p = p, q = q),
                                      h = 1,
                                      opt = 0,
                                      alpha = 1)$pred
   }
   inno_prediction
}

itsmr_onestep_AR <- forecast_prediction(AR_1, p = 1, q = 0)
itsmr_onestep_MA <- forecast_prediction(MA_1, p = 0, q = 1)
itsmr_onestep_ARMA <- forecast_prediction(MA_1, p = 1, q = 1)

MSE(AR_1[-(1:11)], innovation_sim_AR_1[-(1:11)])
MSE(AR_1[-(1:11)], itsmr_onestep_AR[-(1:11)])
plot_mult(AR_1, innovation_sim_AR_1, itsmr_onestep_AR, "AR(1). ts_predict red, forecast green.")

MSE(MA_1[-(1:11)], innovation_sim_MA_1[-(1:11)])
MSE(MA_1[-(1:11)], itsmr_onestep_MA[-(1:11)])
plot_mult(MA_1, innovation_sim_MA_1, itsmr_onestep_MA, "MA(1). ts_predict red, forecast green.")

MSE(ARMA_1[-(1:11)], innovation_sim_ARMA_1[-(1:11)])
MSE(ARMA_1[-(1:11)], itsmr_onestep_ARMA[-(1:11)])
plot_mult(ARMA_1, innovation_sim_ARMA_1, itsmr_onestep_ARMA, "ARMA(1). ts_predict red, forecast green.")

Der $MSE$ der Ein-Schritt-Prognose ist bei allen drei Prozessen für forecast und ts_predict nahe beieinander. Für den ARMA-Prozess ist unser Algorithmus - gemessem am $MSE$ - besser.

Betrachten wir nun im folgenden noch eine Mehr-Schritt Prognose.

itsmr_prediction_1 <- forecast(AR_2, #Time series data
                             M = NULL, #Data model
                             a = arma(AR_2, p=2, q=0), #ARMA model
                             h = 10, #Steps ahead
                             opt = 0, #Display option (0 silent, 1 tabulate, 2 plot and tabulate)
                             alpha = 1)$pred #Level of significance
innovation_sim_1 <- ts_predict(AR_2, 10)
itsmr_prediction_2 <- forecast(MA_2, #Time series data
                             M = NULL, #Data model
                             a = arma(MA_2, p=2, q=0), #ARMA model
                             h = 10, #Steps ahead
                             opt = 0, #Display option (0 silent, 1 tabulate, 2 plot and tabulate)
                             alpha = 1)$pred #Level of significance
innovation_sim_2 <- ts_predict(MA_2, 10)

plot_timeseries(AR_2, itsmr_prediction_1, "AR_2 mit forecast")
plot_timeseries(AR_2, innovation_sim_1, "AR_2 mit ts_predict" )
plot_timeseries(MA_2, itsmr_prediction_2, "MA_2 mit forecast")
plot_timeseries(MA_2, innovation_sim_2, "MA_2 mit ts_predict")

Für eine Mehr-Schritt-Prognose haben wir leider nicht das erwartete Ergebnis erhalten. In der Dokumentation der forecast-Funktion steht zwar, dass die Funktion ebenfalls den gleichen Innovation Algorithmus als Grundlage für die Berechnung benutzt. Es ist jedoch keine gute Übereinstimmung zwischen unserer Vorhersage und der forecast-Vorhersage zu erkennen. Auch für eine Änderung der Werte p und q in forecast ergibt sich kein besseres Ergebnis.

Reale Daten

Da Temperaturdaten eine starke saisonale Komponente haben, eignen sich diese für eine Analyse mit einem ARMA-Modell. Wir analysieren im folgenden Temperaturdaten des Deutschen Wetterdienstes (DWD) für Mannheim: auf stündlicher sowie auf monatlicher Basis.

data_hourly <- as_tibble(read.delim("Luftdaten_stunde_19480101_20201231_05906.txt",
                                    header = T,
                                    sep = ";",
                                    dec = "."))

data_mon <- as_tibble(read.delim("Luftdaten_monat_18810101_20201231_05906.txt",
                                 header = T,
                                 sep = ";",
                                 dec = "."))

data_hourly %>%
  select(time = MESS_DATUM, temp = TT_TU) %>%
  filter((abs(temp) < 50) & (2000010100 <= time) ) -> data_hourly
data_mon %>%
  select(date = MESS_DATUM_BEGINN, temp = MO_TT) %>%
  filter((abs(temp) < 50) & (20100101 <= date) ) %>%
  mutate(date = as.Date(as.character(date), "%Y%m%d"))-> data_mon


plot(data_mon, type = "l",main = "Monatliche Daten")

periodo <- perio((data_mon$temp)-mean(data_mon$temp))
periodo[67] = 0.1 #f(0) = 7e-30 -> stark verzerrter Plot
plot_periodogram(periodo, logscale = T)

Hier sehen wir ein Periodogramm für die um den Mittelwert bereinigten Monatsdaten. Der Peak der jährlichen Periode ist eindeutig zu erkennen (beachte die Logarithmische Skala).

Als nächstes betrachten wir einen kleinen Ausschnitt der Daten und versuchen anhand von 100 Datenwerten eine Vorhersage zu treffen.

plot_predic <- function(series, dl, ...){
  len <- length(series)
  dl[1:100] <- rep(NA,100)
  tbl <- tibble(Series=series, DL_obj = dl, Lag=1:len)

   ggplot(tbl, aes(x=Lag)) + 
    geom_line(aes(y = Series), color = "blue", size=0.4,na.rm=T) + 
    geom_line(aes(y = DL_obj), color= "red", size=0.4,na.rm=T) +
    ggtitle(...)
}


raw_1 <- data_hourly$temp[3700:3800]-mean(data_hourly$temp[3700:3800])
raw_full_1 <- data_hourly$temp[3700:3830]-mean(data_hourly$temp[3700:3830])
predic_DL_1 <- DL_prediction(raw_1,len=30,all=T)
predic_IN_1 <- ts_predict(raw_1,steps=30,all=T)

raw_2 <- data_mon$temp[1:100]-mean(data_mon$temp[1:100])
raw_full_2 <- data_mon$temp[1:130]-mean(data_mon$temp[1:130])
predic_DL_2 <- DL_prediction(raw_2,len=30,all=T)
predic_IN_2 <- ts_predict(raw_2,steps=30,all=T)


acf_plot_1 <- ggplot(tibble(ACF=ACF(raw_1), Lag=seq_along(ACF(raw_1))), aes(x=Lag, y=ACF)) +
                geom_bar(stat="identity") + 
                ggtitle("ACF-Plot")
acf_plot_2 <- ggplot(tibble(ACF=ACF(raw_2), Lag=seq_along(ACF(raw_2))), aes(x=Lag, y=ACF)) +
                geom_bar(stat="identity") + 
                ggtitle("ACF-Plot")

Im Folgenden sehen wir zunächst für die stündlichen Werte neben der Zeitreihe (in blau) je eine Vorhersage für t=101 bis t=131 durch den Durbin-Levinson Algorithmus und dem Innovation-Algorithmus (in rot).

MSE(raw_full_1[-(1:100)], predic_DL_1[-(1:100)])
plot_predic(raw_full_1, predic_DL_1,"Vorhersage Durbin-Levinson")
MSE(raw_full_1[-(1:100)], predic_IN_1[-(1:100)])
plot_predic(raw_full_1, predic_IN_1,"Vorhersage Innovation")

Während die Vorhersage des Durbin-Levinson-Algorithmus stark von den realisierten Werten abweicht, ist die Vorhersage des Innovation-Algorithmus vor allem in den ersten Schritten sehr gut.

Betrachten wir nun die selbe Vorhersagen für die Monatsdaten.

MSE(raw_full_2[-(1:100)], predic_DL_2[-(1:100)])
plot_predic(raw_full_2, predic_DL_2,"Vorhersage Durbin-Levinson")
MSE(raw_full_2[-(1:100)], predic_IN_2[-(1:100)])
plot_predic(raw_full_2, predic_IN_2,"Vorhersage Innovation")

Hier läuft die Vorhersage des Durbin-Levinson schnell gegen den Mittelwert, ohne die Schwankungen zu antizipieren. Im Gegensatz dazu erkennt der Innovation-Algorithmus das periodische Verhalten der Temperaturschwankungen deutlich besser.

In beiden Beispielen liefert der Innovation-Algorithmus eine bessere Vorhersage, als der Durbin-Levinson-Algorithmus. Dies legt die Vermutung nahe, dass es sich hier eher um einen MA-Prozess handelt. Die entsprechenden ACF-Plots bestätigen dies:

acf_plot_1
acf_plot_2

Man erkennt das starke Fallen der ACF mit dem Lag, was typisch für einen MA-Prozess ist. Wie bereits oben erklärt können die Maxima nicht sehr gut vorhergesagt werden. Es lohnt sich also den entsprechenden ACF-Plot anzuschauen um zu entscheiden welche Prediction-Funktion die Vorhersage besser schätzen kann.