In adrian1econ/TimeSeries: What the Package Does (One Line, Title Case)

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

Einleitung

Diese Vignette demonstriert die Fähigkeiten des Pakets TimeSeries, welches folgende Funktionen enthält:

Simulation eines ARMA(p,q)-Prozesses (arma_sim)
Berechnung der Autokovarianzfunktion einer Zeitreihe (acf)
Implementierung des Durbin-Levinson Algorithmus (DL)
Plotten von Zeitreihen (plot_ts)
Implementierung des Innovations Algorithmus (innovation) und Vorhersagen (innovation_prdeiction)
Implementierung des Periodogramms (Periodogram) und des zugehörigen Plots (plot_periodogram)

Dafür wird als erstes das Paket TimeSeries geladen.

library(TimeSeries)
library(ggplot2)

Simulation eines ARMA-Prozesses

Mithilfe der Funktion arma_sim können ARMA-Prozesse simuliert werden. ARMA-Prozesse sind allgemein eine Mischung aus autoregressiven (AR)- und moving average (MA)-Prozessen. Die so simulierten Zeitreihen sind auch die Grundlage der Simulationsstudie, in der die Funktionen des Pakets weiter getestet werden. Die Funktion arma_sim erzeugt eine ARMA-Zeitreihe, welche folgende Formel aus Brookwell et al. (2002) erfüllt: $$X_t-\phi_1X_{t-1}-...-\phi_pX_{t-p}=Z_t+\theta_1Z_{t-1}+...+\theta_qZ_{t-q}.$$ Hier sind die $\phi's$ ($\theta's$), die Paramter welche den autoregressiven (moving average) Prozess charakterisieren. Die $Z's$ sind Rauschterme. Als Standardspezifikation sind diese standardnormalverteilt. In der Funktion arma_sim können diese Paramter ($\phi$'s,$\theta$'s) in den gleichnamigen Inputvariablen jeweils als numerische Vektoren angegeben werden. Es muss mindestens eine der beiden Variablen (phi, theta) spezifiziert werden. Zum zweiten muss die Länge der resultierenden Zeitreihe in der Inputvariable (n) spezifiziert werden.

Der folgende Code simuliert 100 Beobachtungen eines ARMA(2,2)-Prozesses:

set.seed(12345)
sim_obj <- arma_sim(phi = c(0.5,-0.2), 
                    theta = c(0.1,-0.3), 
                    n = 100)
str(sim_obj)

Hier ist auch die Struktur des Outputs zu erkennen. Es wird ein Objekt der Klasse "arma" ausgegeben, welches auf einer Liste beruht. Unter innov sind die Rauschterme der jeweiligen Periode gespeichert. Auch die Inputparamter des ARMA-Prozesse können abgerufen werden. Damit die reultierende Zeitreihe weniger von den Startwerten abhängt, kann eine "Burn-in" Periode bestimmt werden in der sich der Prozess einpendeln kann. Die Länge dieser Periode ist in dem Objekt burnin im Output abgespeichert und kann in der gleichnmaigen Inputvariable spezifiziert werden.

Der folgende Plot zeigt, die oben simulierte Zeitreihe und wird mit der Funltion plot_ts erstellt. Diese Funktion ermöglicht einfache Linienplots von Zeitreihen und reagiert auf "arma"-Objekte.

library(ggplot2)
library(tibble)

plot_ts(sim_obj)

Man kann in diesem Plot deutlich erkennen, dass die Funktion arma_sim standardmäßig Zeitreihen mit einem Mittelwert von Null erzeugt. Ein anderer Mittelwert kann mit der Variable mu spezifiziert werden. Im folgenden ist zu erkennen, dass sich der Mittelwert nach der Eingabe von mu verändert.

set.seed(12345)
sim_obj2 <- arma_sim(phi = c(0.5,-0.2), 
                     theta = c(0.1,-0.3), 
                     mu = 1, 
                     n = 100)
## Mittelwert (mu) auf 1 gesetzt

plot_ts(sim_obj2)

Als letztes ist es in der Funktion möglich die Verteilung der Innovations bzw. auch die Innovations an sich zu spezifizieren. Durch die Inputvariable innov kann eine eigene Reihe von Innovations übergeben werden (Das war beim testen der Funktion nützlich). Die zweite Möglichkeit ist es die Verteilung der Innovations in innov.gen zu spezifizieren, und die Verteilungsparameter in ... zu übergeben. Im folgenden ist ein Plot mit höherer Standardabweichung der Innovations zu sehen. Die Verteilung wird bei der Normalverteilung belassen, könnte jedoch auch geändert werden.

set.seed(12345)
sim_obj3 <- arma_sim(phi = c(0.5,-0.2), 
                     innov.gen = rnorm,  
                     theta = c(0.1,-0.3), 
                     mu = 1, 
                     n = 100, 
                     sd=2)

plot_ts(sim_obj3)

Durch das setzten des "Seeds" ist die Form der Zeitreihe gleich zum vorherigen Plot. Durch die höhere Standardabweihung sind die Ausschläge der Zeitreihe jedoch klar größer.

Autokovarianz-Funktion

Die Autokovarianz-Funktion gibt die Kovarianz der Realisationen eines Prozesses mit sich selbst zu früheren Zeitpunkten zurück. Damit ist sie ein erster Anhaltspunkt zu den zeitlichen Abhängigkeiten des Prozesses. Die Autokovarianz-Funktion ist nach Brookwell et al. (2002) für den Lag $h$ folgendermaßen definiert:

$$\gamma_x(h)=Cov(X_{t+h},X_t) $$ Die Funktion acf gibt standardmäßig die Stichprobenautokovarianzen für den Lag $0$ bis Lag $n-1$ zurück. In der Inputvaribale lag.max kann die Anzahl der zurückgegebenen Autokovarianzen begrenzt werden. Die Funktion acf extrahiert ohne weiter Eingaben die Zeitreihe aus simulierten arma-Objekten.

set.seed(123)

ar_obj <- arma_sim(phi = c(0.9), n = 100)
ar_acf <- acf(ar_obj, lag.max = 20)

ma_obj <- arma_sim(theta = c(0.5), n = 100)
ma_acf <- acf(ma_obj, lag.max = 20)

ar_plot <- ggplot(tibble(ACF=ar_acf, Lag=seq_along(ar_acf)), aes(x=Lag, y=ACF)) +
                geom_bar(stat="identity") + 
                ggtitle("AR(1)")

ma_plot <- ggplot(tibble(ACF=ma_acf, Lag=seq_along(ma_acf)), aes(x=Lag, y=ACF)) +
                geom_bar(stat="identity") + 
                ggtitle("MA(1)")
ar_plot
ma_plot

Anhand der beiden Autokovarianz-Funktionen kann man eine erste Aussage über die Art des Prozesses treffen. Der erste Plot zeigt eine langsam abnehmende Autokovarianz-Funktion, was auf einen AR-Prozess hindeutet. Der zweite Plot zeigt eine abprupt abnehemnde Autokovarianz-Funktion, was eher auf einen MA-Prozess schließen lässt.

Durbin-Levinson Algorithmus

Der Durbin-Levinson Algorithmus kann zur Vorhersage von Zeitreihen verwendet werden. Nach Brookwell et al. (2002) kann der Algorithmus verwendet werden um die $\phi's$ in folgender Formel zu schätzen: $$ P_nX_{n+1}=\phi_{n1}X_n+...+\phi_{nn}X_1. $$ Das heißt der vorgesagte Wert wird als lineare Kombination aller vergangenen Realisiationen der Zeitreihe ausgedrückt. Die $\phi's$ werden rekursiv durch folgende Formeln bestimmt: $$ \phi_{nn}=\Big[ \gamma(n)-\sum \limits_{j=1}^{n-1}\phi_{n-1,j}\gamma(n-j)\Big]v_{n-1}^{-1}, $$ $$\begin{bmatrix}\phi_{n1} \ \vdots \\phi_{n,n-1} \end{bmatrix}= \begin{bmatrix}\phi_{n-1,1} \ \vdots \\phi_{n-1,n-1} \end{bmatrix}- \phi_{nn} \begin{bmatrix}\phi_{n-1,n-1} \ \vdots \\phi_{n-1,1} \end{bmatrix},$$ und $$v_n=v_{n-1}[1-\phi_{nn}^2] $$ Die Startwerte sind $\phi_{11}=\gamma(1)/\gamma(0)$ und $v_0=\gamma(0)$.

Die Funktion DL implementiert den Durbin-Levinson Algorithmus und gibt die geschätzten $\phi's$ zurück. Diese können dann zur Vorhersage verwendet werden. Standardmäßig führt die Funktion die maximale Anzahl von Rekursionen durch (n-1), jedoch kann diese in der Inputvariable p reduziert werden. Damit werden dann auch wenniger $\phi's$ geschätzt. Auf die Auswirkungen von kleinerem p wird in der Simulationsstudie genauer Bezug genommen. Folgender Code zeigt die Anwendung auf einen AR-Prozess, für welche der Algorithmus relativ gut funktioniert. Es wird beim Simulieren eine lange Burn-in Periode gewählt, um zu garantieren, dass der Prozess wenig von den zufälligen Startwerten abhängt.

set.seed(123)
ar_obj <- arma_sim(phi = c(0.6,0.3), n = 100, burnin = 2000)
dl <- DL(ar_obj)

# p=100 (hier: Default)
dl

# p=5
DL(ar_obj, p=5L)

Der Output zeigt für beide Eingaben von p die geschätzten $\phi's$. Die letzten beiden Werte entsprechen jeweils in etwa den Inputparametern des AR-Prozesses. Die restlichen $\phi's$ sind nahe null. Somit wäre auch eine dementsprechend gute Vorhersage für den nächsten Zeitpunkt möglich. Die Vorhersage wird durch das Multiplizieren der $\phi's$ mit der Zeitreihe und anschließendes Summieren bestimmt.

forecast <- sum(ar_obj$arma*dl)
forecast

Innovation Algorithmus

Der Innovation Algorithmus kann auch wie der Durbinson Levinson Algorithmus zur Vorhersage von Zeitreihen benutzt werden. Aber im Gegensatz zum DL-Algorithmus ist der Innovation Algorithmus gut zum Simulieren von MA Prozesse und anstatt von AR Prozesse. Der Algorithmus schätzt die $\theta's$ für den folgenden Algorhitmus (nach Brookwell et al. (2002)): $$ P_0=0. $$ $$ P_{n}=\hat{X}{n+1}=\sum{j=1}^{n}\theta_{n,j}(X_{n+1-j}-\hat{X}{n+1-j}) $$ Die Vorhersage des ${X}{n+1}$ wird durch die $\theta's$ des Innovation Alg, die vorherigen Werte und deren Approximationen bestimmt. Zur Bestimmung der $\theta's$ braucht man zusätzlich die mittlere quadratische Abweichung(Simultan Bestimmung im Algorithmus) $v_i= E(X_{i+1}-\hat{X}{i+1}$ und die ACF-Funktion der Werte $\kappa$. $$v_0=\kappa(0)$$ $$\theta{n,n-k}=v_{k}^{-1}\biggl(\kappa(n-k)-\sum_{j=0}^{k-1} \theta_{k,k-j}\theta_{n,n-j}v_j \biggr).\quad 0\ge k\lt n$$ $$v_n=\kappa(0)-\sum_{j=0}^{n-1}\theta_{n,n-j}^2v_j$$

Die Funktion innovation implementiert den innovation Algorithmus und gibt die geschätzten $\theta's$ und v's zurück. Dieser Algorhithmus wird von der innovation_prediction Funktion, wie oben beschrieben, zur Vorhersage des nächsten Wertes benutzt. Standardmäßig führt die Funktion die maximale Anzahl von Rekursionen durch (n-1), jedoch kann diese in der Inputvariable lag.max reduziert werden. Folgender Code zeigt die Anwendung auf einen MA-Prozess. Es wird beim Simulieren eine lange Burn-in Periode gewählt, um zu garantieren, dass der Prozess wenig von den zufälligen Startwerten abhängt.

set.seed(1234)
ma_obj <- arma_sim(theta = c(0.6,-0.2), n = 500, burnin = 2000)
innov <- innovation(ma_obj)
print(innov[[1]][1:10,1:10])
print(innov[[2]][1:10])

Man sieht das die $\theta_{.,1}$ nahe an 0.6 sind und die $\thetas_{.,2}$ nahe an -0,2 sind. Die restlichen $\theta's$ sind im Vergleich dazu relativ klein. In der Theorie sollten diese gleich Null sein. Dies schafft diese Implemantion nicht. Die Vohersage wird dann durch die innovation_prediction getätigt, wobei hier eine Ein-Schritt Vorhersage getätigt wird.

innovation_prediction(ma_obj)

Innovation Prediction

Diese Funktion benutzt den Innovation Algorithmus um die Vorhersage für das nächsten (bzw. für das n-te) Wert zu bestimmen. Nach Brockwell sind dies die Algorithmen dazu: $$ P_0=0. $$ $$ P_{n}=\hat{X}{n+1}=\sum{j=1}^{n}\theta_{n,j}(X_{n+1-j}-\hat{X}_{n+1-j}. . $$

bzw. zur Berechnung von $P_{n+h-1}$: $$ P_{n+h-1}=\hat{X}{n+h}=\sum{j=h}^{n+h-1}\theta_{n+h-1,j}(X_{n+h-j}-\hat{X}_{n+h-j}). $$

Da aber für die Vorhersage des n-te Wertes $\theta$ bis zu $\theta_{n+h-1,n+h-1}$ berechnet werden muss, werden die ersten n Werte nur zur Berechnung der $\theta's$ benutzt werden.

Die Funktion hat als neben der Zeitreihe ts, noch die optionalen Eingabenwerte steps und lag.max. Wobei steps (Standardwert = 1) angibt welches steps-te Element man vorhersagen will und lag.max (Standardwert: length(ts)-steps), wie viele Elemente der Zeitreihe benutzt werden, wobei die ersten step Schritte immer benötigt werden um die $\theta's$ des letzten Schritt zu berechnen. <<<<<<< HEAD Folgender Beispiel Code zeigt die Benutzung:

innovation_prediction(ma_obj)
innovation_prediction(ma_obj,steps=2)

Periodogramm

Das Periodogramm einer Zeitreihe ${x_1,...,x_n}$ ist eine stichprobenbasierte Funktion, aus der ein nicht-konsistenter Schätzer für die spektrale Leistungsdichte gewonnen werden kann.

Das Periodogramm einer Zeitreihe ${x_1,...,x_n}$ ist definiert als: $$ I_n(\lambda)=\frac{1}{n} \Bigg| \sum \limits_{t=1}^{n}x_{t}e^{-it\lambda} \Bigg| ^2 $$

Das Periodogramm folgt der Idee ein Signal in sinusartige Funktionen mit unkorrelierten Koeffizienten zu zerlegen. $$ x_t=\sum \limits_{j \in F_n}a_j [\cos(\omega_jt) + i \sin(\omega_jt)] $$ mit$\omega_j = \frac{2 \pi j}{n}$, den Fourier Frequenzen der Zeitreihe ${x_1,...,x_n}$, und $j \in F_n={j \in \mathbb{Z} : -\pi < \omega_j \leq \pi }$. Man beachte, dass $F_n$ n Elemente enthält.

Die Folge ${a_j: j\in F_n}$ wird als diskrete Fourier-Transformation der Zeitreihe ${x_1,...,x_n}$ bezeichnet. $$ a_j=\frac{1}{\sqrt{n}}\sum \limits_{t=1}^{n}x_t e^{-it\omega_j} $$

Die Funktion periodogram implementiert das Periodogramm und gibt einen numerischen Vektor zurück, welcher die Ausprägungen des Periodogramms an den Stellen $\omega_j$ ($j \in F_n$) enthält. $$ I(\omega_j)=\frac{1}{n} \Bigg| \sum \limits_{t=1}^{n}x_{t}e^{-it\omega_j} \Bigg| ^2 $$

Folgender Code beschreibt die Anwendung der Funktion in R:

x <-  arma_sim(phi = c(1,-0.6), n = 1000, burnin = 500)
periodo <- periodogram(x)

plot_periodogram(periodo)

Da für die Fourier-Frequenzen gilt $I(-\omega_j) = I(\omega_j)$, bietet es sich an nur eine Seite des Periodogramms in den Graphen aufzunehmen. Die Funktion plot_periodogram führt genau diese graphische Aufbereitung der Ergebnisse durch. Sie stellt die Werte des Periodogramms an den Fourier Frequenzen graphisch dar.

Anhand der graphischen Aufbereitung lassen sich auch die Ergebnisse der Funktion periodogram interpretieren. Wenn eine Zeitreihe glatt (bzw. verwackelt) erscheint, dann werden die Werte des Periodogramms für niedrige (bzw. hohe) Frequenzen im Verhältnis zu seinen anderen Werten groß sein, und man sagt, dass der Datensatz einen Überschuss an niedrigen (bzw. hohen) Frequenz aufweist. Bei einer rein zufälligen Serie sollten alle Sinus- und Cosinuskurven von gleicher Bedeutung sein, und daher wird das Periodogramm zufällig um eine Konstante variieren.

Sollte eine Zeitreihe ein starkes sinusartiges Signal für eine Frequenz aufweise, dann wird das Periodogramm eine Spitze bei dieser Frequenz aufweisen. Falls die Zeitreihe ein starkes nicht-sinusbasiertes Signal für eine Frequenz aufweist, wird das Periodogramm eine Spitze in dieser Frequenz und den Vielfachen aufweisen.