Zurück...

Zeitreihenanalyse-Beispiele mit R

Lineares Modell

Die hier verwendete Beispiel-Datei Zeitreihe_Getreide.csv beinhaltet 191 Bobachtungen und stammt vom Statistischen Bundesamt. Die Datei Zeitreihe_Getreide.csv ist ein Auszug und dient als Beispiel-Datei:

> Getreide <- read.csv2("Zeitreihe_Getreide.csv")
> Getreide
    
Preisindex
1    133.3
2    133.2
3    130.8
4    132.2
5    135.6
6    137.5
7    136.7
8    136.0
9    137.1
10    134.2
11

Über die Funktion ts() wird die Variable Getreide in ein Zeitreihen-Objekt (Time-Series...) gewandelt:

> GetreideIndex <- ts(Getreide$Preisindex, freq = 12, start = 1995)
> GetreideIndex

    Jan  Feb  Mar  Apr  May  Jun  Jul  Aug  Sep  Oct  Nov  Dec
1995 133.3 133.2 130.8 132.2 135.6 137.5 136.7 136.0 137.1 134.2 136.6 137.5
1996 138.2 137.5 138.4 148.2 151.8 146.4 141.4 135.4 130.2 123.7 121.7 122.6
1997 123.4 124.3 125.3 126.9 125.4 123.2 121.7 123.2 120.9 119.3 120.1 120.8
1998 120.2 118.3 117.9 115.3 114.7 113.8 111.4 111.1 107.6 105.6 108.5 108.1
1999 108.1 108.3 109.5 112.2 113.7 115.5 112.5 113.1 110.8 109.7 110.6 111.8
2000 112.8 114.0 115.5 117.8 117.2 116.2 116.8 117.1 116.6 116.6 116.5 116.6
2001 117.5 117.5 116.6 115.8 116.9 117.9 117.6 116.1 114.4 113.5 115.1 115.4
2002 115.4 115.2 114.5 113.7 111.7 109.0 107.9 110.9 110.9 111.6 110.9 110.0
2003 109.3 108.8 107.3 107.4 106.1 106.8 109.6 112.2 119.1 121.9 130.4 131.8
2004 136.2 134.0 132.3 132.4 122.4 119.7 115.2 108.3 104.7 103.6 103.9 99.4
2005 100.1 98.1 97.4 97.0 96.8 99.8 102.0 101.4 102.4 101.0 102.0 101.7
2006 101.9 101.5 101.7 103.2 105.1 105.5 107.5 112.8 117.4 124.1 127.6 129.0
2007 129.0 128.9 128.2 130.0 129.9 133.3 142.2 166.2 171.6 184.7 186.2 193.9
2008 199.0 199.5 200.9 187.7 177.8 180.4 173.8 165.4 154.2 140.9 134.3 129.7
2009 130.3 128.9 124.9 125.5 128.9 127.8 123.4 119.4 112.1 113.8 114.9 119.5
2010 119.0 117.3 118.8 119.0 122.6 129.1 146.2 164.4 167.8 169.0 178.6 187.2
>

Wie Sie oben erkennen können, wurde als Startpunkt 1995 (start = 1995) gewählt. Da jede Beobachtung für einen Monat steht, wurde die Jahreseinteilung auf 12 Monate (freq = 12) festgelegt.

Verschaffen wir uns einen grafischen Überblick:

> plot(GetreideIndex, main = "Getreide-Preisindex", xlab = "Zeit", ylab = "Preisindex")

In obiger Abbildung sind deutliche Schwankungen und Spitzen zu erkennen. Auch wenn ein lineares Modell aufgrund obiger Grafik nicht das optimale Modell ist, fangen wir trotzdem bei unserer Zeitreihenanalyse mit diesem Modell an!

Dazu erstellen wir einen Vektor t über die Anzahl der Beobachtungen (1 bis 192 Monate) in dem Zeitreihenobjekt GetreideIndex:

> t <- 1:length(GetreideIndex)
> t
 [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30
 [31] 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
 [61] 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83  84 85 86 87 88 89 90
 [91] 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120
[121] 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150
[151] 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180
[181] 181 182 183 184 185 186 187 188 189 190 191 192

Das lineare Modell wird über die Parameter GetreideIndex und t geschätzt:

> Linear <- lm(GetreideIndex ~ t)
> summary(Linear)

Call:
lm(formula = GetreideIndex ~ t)

Residuals:
  Min   1Q Median   3Q   Max
-30.616 -12.939 -5.476  6.482 70.258

Coefficients:

       Estimate Std. Error t value Pr(>|t|)  

(Intercept) 115.55472  3.05700 37.800 < 2e-16 ***
t       0.09489  0.02747  3.454 0.00068 ***

---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 21.1 on 190 degrees of freedom
Multiple R-squared: 0.05909,  Adjusted R-squared: 0.05414
F-statistic: 11.93 on 1 and 190 DF, p-value: 0.0006805

Wie schon vermutet, ist die Anpassung für das lineare Modell

Getreide-Preisindex= 115,55 + 0,0948 * t

mit dem Gütemaß Multiple R-squared: 0.05909 nicht überragend!

Schauen wir uns die Anpassung des Modells grafisch an! Dazu verwenden wir die Funktion predict und zur visuellen Beurteilung der Modellgüte lassen wir den Vertrauensbereiches mit P=95% (level = 0.95) mit schätzen:

> Linear_P <- predict(Linear, int="c", level = 0.95)

Zuerst stellen wir wieder die Grafik Getreide-Preisindex ...,

> plot(GetreideIndex, main = "Getreide-Preisindex", xlab = "Zeit", ylab = "Preisindex" )

… dann die Schätzung (Anpassung, fit) des Preisindexes (abgelegt in Linear_P[,1], dargestellt in rot), ...

> lines(ts(Linear_P[,1], freq = 12, start = 1995), col = "red")

… und anschließend den unteren Vertrauensbereich (Linear_P[,2], grün) mit ...

> lines(ts(Linear_P[,2], freq = 12, start = 1995), col = "green", type="p")

... und den oberen Vertrauensbereich (Linear_P[,3], blau) mit ...

> lines(ts(Linear_P[,3], freq = 12, start = 1995), col = "blue", type="p")

... dar:

Die obige Grafik spricht für sich! Das lineare Modell ist nicht geeignet, um Vorhersagen für Zeitreihen wie abgebildet machen zu können.

Lineares Modell, glätten

Vielleicht wird das lineare Modell besser, wenn wir versuchen, die extremen Bobachtungen zu glätten (oder auch filtern genannt). Für das nächste Beispiel wird die Methodik des gleitenden Durchschnitts angewendet. Dazu nutzen wir aus dem Paket zoo die Funktion rollmean():

> Glaetungsfaktor <- 3 # 3 Beobachtungen werden zu einem Wert gemittelt
> Getreide_g3 <- rollmean(Getreide, Glaetungsfaktor)

> length(Getreide$Preisindex)
[1] 192

> length(Getreide_g3)
[1] 190

Nach der Glättung über 3 Beobachtungen ist der Getreide_g3-Vektor erwartungsgemäß um 2 Beobachtungen reduziert.

Führen wir für die geglätteten Getreidedaten die Regressionsanalyse nach bekanntem Modell durch:

> t_g3 <- 1:length(Getreide_g3)
> Linear_g3 <- lm(Getreide_g3 ~ t_g3)
> summary(Linear_g3)

Call:

lm(formula = Getreide_g3 ~ t_g3)

Residuals:
  Min   1Q Median   3Q   Max
-29.743 -12.588 -5.265  7.150 69.952

Coefficients:
       Estimate Std. Error t value Pr(>|t|)  
(Intercept) 115.81944  2.99617 38.656 < 2e-16 ***
t_g3     0.08935  0.02721  3.284 0.00122 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 20.57 on 188 degrees of freedom
Multiple R-squared: 0.05426,  Adjusted R-squared: 0.04923
F-statistic: 10.79 on 1 and 188 DF, p-value: 0.001219

Auch hier geben wir das Modell wieder grafisch aus:

> plot(Getreide_g3, main = "Getreide-Preisindex", sub = "Glättungsfaktor = 3", xlab = "Zeit", ylab = "Preisindex", type = "b", col = "red")
> lines(fitted(Linear_g3), col = "blue")

Die lineare Funktion ist durch das Glätten nicht besser geworden.

Exponentielles Glätten (Holt-Winters)

Eine weitere Glättungsmöglichkeit ist das exponentielle Glätten. Hierbei wird über Parameter gesteuert, mit welchem Gewicht aktuelle Beobachtungen in das Modell einfließen und welches Gewicht den älteren Beobachtungen zugestanden wird. Eine derartige Glättungsfunktion ist die Holt-Winters-Funktion (in R: HoltWinters()).

Bevor wir uns mit dieser Funktion beschäftigen, schauen wir uns die Funktion decompose() an.

Ungeachtet der linearen Modell-Güte ist eine saisonale Komponente in der Zeitreihe Getreide-Preisindex visuell nicht zu erkennen. Deswegen wird erwartet, dass die R-Funktion decompose() die die Zeitreihen in die Trend, Saison und irreguläre Komponente (Restkomponente) zerlegt, kein Resultat liefert:

> Zerlegung <- decompose(Getreide$Preisindex)
Fehler in decompose(Getreide$Preisindex) :
 Zeitreihe hat keine oder weniger als 2 Perioden

Unsere Erwartung wurde bestätigt. Kommen wir zurück zum exponentiellen Glätten unter Verwendung der Holt-Winters-Funktion.

Die Idee, die hinter der exponentiellen Glättung steht, ist besonders für ökonomische Zeitreihen einsichtig: Ist es  sinnvoll, allen Beobachtungen der Zeitreihe das gleiche Gewicht einzuräumen, oder ist es sinnvoller jüngeren Beobachtungen mehr Gewicht als älteren Bobachtungen einzuräumen? Wenn Sie für Ihre Zeitreihe diesem Gedanken zustimmen können, ist die Methodik des exponentiellem Glättens wahrscheinlich die Richtige für Sie!

In dieser Methodik werden die Gewichte der Beobachtungen mit jeder neu hinzukommenden Beobachtung exponentielle „abgewichtet“:

Ändern sich die Beobachtungen langsam, sollte gw groß gewählt werden. Dadurch werden die Beobachtungen etwa gleich gewichtet. Ändern sich die Beobachtungen hingegen schnell, sollte gw klein gewählt werden. Dadurch  liegt das Gewicht bezüglich der Beobachtungen in der Zeitreihe auf den jüngsten Beobachtungen. Für eine große Anzahl Beobachtungen kann in obiger Formel der Term wt unberücksichtigt bleiben, sodass sich als Prognosefunktion für das einfache exponentielle Glätten

ergibt. Dieses Prognoseverfahren wird auch als adaptives Prognoseverfahren bezeichnet. Die Herausforderung dieser Methodik ist das Schätzen des Glättungsfaktor gw!

Schauen Sie sich die Hilfe in R zur Funktion HoltWinters() an, hier wird die Umsetzung des exponentiellen  Glättens beschrieben! Im ersten Anwendungsbeispiel führen wir nur das exponentielle Glätten durch:

> HW <- HoltWinters(GetreideIndex, alpha = 1, beta = F, gamma = F)
> plot(HW)
> plot(HW, main = "Holt-Winters-Glättung", sub = "Exponetielles Glätten: alpha = 1")

Das Argument alpha legt nur den exponentiellen Glättungsfaktor fest. Die R-Ausgabe des Modells macht uns darauf aufmerksam (without trend) und deswegen wird auch nur der letzte Wert a des Modells, vergleichbar mit dem Prognosestartwert oder dem Schnittpunkt der Y-Achse, ausgegeben (Coefficients):

Holt-Winters exponential smoothing without trend and without seasonal component.

Call:
 HoltWinters(x = GetreideIndex, alpha = 1, beta = F, gamma = F)

Smoothing parameters:
 alpha: 1
 beta : FALSE
 gamma: FALSE

Coefficients:
  [,1]
a 187.2

Eine Veränderung von alpha macht sich durch einen deutlich glatteren Darstellungsverlauf bemerkbar (roter Kurvenverlauf):

> HW <- HoltWinters(GetreideIndex, alpha = 0.5, beta = F, gamma = F)
> plot(HW, main = "Holt-Winters-Glättung", sub = "Exponetielles Glätten: alpha = 0,5")

Das dazugehörige Modell:

Holt-Winters exponential smoothing without trend and without seasonal component.

Call:
 HoltWinters(x = GetreideIndex, alpha = 0.5, beta = F, gamma = F)

Smoothing parameters:
 alpha: 0.5
 beta : FALSE
 gamma: FALSE

Coefficients:
   [,1]
a 179.2353

Wie gehen von der Annahme aus, dass unsere Getreide-Zeitreihe frei von einer Session-Komponente ist (gamma = F). Deswegen werden wir das Holt-Winters-Modell nur auf die Trendkomponente b über die Ausprägung des Arguments beta ausdehnen:

> HW <- HoltWinters(GetreideIndex, alpha = 0.5, beta = 0.5, gamma = F)

Holt-Winters exponential smoothing with trend and without seasonal component.

Call:
 HoltWinters(x = GetreideIndex, alpha = 0.5, beta = 0.5, gamma = F)

Smoothing parameters:
 alpha: 0.5
 beta : 0.5
 gamma: FALSE

Coefficients:
    [,1]
a 188.461868
b  7.734264

Die Anpassung des Modells (roter Kurvenverlauf) im Vergleich zu der beobachteten Zeitreihe (schwarzer Kurvenverlauf) sieht wie folgt aus:

Das Vorhersagemodell Holt-Winters ohne Saison-Komponente ist folgende lineare Funktion (die Glättung wurde berücksichtigt, siehe R-Hilfe-Funktion zu Holt-Winters):

Y t+1 = a + b * t
t: Zeitschritt

Machen wir über die predict()-Funktion eine Vorhersage auf Basis obigem Modell für 15 Zeitschritte (Argument n.ahead = 15, Ausgabe des Vertrauensbereiches für P = 95%) in die Zukunft und lassen uns diese grafisch anzeigen:

> HW_P <- predict(HW, 15, prediction.interval = T)
> plot(HW, HW_P, main = "Vorhersage, Holt-Winters-Glättung", sub = "Exponetielles Glätten mit alpha = 0,5 und Trend (beta = 0,5)")

Ab der gestrichelten Linie in obiger Grafik wird die Prognose (rot) und der Vorhersage-Intervall (blau) angezeigt. Die Güte der Prognose ist von der Wahl der Parameter alpha für die Glättung des Niveaus (Parameter a, Startpunkt der Prognose) und beta (Parameter b, Glättung des linearen Trends) abhängig.

Die Schätzung der Prognose lässt sich natürlich auch numerisch ausgegeben (fit: Prognose, upr/lwr: Vorhersage-Intervall):

> HW_P

       fit   upr   lwr
Jan 2011 196.1961 206.7211 185.6712
Feb 2011 203.9304 217.0866 190.7742
Mar 2011 211.6647 228.5128 194.8165
Apr 2011 219.3989 240.7752 198.0226
May 2011 227.1332 253.7074 200.5589
Jun 2011 234.8675 267.2007 202.5342
Jul 2011 242.6017 281.1833 204.0201
Aug 2011 250.3360 295.6056 205.0664
Sep 2011 258.0702 310.4313 205.7092
Oct 2011 265.8045 325.6327 205.9763
Nov 2011 273.5388 341.1878 205.8897
Dec 2011 281.2730 357.0784 205.4677
Jan 2012 289.0073 373.2892 204.7254
Feb 2012 296.7416 389.8072 203.6759
Mar 2012 304.4758 406.6209 202.3308

Hier die grafische Darstellung eines Beispiels für beta = 0,2...

… auf Basis dieser Modell-Parameter:

Holt-Winters exponential smoothing with trend and without seasonal component.

Call:
 HoltWinters(x = GetreideIndex, alpha = 0.5, beta = 0.2, gamma = F)

Smoothing parameters:
 alpha: 0.5
 beta : 0.2
 gamma: FALSE

Coefficients:
    [,1]
a 185.722761
b  7.268521

Und hier die Vorhersage-Werte mit Vorhersage-Intervall:

> HW_P

       fit   upr   lwr
Jan 2011 192.9913 205.3149 180.6677
Feb 2011 200.2598 214.6314 185.8882
Mar 2011 207.5283 224.2902 190.7664
Apr 2011 214.7968 234.2431 195.3506
May 2011 222.0654 244.4523 199.6785
Jun 2011 229.3339 254.8886 203.7792
Jul 2011 236.6024 265.5300 207.6748
Aug 2011 243.8709 276.3594 211.3824
Sep 2011 251.1394 287.3632 214.9157
Oct 2011 258.4080 298.5306 218.2853
Nov 2011 265.6765 309.8527 221.5003
Dec 2011 272.9450 321.3219 224.5681
Jan 2012 280.2135 332.9320 227.4951
Feb 2012 287.4821 344.6773 230.2868
Mar 2012 294.7506 356.5531 232.9481

Der Einfluss, insbesondere des beta-Parameters, ist deutlich zu erkennen. Visuell und durch die dargestellten numerischen Vorhersagedaten wird der Eindruck erweckt, die Güte der Vorhersage ist durch die Wahl beta = 0,2  im Vergleich zu beta = 0,5 erheblich gesteigert worden.

Wird fortgesetzt!

 

 

Seitenanfang

Hat der Inhalt Ihnen weitergeholfen und Sie möchten diese Seiten unterstützen?