knitr::opts_chunk$set( collapse = TRUE, comment = "#>" )
In dieser kurzen Simulationsstudie sollen die Möglichkeiten des Paketes weiter veranschaulicht werden. Dazu wird zuerst wieder das Paket TimeSeries geladen.
library(TimeSeries) library(tibble) library(ggplot2)
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.
Das Modell hat die folgende Form, wobei $\epsilon$ "white noise" ist:
$$ X_t=\phi X_{t-1}+\epsilon_t $$ Die Varianz $\gamma_0$ berechnet sich folgendermaßen:
$$ \gamma(0)=Var(X_t)=Var(\phi X_{t-1}+\epsilon_t)=Var(\phi X_{t-1})+\sigma_{\epsilon}^2=\phi^2 Var(X_{t-1})+\sigma_{\epsilon}^2= \frac{\sigma_{\epsilon}^2}{1-\phi^2} $$ Nun sind auch die weitern Autokovarianzen leicht bestimmbar:
$$ \gamma(1)=Cov(X_{t+1},X_t)=Cov(\phi X_t+\epsilon_t,X_t)=Cov(\phi X_t,X_t)+Cov(\epsilon_{t+1},X_t)=\phi Cov(X_t,X_t)+0=\phi \gamma(0) $$ $$ \gamma(2)=Cov(X_{t+2},X_t)=...=\phi^2 \gamma_0 $$ $$ usw... $$ Nun wird die theoretische ACF (ersten 50 Autokovarianzen) für einen AR(1)-Prozess mit $\phi=0.8$ berechnet. Später wird die Standardnormalverteilung für das Rauschen verwendet. Somit gilt $\sigma_{\epsilon}^2=1$.
gamma_0 <- 1/(1-0.8^2) true_acf <- numeric(50) for(i in 0:49){ gamma_i <- gamma_0*0.8^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("Theoretical ACF: AR(1)") true_acf_plot
Es ist die typische langsam abflachende Autokovarianzfunktion eines AR(1) zu erkennen.
Nun werden Daten mit den oben beschriebenen Eigenschaften simuliert und die ACFs geschätzt. Es werden Zeitreihen unterschiedlicher Länge erzeugt um zu sehen, ob sich die geschätzten ACFs mit größerer Stichprobengröße dem theoretischen ACF annähern. Es werden lange "Burn-in"-Perioden gewählt, damit die Zeitreihe kaum von den zufälligen Startwerten abhängt. Dazu wird eine Funktion erstellt, welche die simulierten ACFs zusammen mit dem theoretischen ACF plottet.
set.seed(12345) ar_obj_50 <- arma_sim(phi = 0.8, n=50, burnin=1000) ar_obj_100 <- arma_sim(phi = 0.8, n=100, burnin=1000) ar_obj_200 <- arma_sim(phi = 0.8, n=200, burnin=1000) ar_obj_500 <- arma_sim(phi = 0.8, n=500, burnin=1000) ar_obj_1000 <- arma_sim(phi = 0.8, n=1000, burnin=1000) ar_obj_5000 <- arma_sim(phi = 0.8, n=5000, burnin=1000) ar <- arma_sim(phi = 0.8, n=50, burnin=1000) acf_plot <- function(ar, ...){ sim_tbl <- tibble(ACF=acf(ar)[1:50],Lag=1:50) sim_tbl$Group <- "Simulation" true_tbl <- tibble(ACF=true_acf,Lag=1:50) true_tbl$Group <- "Theoretical" tbl <- rbind(sim_tbl, true_tbl) acf_plot <- ggplot(tbl, aes(x=factor(Lag), y=ACF, fill=factor(Group))) + geom_bar(stat='identity', position='dodge') + ggtitle(...) + scale_x_discrete(name="Lag", breaks=seq(0,50,5)) acf_plot } acf_plot(ar_obj_50, "50 obs") acf_plot(ar_obj_100, "100 obs") acf_plot(ar_obj_200, "200 obs") acf_plot(ar_obj_500, "500 obs") acf_plot(ar_obj_1000, "1000 obs") acf_plot(ar_obj_5000, "5000 obs")
Der obige Vergleich zeigt deutlich, dass sich die geschätzte ACF dem theoretischen ACF des AR(1)-Prozesses bei höherer Beobachtungszahl annähert. Bei niedriger Beobachtungszahl werden die Autokovarianzen für die niedrigen Lags tendenziell unterschätzt und oszillieren für die höhren Lags deutlich um die null. Für die hohen Beobachtungszahlen ab 1000 können nurnoch kleine Unterschiede festgestellt werden.
In diesem Abschnitt soll untersucht werden ob der Durbin-Levinson Algorithmus gute Vorhersagen für simulierte Daten liefert. Dafür werden zum einen 250 Datenpunkte verschiedener ARMA-Prozesse simuliert. Nach 50 Beobachtungen wird eine "expanding window out-of-sample" Vorhersage durchgeführt. Diese 50 Beobachtungen werden übersprungen, da der Algorithmus eine bestimmte Menge an Daten benötigt um zuverlässig zu funktionieren. Also werden alle Daten bis zum Zeitpunkt t verwendet um die Realisation in t+1 vorherzusagen. Danach wird der Datensatz für die nächste Vorhersage um die Realisation in t+1 erweitert. Für diese Art von Forecast wird eine Funktion definiert.
dl_forecast <- function(x){ x <- x$arma len <- length(x) forecast <- numeric(len) for(t in 0:(len-1)){ if(t<=49){ forecast[t+1] <- NA next} x_temp <- x[1:t] dl <- DL(x_temp) forecast[t+1] <- sum(dl*x_temp) } forecast }
Im folgenden werden beispielhaft Daten zu folgenden Prozessen simuliert: AR(1), AR(2), AR(3), MA(1), Ma(2), ARMA(1,1).
set.seed(123455) dl_sim1 <- arma_sim(phi = 0.9, n=250, burnin=1000) dl_sim2 <- arma_sim(phi = c(0.6,0.2), n=250, burnin=1000) dl_sim3 <- arma_sim(phi = c(0.6,0.5,-0.2), n=250, burnin=1000) dl_sim4 <- arma_sim(theta = 0.5, n=250, burnin=1000) dl_sim5 <- arma_sim(theta = c(0.6,-0.2), n=250, burnin=1000) dl_sim6 <- arma_sim(theta = c(0.6,-0.2,0.7), n=250, burnin=1000) dl_sim7 <- arma_sim(phi=0.5, theta = 0.5, n=250, burnin=1000)
Für die AR-Prozesses ist es leicht zu testen ob der Durbin-Levinson Algorithmus erfolgreich war, da dieser die AR-Coeffizienten schätzt. Beispielsweise bei einem AR(1)-Prozess sollte der Durbin-Levinson für den ersten Coeffizienten den richtigen AR-Coeffizienten schätzen und für die restlichen Coeffizienten Schätzungen gegen Null liefern. Dies wird im folgenden untersucht:
# AR(1) - (0.9) - maximales p tail(DL(dl_sim1),10) # AR(2) - (0.6,0.2) - maximales p tail(DL(dl_sim2),10) # AR(3) - (0.6,-0.1,0.3) - maximales p tail(DL(dl_sim3),10)
Die Coeffizienten stimmen meist ungefähr mit den simulierten Prozessen überein, sodass der Algorithmus später möglicherweise auch sinnvolle Ergbenisse bei der Vorhersage erzielen wird.
Nun werden für alle diese Prozesse die Vorhersagen berechnet.
dl_1 <- dl_forecast(dl_sim1) dl_2 <- dl_forecast(dl_sim2) dl_3 <- dl_forecast(dl_sim3) dl_4 <- dl_forecast(dl_sim4) dl_5 <- dl_forecast(dl_sim5) dl_6 <- dl_forecast(dl_sim6) dl_7 <- dl_forecast(dl_sim7)
Nun sollen die Ergbnisse visualisiert werden. Dazu wird die simulierte Zeitreihe mit den Vorhersagen jeweils in einem Plot dargestellt.
plot_dl_forecasts <- function(series, forecast, ...){ tbl <- tibble(Series=series$arma, Forecast=forecast, Lag=1:250) ggplot(tbl, aes(x=Lag)) + geom_line(aes(y = Series), color = "darkred", size=0.4) + geom_line(aes(y = Forecast), color="steelblue", size=0.4) + ggtitle(...) } plot_dl_forecasts(dl_sim1, dl_1, "AR(1) - (0.9)") plot_dl_forecasts(dl_sim2, dl_2, "AR(2) - (0.6,0.2)") plot_dl_forecasts(dl_sim3, dl_3, "AR(3) - (0.6,-0.1,0.3)") plot_dl_forecasts(dl_sim4, dl_4, "Ma(1) - (0.5)") plot_dl_forecasts(dl_sim5, dl_5, "Ma(2) - (0.6,-0.2)") plot_dl_forecasts(dl_sim6, dl_6, "MA(3) - (0.6,-0.2,0.7)") plot_dl_forecasts(dl_sim7, dl_7, "ARMA(1,1) - (0.5, 0.5)")
Wie erwartet zeigen sich für vor allem für die AR-Prozesse passable Ergebnisse. Trotz recht wenigen Datenpunkten folgen die one-step-ahead Forecats der Zeitreihe zumindest in der Tendenz. Es zeigen sich vor allem bessere Ergebnisse für spätere Zeitpunkte, für die mehr Daten zur Verfügung stehen. Für MA-Prozesse ist theoretisch der Innovations-Alogrithmus besser geeignet. Für den ARMA(1,1)-Prozess sind die Ergebnisse des Durbin-Levinson Algortihmus auch recht zufriedenstellend. Allgemein kann man sagen umso mehr "Struktur" in den Zeitreihen zu finden ist umso besser ist es auch möglich Vorhersagen mit dem Durbin-Levinson Algorithmus zu erzeugen.
Was bei diesem Algorithmus noch auffällt, ist das die meisten geschätzen Parameter theoretisch gegen Null gehen sollten. Wegen dem Rauschen und der niedrigen Beobachtungszahl bei der Berechnung der Autokovarianzen aus hohen Lags ist dies in der Praxis nicht der Fall. Somit könnten möglicherweise die Vorhersagen verbessert werden, wenn man annimmt, dass die letzten 10 Beobachtungen besonders einflussreich für die Vorhersage sein sollten. Dafür werden die Rekursionen beim Durbin-Levinson auf 10 begrenzt, sodass dieser 10 $\phi's$ ausgibt. Diese werden dann für das forecasten verwendet. Dies wird in folgender Funktion umgesetzt:
dl_forecast_10 <- function(x){ x <- x$arma len <- length(x) forecast <- numeric(len) for(t in 0:(len-1)){ if(t<=49){ forecast[t+1] <- NA next} x_temp <- x[1:t] dl <- DL(x_temp, p=10L) forecast[t+1] <- sum(dl*tail(x_temp,10)) } forecast }
Nun werden mit dieser neuen Funktion Forecasts erzeugt.
dl_1 <- dl_forecast_10(dl_sim1) dl_2 <- dl_forecast_10(dl_sim2) dl_3 <- dl_forecast_10(dl_sim3) dl_4 <- dl_forecast_10(dl_sim4) dl_5 <- dl_forecast_10(dl_sim5) dl_6 <- dl_forecast_10(dl_sim6) dl_7 <- dl_forecast_10(dl_sim7)
Danach werden die selben Plots wie oben nochmals mit den neuen Vorhersagen erzeugt.
plot_dl_forecasts(dl_sim1, dl_1, "AR(1) - (0.9) - p=10") plot_dl_forecasts(dl_sim2, dl_2, "AR(2) - (0.6,0.2) - p=10") plot_dl_forecasts(dl_sim3, dl_3, "AR(3) - (0.6,-0.1,0.3) - p=10") plot_dl_forecasts(dl_sim4, dl_4, "Ma(1) - (0.5) - p=10") plot_dl_forecasts(dl_sim5, dl_5, "Ma(2) - (0.6,-0.2) - p=10") plot_dl_forecasts(dl_sim6, dl_6, "MA(3) - (0.6,-0.2,0.7) - p=10") plot_dl_forecasts(dl_sim7, dl_7, "ARMA(1,1) - (0.5, 0.5) - p=10")
Die Plots zeigen, dass durch die Beschränkung der verwendeten Autokovarianzen deutlich bessere Ergebnisse erzielt werden können. Dem zugrunde liegt natürlich eine Annahme über die Ordnung des Prozesses. Meine Interpretation ist das sich viele Schätzungenauigkeiten akkumulieren, wenn so viele $\phi's$ wie Beobachtungen geschätzt werden. Auch die Schätzungen der letzten Autokovarianzen im ACF ist ungenau, da hierfür wenige (im Fall der letzten Autokovarianz nur eine) Beobachtungen vorliegen. Diese Ungenauigkeit schlägt sich dann auch die Vorhersagen durch den Durbin-Levinson Algortihmus durch. Somit dürfte die Wahl der optimalen Anzahl von Rekursionen im Alogrithmus von der Einschätzung der Länge des Prozesses und Anzahl der Datenpunkte abhängen.
In diesem Abschnitt soll untersucht werden, ob der Innovation Algorithmus gute Vorhersagen für simulierte Daten liefert. Es wir dafür dasselbe Verfahren und die selben Funktionen wie beim DL Algorithmus angewendet, damit die beiden vergleichbar sind. Wir benutzen hier für die innovation_prediction Funktion.
theta_1 <- innovation(dl_sim1) theta_4<- innovation(dl_sim4) # Ar(1)-Prozess print(theta_1$theta[5:15,1:5]) # Ma(1)-Prozess print(theta_4$theta[5:15,1:5])
Doch zuerst schauen wir uns die Thetas des MA(1) und des Ar(1) Prozesses des Innovation Algorithmus an. Beim Ma(1) erkennt man, dass die erste Spalte nahe an unseren theta = 0,5 ist und die restlichen Werte Betrags kleiner als 0,1 sind. In der Theorie gibt der Innovation Algorithmus sogar den MA(1)-Prozess wieder. Dies bekommt der Algorithmus nicht hin, da die unsere Werte noch ein Rauschen haben und oder es ist der Maschinenungenauigkeit verschuldet. Zumindest ist $\hat{X}{n+1}\approx\theta{n,1}(X_{n}-\hat{X}_{n})$, was nahe am Prozess liegt und somit soll der Algorithmus ein guter Schätzer dafür sein.
Beim AR(1) sind mindestens die fünf Spalte nicht Betrags kleiner als 0,1, weshalb dort $\hat{X}_{n+1}$ von mehr Variabeln stärker abhängt, was das Vorhesagen ungenauer machen sollte.
Nun überprüfen durch die Simulation wie oben, ob unsere Werte mit der Theorie übereinstimmen.
in_sim1 <- in_sim2 <- in_sim3 <- in_sim4 <- in_sim5 <- in_sim6 <- in_sim7 <- rep(NA,50) for(i in 49:(length(dl_sim1$arma)-1)) { in_sim1[i+1] <- innovation_prediction(dl_sim1, lag.max = i) in_sim2[i+1] <- innovation_prediction(dl_sim2, lag.max = i) in_sim4[i+1] <- innovation_prediction(dl_sim4, lag.max = i) in_sim5[i+1] <- innovation_prediction(dl_sim5, lag.max = i) in_sim7[i+1] <- innovation_prediction(dl_sim7, lag.max = i) }
plot_in_forecasts <- function(series, forecast, ...){ tbl <- tibble(Series=series$arma, Forecast=forecast, Lag=1:250) ggplot(tbl, aes(x=Lag)) + geom_line(aes(y = Series), color = "darkred", size=0.4) + geom_line(aes(y = Forecast), color="darkgreen", size=0.4) + ggtitle(...) } plot_in_forecasts(dl_sim1, in_sim1, "AR(1) - (0.9)") plot_in_forecasts(dl_sim2, in_sim2, "AR(2) - (0.6,0.2)") plot_in_forecasts(dl_sim4, in_sim4, "Ma(1) - (0.5) - p=10") plot_in_forecasts(dl_sim5, in_sim5, "Ma(2) - (0.6,-0.2)") plot_in_forecasts(dl_sim7, in_sim7, "ARMA(1,1) - (0.5, 0.5)")
Allgemein kann man erkennen, dass der Algorithmus nach 100 Lags deutlich verbessert, aber bei den Spitzen meist noch zu niedrig schätzt. Das kann daran liegen, dass er durch mehr Werte genauer wird und die Extreme sind meist die unwahrscheinlichsten Ereignisse, was sie schwieriger zum Schätzen macht. Die Verbesserung zwischen MA und AR Werte ist nicht zu sehen. Dadurch dass der AR-Prozess deutlich weniger Spitzen hat, kann der Innovation Algorithmus sogar den Prozess besser vorhersagen. Auch der ARMA(1,1) Prozess lässt sich mit derselben Begründung gut durch den Algorithmus vorhersagen. Schlussendlich vergleichen wir noch beide Algorithmen.
plot_comp_forecasts <- function(series, forecast1, forecast2, ...){ tbl <- tibble(Series=series$arma[200:250], Forecast1=forecast1[200:250], Forecast2=forecast2[200:250], Lag=200:250) ggplot(tbl, aes(x=Lag)) + geom_line(aes(y = Series), color = "darkred", size=0.7) + geom_line(aes(y = Forecast1), color="steelblue", size=0.4) + geom_line(aes(y = Forecast2), color="darkgreen", size=0.4) + ggtitle(...) } plot_comp_forecasts(dl_sim1,dl_1,in_sim1, "AR(1) - (0.9)") plot_comp_forecasts(dl_sim4,dl_4,in_sim4, "MA(1) - (0.5)")
Man sieht in der näheren Analyse, dass beide Algorithmen nur die Tendenzen von den jeweiligen Graphen hinbekommen. Dies zeigt, wie schwer schon allein bei so einfachen Prozessen die Vorhersage von Werten ist. Man kann erkennen, dass der Innovation Algorithmus besser beim Ma Prozess ist und DL Algorithmus beim AR Prozess ist, da er besser die Extreme vorhersagt und der DL Algorithmus sich nur in der Mitte bewegt. Aber wie in der Analyse zum Innovation Algorithmus angemerkt wurde, ist dieser ein erstaunlich guter Schätzer für die AR Prozesse und somit lassen sich nicht wirklich herrausdeuten, ob der DL Algorithmus der besser Schätzer ist. Dies kann an einfach an unseren simulierten Werten liegen.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.