Ekstrapolacja szeregów czasowych
Marcin Orchel
1 Wstęp
Sezonowość to zależność korelacyjna rzędu k między i-tym elementem szeregu a (i − k)- tym elementem. Miarą sezonowości jest autokorelacja.
Usunięcie składowej k poprzez różnicowanie szeregu – przekształcenie każdego i-tego elementu szeregu na jego różnicę z (i − k)-tym elementem.
Szeregi stacjonarne na potrzeby ARIMA, czyli szeregi o stałych w czasie średniej, wariancji i autokorelacji.
2 Zadania
2.1 Zadania na 3.0
Napisać skrypt w R. W skrypcie
• Przetestować wygładzanie szeregu za pomocą średniej ruchomej
• Przetestować transformację Boxa-Coxa logarytmiczną (λ = 0) i pierwiastkową (λ = 0, 5)
• Przetestować różną szerokość okna wygładzania oraz różne metody: simple, Trian- gular, Exponential Simple, Exponential Modified, Cumulative.
• Dokonać ekstrapolacji modelem ARIMA oraz modelem lm na wygenerowanym szeregu czasowym z sezonowością.
• Porównać jakość ekstrapolacji dla różnych modeli.
• Przetestować różne wartości parametrów tych modeli.
• Narysować na jednym wykresie punkty treningowe i funkcje ekstrapolacji dla wy- branych modeli.
• Sprawdzić czy zlogarytmowanie obserwacji pozwala na usunięcie trendu.
• Sprawdzić czy gdy w szeregu czasowym znajdują się obserwacje odstające, lepiej użyć mediany niż średniej kroczącej.
• Dodać komentarz do skryptu opisujący krótko na czym polegają użyte metody oraz wnioski z badań.
Wskazówki:
• opcjonalnie można przetestować najpierw model średniej kroczącej ma https://
www.rdocumentation.org/packages/forecast/topics/ma
• można sprawdzić dekompozycję szeregu na różne składniki takie jak sezonowość po- leceniem decomposehttps://www.rdocumentation.org/packages/stats/topics/
decompose
• jeśli szereg czasowy jest stacjonarny to nie ma potrzeby różnicowania, p-value to graniczny poziom istotności. Im mniejsza tym bardziej prawdopodobna jest hipoteza alternatywna, że szereg czasowy jest stacjonarny
• transformacja Boxa-Coxa przekształcenie potęgowe stosujemy gdy wariancja wa- riancja wzrasta lub maleje, gdy w szeregu występują obserwacje odstające, oraz aby uzyskać rozkład danych zbliżony do normalnego
• różnicowanie w celu przekształcenia do postaci stacjonarnej (np. wyeliminowanie trendu długoterminowego, sezonowości)
Wskazówki do R:
• przykładowe szeregi czasowe można znaleźć w pakiecie forecast lub standardowo w pakiecie MASS i datasets jak w konspekcie Wprowadzenie do R
• sprawdzenie czy szereg czasowy jest stacjonarny, funkcja Box.testhttps://stat.
ethz.ch/R-manual/R-devel/library/stats/html/box.test.htmllub funkcje adf.test https://www.rdocumentation.org/packages/aTSA/topics/adf.testi kpss.test https://www.rdocumentation.org/packages/tseries/topics/kpss.test
• https://www.rdocumentation.org/packages/forecast/topics/forecast
• https://stat.ethz.ch/R-manual/R-devel/library/stats/html/arima.html
• https://www.rdocumentation.org/packages/forecast/topics/ets
• https://www.rdocumentation.org/packages/forecast/topics/auto.arima
• https://www.rdocumentation.org/packages/forecast/topics/tslm
• https://www.rdocumentation.org/packages/forecast/topics/forecast.lm
• https://www.rdocumentation.org/packages/forecast/topics/BoxCox
• https://www.rdocumentation.org/packages/xts/topics/xts
• https://stat.ethz.ch/R-manual/R-devel/library/stats/html/ts.html
• https://stat.ethz.ch/R-manual/R-devel/library/stats/html/tsdiag.html
• http://search.r-project.org/library/zoo/html/read.zoo.html
• https://stat.ethz.ch/R-manual/R-devel/library/stats/html/stl.html
• https://stat.ethz.ch/R-manual/R-devel/library/stats/html/lag.plot.html
• https://stat.ethz.ch/R-manual/R-devel/library/stats/html/spectrum.html
• https://stat.ethz.ch/R-manual/R-devel/library/stats/html/box.test.html
• https://stat.ethz.ch/R-manual/R-devel/library/stats/html/HoltWinters.
htmlorazhttps://www.rdocumentation.org/packages/forecast/topics/forecast.
HoltWinters
• https://www.rdocumentation.org/packages/SparseM/topics/slm.fit
• https://www.rdocumentation.org/packages/tframe/topics/tfwindow
• https://www.rdocumentation.org/packages/TTR/topics/GD
• https://www.rdocumentation.org/packages/forecast/topics/wineind
• https://www.rdocumentation.org/packages/forecast/topics/accuracy
• https://www.rdocumentation.org/packages/smatr/topics/sma Przykłady
• x = seq(from=1, to=50, by=0.1)
y = x + sin(x) + rnorm(length(x), 0, 1) zbior = cbind(x, y)
plot(x,y)
zbior.ts <- ts(y, start=c(1), frequency=length(x)/8.0) zbior.ts.dekompozycja = decompose(zbior.ts, "additive") plot(zbior.ts.dekompozycja)
• informacje o szeregu czasowym library(forecast)
zbior.ts <- ts(zbior, start=c(2009, 1), frequency=12) start(zbior.ts)
end(zbior.ts)
frequency(zbior.ts) zbior.ts
plot(zbior.ts) par(mfrow=c(2,1))
monthplot(zbior.ts, main="Miesieczny")
seasonplot(zbior.ts, year.labels=TRUE, col=rainbow(5), main="Sezonowy")
• wykres funkcji autokorelacji Acf(zbior.ts)
• dekompozycja szeregu na trend, sezonowość i losowe fluktuacje
zbior.ts.dekompozycja.add <- decompose(zbior.ts, "additive") plot(zbior.ts.dekompozycja.add)
• Transformacja Boxa-Coxa par(mfrow=c(3,1))
plot(zbior.ts, main="Dane", xlab="", ylab="") grid()
plot(BoxCox(zbior.ts, lambda=0.5), main="Box-Cox 1", xlab="", ylab="")
grid()
plot(BoxCox(zbior.ts, lambda=0), main="Box-Cox 2", xlab="", ylab="")
grid()
• logarytmowanie danych
tsZbior <- ts(zbior, frequency=12, start=c(1981, 1)) logTsZbior <- log(tsZbior)
par(mfrow=c(1,2))
plot(tsZbior, lwd=4, xlab="Years", ylab="Y1", main="Main") plot(logTsZbior, lwd=4, xlab="Years", ylab="Y1", main="Main po
zlogarytmizowaniu")
• wygładzanie library(TTR)
tsZbiorSmooth <- SMA(zbior.ts, n=5) par(mfrow=c(2,1))
plot(zbior.ts, lwd=4, xlab="Indeks", ylab="Y1", main="Main") plot(tsZbiorSmooth, lwd=4, xlab="Indeks", ylab="Y1",
main="Wygladzony")
• wygładzanie medianą par(mfrow = c(3, 1))
tsZbiorSMooth3 <- smooth(zbior) tsZbiorSmooth5 <- runmed(zbior, 5)
plot(tsZbior, lwd=4, xlab="Indeks", ylab="Y1", main="Main") plot.ts(tsZbiorSmooth3, lwd=4, xlab="Indeks", ylab="Y1",
main="Wygladzony 1")
plot.ts(tsZbiorSmooth5, lwd=4, xlab="Indeks", ylab="Y1", main="Wygladzony 2")
• różnicowanie szeregów czasowych diffTsZbior <- diff(zbior.ts)
plot(zbior.ts, lwd=4, xlab="Years", ylab="Y1", main="Main") plot(diffTsZbior, lwd=4, xlab="Years", ylab="Y1", main="Main")
• różnicowanie szeregów czasowych 2 tsdisplay(diff(zbior.ts, 12))
• sprawdzenie stacjonarności szeregu czasowego library(tseries)
Box.test(zbior, lag=10, type="Box-Pierce") Box.test(zbior, lag=10, type="Ljung-Box") adf.test(zbior, alternative="stationary") kpss.test(zbior)
• podział na zbiór treningowy i testowy
zbior.learn <- window(zbior.ts, end=c(6)) start(zbior.learn)
end(zbior.learn)
zbior.test <- window(zbior.ts, start=c(6.001)) start(zbior.test)
end(zbior.test) length(zbior.test)
ts.plot(zbior.learn, zbior.test, col=1:2, lty=c(1,2)) abline(v=time(zbior.test)[1], lty=3)
legend("topleft", legend=c("training", "test"), col=1:2, lty=c(1,2))
• model ARIMA
model.ARIMA1 <- auto.arima(zbior.learn) summary(model.ARIMA1)
tsdiag(model.ARIMA1)
• ręczny wybór modelu
tsdisplay(diff(zbior.learn, 12), lag.max=12)
model.ARIMA2 <- Arima(zbior.learn, order=c(1,0,0), season=c(0,1,0)) summary(model.ARIMA2)
tsdiag(model.ARIMA2)
• model dekompozycji
model.tslm <- tslm(zbior.learn ~ trend + season) summary(model.tslm)
Acf(residuals(model.tslm))
Box.test(residuals(model.tslm), type="Ljung-Box", lag=20)
• model ETS (Exponential Smoothing) model.ets <- ets(zbior.learn) summary(model.ets)
tsdiag(model.ets)
• prognozowanie
horyzont <- length(zbior.test)
snaive.prognozy <- snaive(zbior.learn, h=horyzont) ARIMA1.prognozy <- forecast(model.ARIMA1, h=horyzont) ARIMA2.prognozy <- forecast(model.ARIMA2, h=horyzont) tslm.prognozy <- forecast(model.tslm, h=horyzont) ets.prognozy <- forecast(model.ets, h=horyzont) par(mfrow=c(2,2))
plot(snaive.prognozy)
lines(zbior.test, col="red", lty=2)
legend("topleft", legend=c("historyczne", "prognozy"), lty=c(2,1), col=c("red", "blue"))
grid()
plot(ARIMA2.prognozy)
lines(zbior.test, col="red", lty=2) grid()
plot(tslm.prognozy)
lines(zbior.test, col="red", lty=2) grid()
plot(ets.prognozy)
lines(zbior.test, col="red", lty=2) grid()
• porównanie jakości prognoz
kryteria <- c("RMSE", "MAE", "MAPE", "MASE") accuracy(snaive.prognozy, zbior.test)[, kryteria]
accuracy(ARIMA1.prognozy, zbior.test)[, kryteria]
accuracy(ARIMA2.prognozy, zbior.test)[, kryteria]
accuracy(tslm.prognozy, zbior.test)[, kryteria]
accuracy(ets.prognozy, zbior.test)[, kryteria]
