08.04

2016

Prognozowanie szeregów czasowych o złożonej sezonowości na przykładzie danych dotyczących sprzedaży benzyny w USA

Autor: adam

Wiele spotykanych w praktyce szeregów czasowych cechuje występowanie niestandardowych wzorców sezonowych (sezonowości). Typowy przykład takiej sytuacji to szeregi o ułamkowym (niecałkowitym) okresie, np. szeregi zawierające tygodniowe wielkości produkcji lub wielkości sprzedaży wykazują często sezonowość roczną o okresie $365.25/7 \approx 52.18$ .

W takim przypadku, klasyczne modele wykorzystywane do analizy i prognozowania szeregów czasowych, które nie uwzględniają poprawnie specyficznych własności danych, mogą nie być efektywne. Pojawia się więc konieczność stosowania metod dedykowanych szeregom o niestandardowej sezonowości.

Dzisiejszy wpis przedstawia interesujące case study przygotowane przez Piotra Kopszaka, dotyczące prognozowania wielkości sprzedaży benzyny w USA (dane tygodniowe). W analizie szczegółowo porównana została m.in. skuteczność klasycznych oraz nowszych modeli, opracowanych z myślą o występującej w danych nietypowej sezonowości.

Na koniec zaznaczmy jeszcze, że poniższy artykuł (z założenia) nie jest tutorialem. który ma nauczyć od podstaw stosowania metod analizy i prognozowania szeregów czasowych. Osoby zainteresowane poznaniem podstaw metodologicznych, szczegółów dotyczących określonych metod analizy szeregów czasowych czy też konkretnych algorytmów dostępnych w pakiecie R odsyłamy do naszej książki: Analiza i prognozowanie szeregów czasowych. Praktyczne wprowadzenie na podstawie środowiska R, A. Zagdański, A. Suchwałko, PWN 2015.

1 Opis danych

Przeprowadzimy analizę danych tygodniowych dotyczących wielkości sprzedaży benzyny w USA (w baryłkach) od lutego 1991 do lipca 2005, łącznie 744 wartości. Dane pochodzą z raportów EIA (U.S. Energy Information Administration), a w wersji wykorzystanej w analizie zostały pobrane ze strony prof. R.Hyndmana http://robjhyndman.com/papers/complex-seasonality/. Analizowany szereg czasowy przedstawiono na wykresie 1.

Dane te są interesujące pod kątem badania złożonej sezonowości, ponieważ występuje w nich ułamkowy okres, związany z niecałkowitą liczbą tygodni w roku (długość okresu wynosi $365.25/7=$ 52.1785714). Taka specyfika danych wymaga stosowania adekwatnych metod konstrukcji prognoz, opracowanych z myślą o, występujących w analizowanym szeregu czasowym, niestandardowych wzorcach sezonowych (sezonowości).

2 Cel analizy

Porównamy prognozy otrzymane za pomocą kilku konkurencyjnych metod/modeli, w szczególności w analizie uwzględnione będą:

  • regresja harmoniczna (w kilku wariantach) – modelowanie złożonej sezonowości z wykorzystaniem narzędzi analizy fourierowskiej (analizy spektralnej),
  • BATS i TBATS – modele wygładzania wykładniczego dedykowane szeregom czasowym o złożonej sezonowości, uwzględniające: transformację Boxa-Coxa, zakłócenia ARMA oraz odpowiednie modelowanie składowej trendu (Trend) i sezonowości (Seasonality),
  • sARIMA – klasyczny sezonowy model ARIMA,
  • STL+ETS – połączenie dekompozycji STL (do wyznaczenia składowej sezonowej) oraz modelu ETS (Error-Trend-Seasonal) do konstrukcji prognozy na podstawie danych odsezonowanych.

Modele sARIMA oraz STL+ETS potraktowane zostaną jako metody referencyjne, ponieważ nie pozwalają one na uwzględnienie różnych aspektów złożonej sezonowości, w tym niecałkowitego okresu, który występuje w analizowanym szeregu czasowym.

Zbadamy, jak zmienia się dokładność prognoz w zależności od długości horyzontu prognozy oraz porównamy szerokości skonstruowanych przedziałów predykcyjnych (ang. prediction intervals). Sprawdzimy także poprawność dopasowania poszczególnych modeli do danych, na bazie których były zbudowane.

Wszystkie analizy zostały przeprowadzone w pakiecie $R$ w wersji 3.2.4.

Rysunek 1. Sprzedaż benzyny w USA (w tysiącach baryłek na dzień). Dane tygodniowe z okresu od lutego 1991 do lipca 2005.

Rysunek 1. Sprzedaż benzyny w USA (w tysiącach baryłek na dzień). Dane tygodniowe z okresu od lutego 1991 do lipca 2005.

3 Wstępna analiza

Jak widać na wykresie 1 w danych występuje wyraźny liniowy trend wzrostowy. W pierwszej kolejności usuniemy go z danych, aby móc dokładniej zbadać występujące efekty sezonowe. Do estymacji trendu wykorzystamy klasyczną metodę najmniejszych kwadratów.

Sprawdzimy teraz, wykorzystując test g-Fishera, czy w danych, z których wyeliminowaliśmy trend długoterminowy, obecny jest jakikolwiek komponent okresowy. Hipoteza zerowa zakłada, że dane są białym szumem, natomiast według hipotezy alternatywnej w danych występuje komponent okresowy o nieokreślonej częstotliwości. Wyniki testu przedstawiamy w tabeli 1:


Tabela 1. Wyniki testu g-Fishera na obecność składnika okresowego o nieznanej częstotliwości.
p-wartość statystyka testowa $\xi$
0.000 75.425

Jak widać, na poziomie istotności $\alpha = 0.05$ odrzucamy hipotezę zerową o braku okresowości w danych.

Rysunek 2. Periodogram dla danych po usunięciu liniowego trendu

Rysunek 2. Periodogram dla danych po usunięciu liniowego trendu

Na wykresie 2 przedstawiamy periodogram analizowanych danych. Widoczne są dwa największe piki odpowiadające częstościom $\nu_{14}=\omega_{14} / {2\pi}=14/744=$ 0.019, $\nu_{15}=\omega_{15}/{2\pi}=15/744=$ 0.02. Wykorzystując związek pomiędzy częstotliwością $\omega$ i okresem $T$, tzn. $\omega=2\pi/T$, otrzymujemy wartości $T_{14}\approx$ 53.57, $T_{15}\approx$ 50, gdzie $T_j$ jet okresem odpowiadającym $j$-tej częstości fourierowskiej. Jest to wyraźny argument świadczący o występowaniu sezonowości o okresie zbliżonym do średniej liczby tygodni w roku, która wynosi 52.178. Dodatkowo, na wykresie periodogramu możemy zaobserwować tzw. ,,przeciek” widmowy (ang. spectral leakage) – zjawisko polegające na przypisaniu wysokich wartości $I(\omega)$ kilku częstościom bliskim rzeczywistej częstości (odpowiadającej okresowości rocznej). W tym przypadku, polega ono na występowaniu obok siebie pików dla $\omega_{14}$ oraz $\omega_{15}$.

Możemy również zauważyć (Rys.2) duże wartości periodogramu dla częstości $\omega_{43}$ oraz $\omega_{86}$. Pierwszej z nich odpowiada okres o długości około 17 tygodni, natomiast drugiej okres o długości około 8 tygodni. Może to być odzwierciedlenie pewnych wahań sezonowych występujących w branży paliwowej, o okresach krótszych niż rok. Wykres 3 przedstawia autokorelogram analizowanych danych. Również on wyraźnie wskazuje na występowanie sezonowości rocznej, o okresie zbliżonym do 52.

Rysunek 3. Autokorelogram dla danych po usunięciu liniowego trendu

Rysunek 3. Autokorelogram dla danych po usunięciu liniowego trendu

4 Porównanie modeli

W dalszej części analizy porównamy dokładność prognoz skonstruowanych na podstawie różnych metod/modeli, uwzględniających regularności występujące w naszych danych. W celu oceny dokładności prognoz dokonujemy podziału na zbiór uczący (ang. training set) o długości 594 (a więc odpowiadający okresowi od początku 1991 r. do połowy 2002 r.), na podstawie którego będziemy budowali nasze modele oraz zbiór testowy (ang. test set) o długości 150, który wykorzystamy do porównania dokładności skonstruowanych prognoz.

Na początek zbadamy jakość dopasowania poszczególnych modeli do zbioru uczącego, a następnie przyjrzymy się prognozom otrzymanym na bazie kolejnych modeli, porównując w szczególności miary błędów predykcji oraz szerokości przedziałów predykcyjnych. Sprawdzimy również wykorzystując test Diebolda-Mariano istotność różnic dokładności prognoz dla poszczególnych metod. Ponieważ każda z rozważanych metod uwzględnia możliwość wystąpienia trendu, modele będą dopasowane do wyjściowych danych, a więc zawierających także trend długoterminowy (trend liniowy).

4.1 Regresja harmoniczna

4.1.1 Zakłócenia nieskorelowane (biały szum)

Na początek dopasujemy do naszych danych model regresji harmonicznej z automatycznie zidentyfikowanymi (na podstawie testu statystycznego) istotnymi składowymi okresowymi (harmonikami), przyjmując także najprostszy możliwy model dla zakłócenia losowego, tzn. biały szum (ang. white noise).

W tabeli 2 przedstawiamy statystycznie istotne częstości, które zostały wybrane do modelu.


Tabela 2. Statystycznie istotne częstości fourierowskie
$\omega_{14}$ 0.12
$\omega_{15}$ 0.13
$\omega_{43}$ 0.36

Dopasowany model zawiera więc trzy pary funkcji sinus i cosinus, tzn. ma postać $$X_t = \sum_{i=1}^3 A_{k_i}cos(\omega_{k_i}t) + B_{k_i}sin(\omega_{k_i}t) + Z_t,$$ gdzie $A_{k_i}, B_{k_i}$ są estymowanymi współczynnikami odpowiadającymi harmonikom o częstościach $\omega_{k_i}$, $k_i \in\{14,15,43\}$$Z_t$ jest gaussowskim białym szumem. Na wykresie 4 przedstawiamy dopasowanie tego modelu do danych.

Rysunek 4. Dopasowanie modelu regresji harmonicznej z nieskorelowanymi zakłóceniami.

Rysunek 4. Dopasowanie modelu regresji harmonicznej z nieskorelowanymi zakłóceniami.

Zastosowanie kilku harmonik sprawia, że dopasowany model wykazuje tendencję do zmiany amplitudy w czasie, co niekoniecznie jest pożądanym zjawiskiem (widać to zwłaszcza jeśli porówna się dopasowanie do danych dla ostatniego okresu). W tabeli 3 przedstawiamy miary dopasowania modelu do zbioru uczącego.


Tabela 3. Miary błędów dopasowania modelu regresji harmonicznej do zbioru uczącego
RMSE MAE MAPE MASE ACF1
298.18 238.34 3.04 0.70 0.02

Kryteria wskazują na względnie dobre dopasowanie modelu do danych. Średni bezwzględny procentowy błąd (MAPE) jest na poziomie 3%, a prognozy wewnątrz zbioru uczącego są dokładniejsze niż średnia prognoza naiwna (błąd MASE<1).

4.1.2 Zakłócenia typu ARMA

Dla porównania dopasujemy teraz inny, ogólniejszy wariant modelu regresji harmoniczny, zakładający, że szereg błędów (zakłóceń) jest stacjonarnym modelem typu ARMA (dopuszczamy więc obecność krótkoterminowej korelacji czasowej w analizowanym szeregu). Na wykresie 5 znajduje się wykres dopasowania modelu do danych. Dla reszt dopasowany został model ARMA(1,2) (rzędy: $p=1, q=2$).

Rysunek 5. Dopasowanie modelu regresji harmonicznej z zakłóceniami (błędami) typu ARMA.

Rysunek 5. Dopasowanie modelu regresji harmonicznej z zakłóceniami (błędami) typu ARMA.

Na wykresie 6 przedstawiamy reszty z modelu regresji harmonicznej (zakładającego białoszumowość reszt) oraz ich autokorelogramy (ACF i PACF).

Rysunek 6. Wykres reszt pochodzących z modelu regresji harmonicznej (reszty nieskorelowane).

Rysunek 6. Wykres reszt pochodzących z modelu regresji harmonicznej (reszty nieskorelowane).

Na wykresie 7 przedstawiamy reszty z modelu regresji harmonicznej, w przypadku którego zakładamy, że szereg reszt (zakłóceń) jest szeregiem typu ARMA.

Rysunek 7. Wykres i autokorelogram reszt z modelu regresji harmonicznej z zakłóceniami typu ARMA.

Rysunek 7. Wykres i autokorelogram reszt z modelu regresji harmonicznej z zakłóceniami typu ARMA.

Wykres 6 pokazuje, że nie udało się uzyskać końcowych reszt będących białym szumem: na wykresie funkcji autokorelacji (ACF) widoczne są istotne autokorelacje, świadczące o tym, że założenie o białoszumowości reszt nie jest spełnione. Z tego powodu warto rozważyć ogólniejszy model ARMA (uwzględniający obecność autokorelacji w szeregu reszt/zakłóceń). W wypadku reszt przedstawionych na wykresie 7są one bardziej zbliżone do białego szumu (szeregu nieskorelowanego): liczba statystycznie istotnych autokorelacji próbkowych zmniejszyła się w porównaniu do Rys.6.

Jak widać w tabeli 4, przyjęcie ogólniejszego założenia o postaci reszt wpłynęło także na nieznaczną poprawę jakości dopasowania modelu do danych (w porównaniu do tabeli 3).


Tabela 4. Miary błędów dopasowania modelu regresji harmonicznej do zbioru uczącego
RMSE MAE MAPE MASE ACF1
295.09 236.43 3.01 0.70 -0.00

4.1.3 Sezonowość roczna i zakłócenia typu ARMA

Rozważając różne warianty regresji harmonicznej, zbadamy jeszcze model, w którym występuje jedynie główna składowa okresowa, tzn. składowa harmoniczna o okresie odpowiadającym sezonowości rocznej równym $T=365.25/7=$ 52.1785714. Model ten przyjmuje zatem postać $$X_t = A\cos(\omega t) + B\sin(\omega t) + Z_t,$$ gdzie $\omega=2\pi/T$.

Dopasowanie modelu do danych prezentujemy na wykresie 8, natomiast na wykresie 9 przedstawione są reszty z modelu.

Rysunek 8. Dopasowanie modelu regresji harmonicznej z sezonowością roczną oraz błędami (zakłóceniami) ARMA.

Rysunek 8. Dopasowanie modelu regresji harmonicznej z sezonowością roczną oraz błędami (zakłóceniami) ARMA.

Rysunek 9. Wykres i autokorelogram reszt z modelu zakładającego sezonowość roczną oraz postać ARMA zakłócenia losowego.

Rysunek 9. Wykres i autokorelogram reszt z modelu zakładającego sezonowość roczną oraz postać ARMA zakłócenia losowego.

W tabeli 5 przedstawiamy wartości kryteriów oceniających jakość dopasowania modelu.


Tabela 5. Miary błędów dopasowania modelu regresji harmonicznej do zbioru uczącego
RMSE MAE MAPE MASE ACF1
289.87 233.16 2.97 0.69 -0.01

Zastosowane podejście nie poprawiło losowości reszt (widoczna istotna autokorelacja na Rys.9), natomiast w przeciwieństwie do poprzednich modeli uzyskaliśmy stałą amplitudę składowej sezonowej. Jakość dopasowania modelu (tab. 5) jest jednak nieznacznie lepsza w porównaniu do poprzednich wariantów modelu regresji harmonicznej.

4.2 Model TBATS

Parametry dopasowanego modelu TBATS przedstawiono w tabeli 6.


Tabela 6. Parametry dopasowanego modelu TBATS
$\alpha$ $\beta$ $\phi$ $\gamma_1^{(1)}$ $\gamma_2^{(1)}$ p q
0.05 0.00065 1.00 -0.0042322 0.00121 0 0

Parametry $\alpha$$\beta$ są parametrami wygładzającym, decydującym o tym jak wraz z czasem zmienia się wpływ obserwacji historycznych na wartość poziomu i trendu. $\phi$ oznacza parametr tłumienia trendu (ang. damping parameter). Parametry $\gamma_1^{(1)}$ oraz $\gamma_2^{(1)}$ są parametrami wygładzającymi w równaniach rekurencyjnych opisujących dynamikę zmian składowej sezonowej. Wartości $p=0, q=0$ oznaczają natomiast, że losowe zakłócenia modelowane są jako biały szum. Na wykresie 10 prezentujemy dopasowanie modelu do danych, a na wykresie 11 przedstawiamy reszty z modelu.

Rysunek 10. Dopasowanie modelu TBATS do zbioru uczącego.

Rysunek 10. Dopasowanie modelu TBATS do zbioru uczącego.

Rysunek 11. Reszty z modelu TBATS.

Rysunek 11. Reszty z modelu TBATS.

W tabeli 7 przedstawiamy wartości błędów dopasowania modelu TBATS do zbioru uczącego.


Tabela 7. Miary błędów dopasowania modelu TBATS do zbioru uczącego
RMSE MAE MAPE MASE ACF1
264.29 207.92 2.64 0.61 -0.18

Zauważmy, że wartości błędów dopasowania są niższe niż w przypadku regresji harmonicznej. Dokładność dopasowania modelu TBATS jest zatem lepsza, choć wykres reszt (Rys.11), podobnie jak w przypadku regresji harmonicznej, świadczy o pewnej niedoskonałości tego modelu (widoczne pojedyncze istotne korelacje).

4.3 Model BATS

Przeanalizujemy teraz dopasowanie modelu BATS do danych. Ponieważ model BATS pozwala na uwzględnienie jedynie sezonowości o okresie całkowitym (a nie ułamkowym), przyjmujemy że w naszych danych występuje składowa sezonowa okres o długości 52. Zbudowany model ma parametry przedstawione w tabeli 8.


Tabela 8. Parametry w dopasowanym modelu BATS
$\alpha$ $\beta$ $\gamma$ $\phi$ p q
0.089429 -0.000022 -0.122968 1.0 0 0

Parametr $\alpha$ jest parametrem wygładzającym dla poziomu a $\beta$ parametrem wygładzającym dla trendu. W tym przypadku nie stosujemy tłumienia trendu ($\phi=1$). Parametr $\gamma$ jest parametrem wygładzania w równaniu opisującym zmienność sezonową (ponieważ jest ono mniej skomplikowane niż w przypadku modelu TBATS, występuje tylko jeden parametr). Podobnie jak w modelu TBATS, w przypadku modelu BATS reszty są modelowane jako biały szum. Na wykresie 12 przedstawiamy dopasowanie modelu do zbioru uczącego, natomiast na wykresie 13 przedstawiamy reszty (błędy dopasowania na zbiorze uczącym).

Rysunek 12. Dopasowanie modelu BATS do zbioru uczącego.

Rysunek 12. Dopasowanie modelu BATS do zbioru uczącego.

Rysunek 13. Reszty z modelu BATS.

Rysunek 13. Reszty z modelu BATS.

W tabeli 9 przedstawiamy miary błędów dopasowania modelu BATS do zbioru uczącego.


Tabela 9. Miary błędów dopasowania modelu BATS do zbioru uczącego
RMSE MAE MAPE MASE ACF1
259.93 205.24 2.62 0.60 -0.12

Możemy zauważyć, że dokładność dopasowania modelu BATS jest porównywalna jak w przypadku modelu TBATS i lepsza niż w przypadku wszystkich wariantów regresji harmonicznej. W tym przypadku również możemy mieć pewne zastrzeżenia co do braku korelacji końcowych reszt z modelu, aczkolwiek zachowanie reszt jest bliższe białoszumowości niż w przypadku poprzednich modeli.

4.4 Model sezonowy ARIMA (sARIMA)

Do naszych danych dopasujemy teraz klasyczny model sARIMA (sezonowy model ARIMA). W przypadku tego modelu również zaokrąglamy okres składowej sezonowej do wartości 52. Model sARIMA nie może być oczywiście traktowany jako metoda dedykowana prognozowaniu szeregów czasowych o złożonej sezonowości i wykorzystujemy ją głównie jako metodę referencyjną. Jako optymalny wybrany został model $\mathrm{sARIMA}(0,1,1)(0,0,2)_{52}$ z dryfem. Jego parametry przedstawiamy w tabeli 10


Tabela 10. Współczynniki w dopasowanym modelu sARIMA
$\theta_1$ $\Theta_1$ $\Theta_2$ $c$
-0.81 0.19 0.13 3.29

Istotny jest fakt, że zmienność sezonowa jest w całości modelowana przez stacjonarną część modelu sARIMA (nie występuje różnicowanie sezonowe). Na wykresie 14 przedstawiamy wykres dopasowania modelu do zbioru uczącego, natomiast na wykresie 15 przedstawiamy reszty z dopasowanego modelu.

Rysunek 14. Dopasowanie modelu sARIMA do zbioru uczącego.

Rysunek 14. Dopasowanie modelu sARIMA do zbioru uczącego.

Rysunek 15. Reszty z modelu sARIMA.

Rysunek 15. Reszty z modelu sARIMA.

W tabeli 11 przedstawiamy miary dokładności dopasowania modelu sARIMA do zbioru uczącego.


Tabela 11. Miary błędów dopasowania modelu sARIMA do zbioru uczącego
RMSE MAE MAPE MASE ACF1
316.56 248.46 3.17 0.73 -0.07

Widzimy, że model sARIMA jest dużo gorzej dopasowany niż poprzednie (konkurencyjne) modele. Świadczy o tym m.in. autokorelogram reszt (wykres 15), wskazujący na wyraźną pozostałość wzorców okresowych. Ponadto miary dokładności dopasowania modelu świadczą o najniższej dokładności z dotychczas omawianych modeli.

4.5 Model oparty na dekompozycji (STL+ETS)

Ostatnią omawianą metodą, będzie połączenie dekompozycji STL oraz modelu ETS (w skrócie: STL+ETS). Z uwagi na wymagania modelu, również w tym wypadku stosujemy zaokrąglenie okresu sezonowości do wartości całkowitej równej 52. Model opisujący szereg po usunięciu składowej sezonowej to ETS(A,A,N) (a więc model zakładający występowanie addytywnej składowej losowej oraz addytywnego trendu). Na wykresie 16 przedstawiamy wykres dopasowania modelu STL+ETS do zbioru uczącego, natomiast na wykresie 17 przedstawiamy reszty z tego modelu.

Rysunek 16. Dopasowanie modelu STL+ETS do zbioru uczącego.

Rysunek 16. Dopasowanie modelu STL+ETS do zbioru uczącego.

Rysunek 17. Reszty z modelu STL+ETS.

Rysunek 17. Reszty z modelu STL+ETS.

W tabeli 12 przedstawiamy miary dokładności dopasowania modelu STL+ETS do zbioru uczącego.


Tabela 12. Miary błędów dopasowania modelu STL+ETS do zbioru uczącego
RMSE MAE MAPE MASE ACF1
265.88 211.60 2.69 0.62 -0.14

Jak widać, dokładność dopasowania jest jedną z najwyższych z dotychczas otrzymanych. Może wynikać to z faktu, że dekompozycja STL, jako metoda nieparametryczna, ma stosunkowo dużą zdolność do prawidłowego modelowania nieregularnego wzorca sezonowego (mimo występowania niecałkowitego okresu).

5 Porównanie dokładności prognoz

Porównamy teraz dokładności prognoz dla zbioru testowego, skonstruowanych na podstawie dopasowanych modeli. Przeanalizujemy m.in. jak zmienia się trafność prognoz w zależności od długości horyzontu (będziemy rozważali horyzont od 10 do 150 tygodni). Porównamy również otrzymane przedziały predykcyjne oraz złożoność obliczeniową (tzn. czas potrzebny na dopasowanie modeli do danych oraz wyznaczenie prognoz).

W porównaniu uwzględnione będą dwa warianty regresji harmonicznej. Przyjmujemy następujące oznaczenia:

  • regresja harmoniczna I — wybór składowych harmonicznych (istotnych częstości) oparty jest na odpowiednim teście statystycznym. Dodatkowo model zakłada brak korelacji (białoszumowość) zakłóceń.
  • regresja harmoniczna II — model z pojedynczą składową okresową, odpowiadającą sezonowości rocznej o okresie $365.25/7=$ 52.1785714 oraz losowymi zakłóceniami typu ARMA.

Prognozy przedstawimy najpierw na wspólnym wykresie (wykres 18).

Rysunek 18. Zestawienie prognoz skonstruowanych na bazie omawianych modeli.

Rysunek 18. Zestawienie prognoz skonstruowanych na bazie omawianych modeli.

Na wykresie 19 przedstawiamy zestawienie prognoz otrzymanych na bazie poszczególnych modeli, porównanych z wartościami ze zbioru testowego.

Rysunek 19. Porównanie prognoz otrzymanych na bazie dopasowanych modeli z rzeczywistymi wartościami ze zbioru testowego.

Rysunek 19. Porównanie prognoz otrzymanych na bazie dopasowanych modeli z rzeczywistymi wartościami ze zbioru testowego.

Jak widać, dopasowanie pojedynczej okresowości (regresja harmoniczna II) pozwoliło na lepsze odwzorowanie sezonowości rocznej. Gorsze dopasowanie modelu z wieloma harmonikami (regresja harmoniczna I) może być związane z przeuczeniem modelu: poza dwiema częstościami odpowiadającymi okresowi rocznemu ($\omega_{13}$ oraz $\omega_{14}$), do modelu została wybrana również częstość $\omega_{43}$ niemająca związku z okresem rocznym, a której obecność mogła wpłynąć negatywnie na dokładność prognoz. Ponadto, przyjęcie ogólniejszego modelu dla składnika losowego/zakłócenia (model ARMA zamiast białego szumu) nie miało dużego wpływu na postać prognozy, szczególnie w przypadku dłuższego horyzontu.

Analizując wyniki otrzymane na bazie modelu TBATS, możemy zauważyć, że model ten prowadzi do znacznie dokładniejszych prognoz w porównaniu z regresją harmoniczną. Co więcej, prognozy skonstruowane na bazie TBATS wydają się bardziej trafne (dla wszystkich rozważanych horyzontów) także w porównaniu do pozostałych metod/modeli.

Z kolei na wykresie przedstawiającym prognozy skonstruowane na bazie modelu BATS można zauważyć wpływ zaokrąglenia długości okresu do liczby całkowitej: szczególnie w dłuższym horyzoncie skrócenie długości okresu z 52.1785714 do 52 powoduje minimalne ,,wyprzedzanie” rzeczywistych wartości szeregu przez prognozy.

W przypadku prognoz skonstruowanych na bazie modelu sARIMA, największym mankamentem jest niezadowalające odzwierciedlenie w skonstruowanej prognozie zachowań sezonowych, które występują w analizowanym szeregu (prognozy wraz ze wzrostem horyzontu są coraz bliższe linii trendu). Jest to prawdopodobnie związane z faktem modelowania sezonowości jedynie przez część stacjonarną modelu (brak różnicowania sezonowego).

W prognozach wyznaczonych na bazie metody STL+ETS, podobnie jak w przypadku metody BATS, widzimy efekty dopasowania okresu krótszego niż rzeczywisty: na przykład ostatni duży spadek oraz pik w połowie roku 2005 w zbiorze testowym są ,,wyprzedzone” przez prognozę.

Wszystkie rozważane metody pozwalają również na konstrukcję przedziałów predykcyjnych (ang. prediction intervals), dzięki którym możemy m.in. uzyskać dodatkowe informacje na temat niepewności towarzyszącej wyznaczonym prognozom punktowym. Na wykresie 20 przedstawione są 95% i 80% przedziały predykcyjne skonstruowane na podstawie poszczególnych modeli.

Rysunek 20. Prognozy wraz z przedziałami predykcyjnymi skonstruowane na bazie poszczególnych modeli.

Rysunek 20. Prognozy wraz z przedziałami predykcyjnymi skonstruowane na bazie poszczególnych modeli.

Widać, że dla wszystkich metod, poza modelami TBATS i sARIMA, otrzymujemy przedziały predykcyjne podobnej postaci. W przypadku sARIMA i TBATS widoczny jest wyraźny (i zgodny z oczekiwaniem) wzrost szerokości przedziałów wraz ze wzrostem horyzontu prognozy, świadczący o większej niepewności prognoz długoterminowych. W dalszej części tego podrozdziału przeanalizujemy jeszcze nieco dokładniej własności skonstruowanych przedziałów predykcyjnych (tzn. porównamy ich szerokość oraz pokrycie empiryczne).

Przyjrzyjmy się teraz kryteriom błędów predykcji w zależności od zastosowanej metody i długości horyzontu prognozy. Na wykresie 21 przedstawiamy zależność błędu MASE (ang. Mean Absolute Scaled Error) oraz RMSE (ang. Root Mean Squared Error) od horyzontu prognozy h. W tabeli 13 przedstawione są natomiast dokładne wartości kryteriów MASE, RMSE, MAE oraz MAPE dla najdłuższego horyzontu, tzn. h=150 tygodni.

Rysunek 21. Wykres błędów MASE i RMSE w zależności od horyzontu prognozy.

Rysunek 21. Wykres błędów MASE i RMSE w zależności od horyzontu prognozy.


Tabela 13. Kryteria dokładności prognoz dla horyzontu $h$ = 150
RMSE MAE MAPE MASE
Regresja harmoniczna I 324.24 273.09 3.05 0.80
Regresja harmoniczna II 224.89 177.15 1.98 0.52
sARIMA 298.98 245.59 2.78 0.72
TBATS 212.12 168.94 1.88 0.50
BATS 238.36 189.39 2.11 0.56
STL+ETS 232.34 184.04 2.04 0.54

Z wykresów błędów MASE i RMSE (Rys.21) wynika, że na bazie metody TBATS otrzymujemy najdokładniejsze prognozy praktycznie dla wszystkich rozważanych horyzontów czasowych. Nieco gorzej wypadają modele: regresja harmoniczna II, BATS oraz STL+ETS, które charakteryzują się bardzo podobną dokładnością. Najniższą dokładność prognoz otrzymujemy na podstawie modelu regresji harmonicznej I oraz sARIMA. Dodatkowo, w przypadku sARIMA widoczny jest również bardzo szybki spadek dokładności prognoz, nawet dla umiarkowanych horyzontów (dłuższych niż pół roku).

W ocenie poprawności prognoz bardzo ważną rolę odgrywa ocena losowości reszt/błędów predykcji (czyli różnic pomiędzy wartością rzeczywistą a prognozą). Na wykresie 22 przedstawimy autokorelogramy błędów predykcji dla poszczególnych modeli.

Rysunek 22. Autokorelogramy reszt (błędów prognozy na zbiorze testowym) dla poszczególnych modeli.

Rysunek 22. Autokorelogramy reszt (błędów prognozy na zbiorze testowym) dla poszczególnych modeli.

Widzimy, że błędy prognoz różnią się między sobą pod względem obecności lub braku korelacji. Najmniej korzystnie przedstawiają się autokorelogramy dla błędów prognozy na bazie modeli regresji harmonicznej I oraz sARIMA – widzimy, że wiele informacji dotyczącej zmienności sezonowej pozostało w resztach (nie zostało wyjaśnione przez dopasowany model). Dla pozostałych modeli wyniki są zbliżone: otrzymujemy błędy prognozy bliskie białemu szumowi, z niewieloma istotnymi autokorelacjami.

Na podstawie analizy błędów prognoz możemy także określić czy różnice w jakości prognoz są statystycznie istotne. Posłużymy się do tego testem Diebolda – Mariano. W tabeli 14 przedstawiamy p-wartości dla hipotezy zerowej o braku istotnych różnic w dokładności prognoz.


Tabela 14. p-wartości dla testu Diebolda-Mariano.
regresja harmon.II TBATS BATS sARIMA STL+ETS
regresja harmon.I 0.000000 0.000000 0.000073 0.038196 0.000006
regresja harmon.II 0.163851 0.343552 0.000015 0.442894
TBATS 0.027452 0.000003 0.045025
BATS 0.003036 0.482029
sARIMA 0.000606

Jak widać,w przypadku modeli regresji harmonicznej I oraz sARIMA otrzymujemy prognozy istotnie gorsze od prognoz na bazie pozostałych metod (p-wartości są wyraźnie mniejsze od poziomu istotności 0.05). Pozostałe p-wartości dla testu D-M zdają się potwierdzać wyniki zawarte w tabeli 13. W szczególności, dokładność prognoz dla metody TBATS została uznana za statystycznie równoważną dokładności prognoz skonstruowanych na podstawie modelu regresji harmonicznej II oraz lepszą (na poziomie istotności 0.05) niż dokładność BATS i STL+ETS.

6 Porównanie przedziałów predykcyjnych

Na zakończenie porównamy również średnie długości 95% przedziałów predykcyjnych dla poszczególnych metod (patrz tabela 15) oraz empiryczne prawdopodobieństwa pokrycia (patrz tabela 16).


Tabela 15. Średnia długość przedziałów predykcyjnych dla horyzontu prognozy $h=150$.
Średnia szerokość przedziałów predykcyjnych
regresja harmoniczna I 1168.86
regresja harmoniczna II 1146.28
sARIMA 2451.06
TBATS 1997.96
BATS 1271.43
STL+ETS 1099.61

Najszersze przedziały predykcyjne otrzymujemy na bazie modelu sARIMA. Drugie miejsce (pod względem szerokości) zajmują przedziały wyznaczone dla modelu TBATS. Średnie szerokości dla pozostałych metod są zbliżone.


Tabela 16. Empiryczne prawdopodobieństwa pokrycia dla horyzontu prognozy $h=150$ i 95$\%$ przedziałów predykcyjnych.
Empiryczne prawdopodobieństwo pokrycia
regresja harmoniczna I 0.93
regresja harmoniczna II 1.00
sARIMA 1.00
TBATS 1.00
BATS 0.99
STL+ETS 0.99

Przedstawimy na wykresie 23 zależność szerokości przedziałów predykcyjnych od horyzontu prognozy $h$. W tabeli 17 przedstawiamy natomiast zależność empirycznego prawdopodobieństwa pokrycia od horyzontu prognozy $h$.

Rysunek 23. Średnia szerokość przedziałów predykcyjnych z zależności od horyzontu prognozy h.

Rysunek 23. Średnia szerokość przedziałów predykcyjnych z zależności od horyzontu prognozy h.


Tabela 17. Empiryczne prawdopodobieństwa pokrycia dla 95$\%$ przedziałów predykcyjnych, zależności od długości horyzotu $h$.
10 30 50 75 100 150
regresja harmoniczna I 1.00 1.00 0.94 0.92 0.92 0.93
regresja harmoniczna II 1.00 1.00 1.00 1.00 1.00 1.00
sARIMA 1.00 1.00 1.00 1.00 1.00 1.00
TBATS 1.00 1.00 1.00 1.00 1.00 1.00
BATS 1.00 0.97 0.98 0.99 0.99 0.99
STL+ETS 1.00 1.00 1.00 1.00 1.00 0.99

W przypadku krótkiego i średniego horyzontu prognozy ($h\leq 75$), średnie szerokości przedziałów predykcyjnych na bazie wszystkich metod, poza metodą sARIMA, są zbliżone do siebie. Dla dłuższych horyzontów dosyć gwałtownie wzrasta szerokość przedziałów predykcyjnych wyznaczonych na bazie modelu TBATS. Empiryczne prawdopodobieństwa pokrycia są wysokie i bliskie $100\%$. Jedynie w przypadku modelu regresji harmonicznej I, dla prognoz o dłuższym horyzoncie pokrycie empiryczne jest niższe od nominalnego poziomu ufności równego 0.95.

7 Porównanie złożoności obliczeniowej metod

Ostatnim z aspektów analizowanych metod na który zwrócimy uwagę, jest czas potrzebny na przeprowadzenie obliczeń. Jest to o tyle istotne, że w przypadku długich szeregów czasowych niektóre modele mogą potrzebować długiego czasu na konstrukcję. Czas wykorzystany przez każdą z metod zaprezentowany jest w tabeli 18.


Tabela 18. Czas obliczeń (w sekundach) dla poszczególnych metod
czas[s]
Regresja harmoniczna 0.28
TBATS 2.75
BATS 0.90
sARIMA 21.64
STL+ETS 0.00

Widoczna jest wyraźna różnica między modelem sARIMA a pozostałymi, wynika ona jednak z zastosowania automatycznego wyboru modelu na podstawie (złożonego obliczeniowo) przeglądu zupełnego. Drugą metodą, której czas wykonania odbiega od pozostałych, jest metoda TBATS, jednak różnica między nią a pozostałymi metodami nie jest już tak duża.

8 Podsumowanie

  • Powyższa analiza pozwoliła zaobserwować pewną przewagę modeli TBATS oraz BATS (dedykowanych złożonej sezonowości) nad pozostałymi modelami. Dla prawie wszystkich horyzontów czasowych dokładności prognoz skontruowanych na podstawie tych modeli były jednymi z najwyższych (dla metody TBATS najwyższe). Dla niektórych horyzontów klasyczny model regresji harmonicznej z ustaloną częstością pozwalał jednak na wyznaczenie nieco dokładniejszych prognoz niż model BATS (mniejsze wartości błędów MASE i RMSE).
  • Porównanie między modelami TBATS oraz BATS wskazuje model TBATS jako lepszy pod kątem dokładności prognoz.

  • Pewnym zaskoczeniem może być wysoka dokładność, jaką charakteryzowały się również prognozy wyznaczone na bazie modelu STL+ETS, który podobnie jak model BATS nie uwzględnia poprawnie ułamkowej długości okresu odpowiadającego sezonowości rocznej. Dokładność prognoz na bazie modelu STL+ETS była nieco gorsza od prognoz dla modelu TBATS, ale porównywalna do modeli BATS i regresji harmonicznej II.

  • Dla krótkich horyzontów, prognozy na bazie wszystkich modeli, poza I-szym wariantem regresji harmonicznej oraz sARIMA, charakteryzują się podobną dokładnością.

  • Spośród analizowanych modeli najmniej dokładne prognozy otrzymujemy na bazie (wybranego automatycznie) sezonowego wariantu modelu ARIMA (model sARIMA) oraz wariantu regresji harmonicznej z automatycznym wyborem statystycznie istotnych częstości fourierowskich (regresja harmoniczna I). W pierwszym przypadku spowodowane było to prawdopodobnie nieuwzględnieniem w modelu różnicowania sezonowego, przez co zależności sezonowe były "opisane" jedynie przez stacjonarną część modelu (wahania sezonowe nie były prawidłowo odzwierciedlone w prognozie). W drugim przypadku, uwzględnienie w modelu dodatkowych częstości (niezwiązanych z występującą w danych sezonowością roczną) wpłynęło na znaczne obniżenie dokładności prognoz. Oba te przypadki wskazują na typowe problemy, z którymi można się spotkać stosując automatyczne metody wyboru optymalnego modelu.

  • Porównanie szerokości i empirycznych prawdopodobieństw pokrycia dla skonstruowanych przedziałów predykcyjnych nie wskazuje żadnej konkretnej metody jako najlepszej. Negatywnie wyróżniają się jedynie: model sARIMA (największa szerokość przedziałów predykcyjnych) oraz model regresji harmonicznej I (niskie pokrycie empiryczne dla dłuższych horyzontów prognozy).

Spróbuj ponownie