In leopopeo/Time-Series: Analysen von Zeitreihen

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

Einleitung

In dieser Vignette werden die Funktionen des Packages TimeSeries ausführlich erklärt. Das Package enthält Implementationen einiger wichtiger Algorithmen zur Analyse und Vorhersage von sogenannten ARMA-Prozessen aus Introduction to Time-Series Analysis and Forecasting (Brockwell et. al, 2002):

• ARMA Generator für AR(p) und MA(q) Prozesse arma_sim
• Autokovarianz-Funktion für Zeitreihen ACF
• Durbin-Levinson-Algorithmus für Zeitreihen DLA
• Innovations Algorithmus innovation und Vorhersagen mit diesem ts_predict
• Periodogramm Algorithmus für Zeitreihen perio
• Plot-Funktion für Zeitreihen plot_timeseries und Periodogramme plot_periodogram

Um die Funktionen nutzen zu können, muss zuerst das Paket TimeSeries geladen werden. Die anderen Pakete werden zu Demonstrationszwecken geladen.

library(TimeSeries)
library(ggplot2)
library(tibble)

ARMA Generator

Ein ARMA(p,q)-Prozess (Autoregressive Moving Average) ist ein stochastischer Prozess der folgendem Modell genügt: ${X_t}$ ist stationär und für alle $t$ gilt:

$$ X_t = \sum_{i= 1}^{p} \phi_i X_{t-i} + \sum_{j = 1}^{q} \theta_j Z_{t-j} + Z_t, \qquad\ Z_t \sim \mathcal{N}(0, \sigma^2) $$ und die Polynome $1 + \phi_1X + ... + \phi_p X^p$ und $1+ \theta_1X + ... + \theta_q X^q$ sind teilerfremd (Brockwell et. al (2002), S. 83).

Hier zunächst einige Beispiele für ARMA-Prozesse.

1. Als erstes sehen wir einen ARMA(1,1)-Prozess mit einer kleinen und einer großen Standartabweichung. Je größer die Standartabweichung ist, desto größer sind die Schwankungen des Rauschterms um den Mittelwert. Dieser ist bei uns als Defaultwert mit 0 festgelegt. Der Plot nutzt dabei die Funktion plot_timeseries und erzeugt einen einfachen Linienplot der Zeitreihe.

set.seed(3)
arma_klein <- arma_sim(phi = 0.7, theta=0.3,sd=0.1,I=100)
plot_timeseries(arma_klein)

set.seed(3)
arma_groß <- arma_sim(phi = 0.7, theta=0.3,sd=1,I=100)
plot_timeseries(arma_groß)

1. Hier sehen wir einen AR(1)-Prozess mit einem kleinen und einem großen Phi-Wert. Man erkennt, dass der Phi-Wert mit Glattheit der Zeitreihe zusammenhängt: je größer der Wert, desto glatter ist die Funktion. Gleiches gilt für einen MA(1)-Prozess mit großen und kleinen Theta-Werten.

set.seed(3)
arma_klein <- arma_sim(phi = 0.3,sd=1,I=100)
plot_timeseries(arma_klein)

arma_groß <- arma_sim(phi = 0.9,sd=1,I=100)
plot_timeseries(arma_groß)

Autokovarianz-Funktion

Der Begriff der Autokovarianz stammt aus der Signalverarbeitung und gibt die Kovarianz einer Funktion bzw. eines Signals mit sich selbst zu einem früheren Zeitpunkt an. Die Funktion drückt also aus, wie viel Ähnlichkeit die um eine gewisse Zeitdifferenz verschobene Folge mit der ursprünglichen Folge besitzt. Die Autokovarianz-Funktion ist nach Brockwell et al. (2002) (S. 19) für eine Zeitdifferenz h wie folgt definiert:

$$\gamma_x(h):=COV(X_{t+h},X_t) $$ Da die obige Funktion jedoch nur funktioniert, wenn das der Zeitreihe zugrunde liegende Modell bekannt ist arbeitet unsere Funktion nur mit den gemessenen Daten. Die Funktion ACF gibt also für eine Zeitreihe $x_{1},...,x_{n}$ die geschätzte Sichtprobenkovarianz für die Zeitdifferenz lag = $0$ bis lag = $n-1$ an. Zusätzlich zur Input Zeitreihe kann die lag-Variable individuell angegeben werden. Es kann also die Anzahl der zurückgegebenen Autokovarianzen gewählt werden. Dafür wird folgende Formel genutzt:

$$ \hat{\gamma}(h) := n^{-1} \sum_{t=1}^{n-|h|} (x_{t+|h|}-\bar{x})(x_t-\bar{x})$$ ````r set.seed(3) ar <-arma_sim(phi = c(0.2,0.6),I=100,sd=0.01)

ar <- arima.sim(n = 100,list(ar = c(0.2,0.6)))

ar_acf <- ACF(ar, lag = 20) ma <- arma_sim(theta = c(0.5,0.9),I=100,sd=0.01)

ma <- arima.sim( n = 100,list(ma = c(0.5,0.9)))

ma_acf <- ACF(ma, lag = 20) ar_plot <- ggplot(tibble(ACF=ar_acf, Lag=seq_along(ar_acf)), aes(x=Lag, y=ACF)) + geom_bar(stat="identity") + ggtitle("AR(2)") ma_plot <- ggplot(tibble(ACF=ma_acf, Lag=seq_along(ma_acf)), aes(x=Lag, y=ACF)) + geom_bar(stat="identity") + ggtitle("MA(2)") ar_plot ma_plot

Anhand der beiden Autokovarianz-Funktionen können wir eine erste Aussage über die Prozessart 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 auf einen MA-Prozess schließen lässt.

## Durbin-Levinson Algorithmus

Der Durbin-Levinson Algorithmus kann für Vorhersagen einer Zeitreihe genutzt werden. Der Algorithmus benutzt dabei eine Linearkombination aller vorherigen Zeitreihen-Punkte, um einen neuen Wert vorherzusagen. Dies wird durch die folgende Gleichung nach Brookwell et al. (2002) ausgedrückt:

$$\begin{align} 
P_nX_{n+1}=\phi_{n1}X_n+...+\phi_{nn}X_1.
\end{align}$$

Die entprechenden $\phi$ werden durch den folgenden Algorithmus rekursiv 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]. $$

Für die Startwerte gilt $\phi_{11}=\gamma(1)/\gamma(0)$ und $v_0=\gamma(0)$.

Die Funktion *DLA* ist die Implementation des Durbin-Levinson Algorithums und gibt die geschätzten $\phi$ zurück. Die Funktion benutzt standartmäßig die maximale Anzahl an Rekursionen, also n-1. Dieser Wert kann individuell mit dem *len*-Parameter gesetzt werden. Die erhaltenen $\phi$ können nun für eine Vorhersage benutzt werden. \n
Im folgenden Beispiel wird die Funktion auf AR-generierte Daten angewendet:

```r
set.seed(2)
ar <- arma_sim(phi = c(0.2,0.6),I=100,sd=0.01)
res <- DLA(ar)

# len=100 (hier: Default)
res

# len=5
DLA(ar, len=5)

Der Output zeigt für die beiden angegebenen Rekursionen die geschätzten $\phi$ an. Die ersten beiden Werte entsprechen hier nahezu dem angegebenen Parameter phi in arma_sim. Für die Vorhersage müssen nun die $\phi_{n1},...,\phi_{nn}$ mit der Zeitreihe (außgenommen dem letzten Wert, da die maximale Rekusionszahl n-1 ist) multipliziert und anschließend summiert werden.

res <- sum(ar*rev(res))
res

Dies wird nun auch im DL_prediction Algorithmus realisiert.

Durbin-Levinson Prediction

Die Funktion DL_predict nutzt obig beschriebenes Verfahren, um Vorhersagen für eine Zeitreihe basierend auf dem Durbin-Levinson Algorithmus zu treffen. Der Eintrag $X_{n+1}$ der Zeitreihe $X_1, ... , X_{n}$ wird dabei jeweils durch:

$$ X_{n+1}=\phi_{n1}X_n+...+\phi_{nn}X_1. $$ bestimmt. Mit Hilfe der Funktion plot_timeseries können wir als zweites Argument die Vorhersage übergeben und bekommen einen Linienplot ausgegeben. Dies sind für einen AR(1)-Prozess folgendermaßen aus:

set.seed(13)
arma <- arma_sim(phi=0.5,sd=1,I=150)
predict <- DL_prediction(arma,len=15)
plot_timeseries(arma,predict)

Auch für einen ARMA(2,2)-Prozess kann die Funktion DL_prediction verwendet werden.

set.seed(13)
arma <- arma_sim(phi=c(0.3,0.5),theta=c(0.2,0.4),sd=1,I=150)
predict <- DL_prediction(arma,len=15)
plot_timeseries(arma,predict)

Innovation Algorithmus

Der Innovation Algorithmus ist eine zweite Möglichkeit, um Werte der Zeitreihe vorherzusagen. Sei $\mathbf{\hat X_n} = (X_1,\ P_1 X_2,\ ...\ ,\ P_{n-1}X_n)$. Dann gibt es eine Matrix $\mathbf{ \Theta_n}$, sodass $\mathbf{\hat X_n} = \mathbf{ \Theta_n} (\mathbf{X_n - } \mathbf{\hat X_n})$. Es gilt also $$\begin{equation} \hat{X}{n+1}= \begin{cases}0, & \text { if } n=0, \ \sum{j=1}^{n} \theta_{n j}\left(X_{n+1-j}-\hat{X}_{n+1-j}\right), & \text { if } n=1,2, \ldots,\end{cases} \end{equation}$$ (ebd. S.72)

Wir können nun die Koeffizienten $\theta_{n,i}$ rekursiv berechnen, wobei $\kappa(i,j)$ die Kovarianz ist.

$$ \begin{aligned} v_{0} &= \kappa(1,1), \ \theta_{n, n-k} &= v_{k}^{-1}\left(\kappa(n+1, k+1)-\sum_{j=0}^{k-1} \theta_{k, k-j} \theta_{n, n-j} v_{j}\right), \quad 0 \leq k<n,\ v_{n} &= \kappa(n+1, n+1)-\sum_{j=0}^{n-1} \theta_{n, n-j}^{2} v_{j} \end{aligned} $$ (ebd. S.73). Die Funktion innovation(X, lag = n) berechnet auf diese Weise $\mathbf\Theta_n$, welches wir für die Vorhersage mit ts_predict verwenden. Mit dem Parameter lag $<$ length(X) kann das $n$ gewählt werden, bis zu welchem $\theta_{n,i}$ berechnet wird.

X = arma_sim(theta = c(0.5,0.9),I=5,sd=0.01)
innovation(X, lag = 4)

Innnovation Prediction

Obige Gleichung können wir nun auch verwenden, um eine $h$-Schritt-Prognose des Prozesses zu machen. Sei dazu $m$ die Zeitreihenlänge. Dann berechnen wir mit $\hat X_{m+h}$ eine h-Schritt Prognose. Die Ausgabe von ts_predict ist per Default der Vektor $(\hat X_{m+1},..., \hat X_{m+h})$ der Länge $h$ mit den Prognosewerten. Wird all = TRUE gesetzt, wird der gesamte Vektor $\mathbf{\hat X_{m+h}}$ ausgegeben. Die Vorhersage mit dem Innovation-Algorithmus eignet sich besonders gut für $MA(q)$-Prozesse mit kleinem $q$. Hier wird die ts_predict-Funktion auf einen $MA(2)$-Prozess angewendet:

I = 50
h = 10
X = arma_sim(theta = c(0.5,0.4),I = I,sd=0.01)
X_hat = ts_predict(X, h)
plot_timeseries(X, pred = X_hat)

Periodogramm

Das Periodogramm einer Zeitreihe $(X_1,...,X_n)$ ist deren Spektraldichte. Die Definition ist dabei nach Brockwell et al. (2002) (S.123):

$$ I_n(\lambda)=\frac{1}{n} \Bigg| \sum \limits_{t=1}^{n}x_{t}e^{-it\lambda} \Bigg| ^2 $$ Dabei steht der Koeffizient $x_t$ für eine Kombination trigonometrischer Funktionen:

$$ x_t=\sum \limits_{k \in F_n}a_k [\cos(\omega_kt) + i \sin(\omega_kt)] $$ wobei $F_n$ für die Fourierfrequenzen der Zeitreihe im Intervall $(-\pi,\pi]$ steht und für $\omega_k$ gilt: $$\omega_k=\frac{2\pi k}{n}, \space\space\space\space\space\space\space\space k = [-\frac{n-1}{2},...,\frac{n}{2}]$$ Die Folge ${a_k}$ ist die diskrete Fouriertransformation der Zeitreihe und berechet sich durch:

$$ a_k=\frac{1}{\sqrt{n}}\sum \limits_{t=1}^{n}x_t e^{-it\omega_k} $$ Die Funktion perio implementiert nun den Periodogramm-Algorithmus und gibt für die $\omega$ im oben definierten Intervall den Wert $I_n(\lambda)$ zurück. Der folgende Code veranschaulicht dies:

x <- arma_sim(phi = c(0.2,-0.6),I=1000,sd=0.01)
periodo <- perio(x)
plot_periodogram(periodo) 

Quellen

Brockwell, Peter J. und Davis, Richard A. (2002): Introduction to Time Series and Forecasting (2nd Edition), Springer: New York (u.a.).

leopopeo/Time-Series documentation built on Dec. 21, 2021, 10:42 a.m.