24.03

2015

Goście hotelowi na Dolnym Śląsku — analizy i prognozy

Autor: adam

Tym razem wpis dotyczy praktycznej analizy szeregów czasowych. Powstał podczas stażu Hanny Loch i Piotra Kopszaka i jest efektem naszej wspólnej pracy.

Na warsztat wzięliśmy dane (z ostatnich kilku lat) dotyczące liczby turystów korzystających z
noclegów na terenie Dolnego Śląska. Przyjrzymy się podstawowym własnościom analizowanego szeregu, takim jak tendencja długoterminowa (trend) i wahania sezonowe (sezonowość). Porównamy dokładność różnych metod stosowanych do konstrukcji prognoz. Zastanowimy się także czy i na jakie pytania biznesowe można by odpowiedzieć przeprowadzając taką analizę.

Artykuł skoncentrowany jest na aspektach praktycznych i nie wymaga od czytelnika znajomości wykorzystywanej metodologii. Powinien przydać się m.in. tym osobom, które wcześniej nie zajmowały się praktyczną analizą i prognozowaniem szeregów czasowych, a które chciałyby poznać potencjalne korzyści związane ze stosowaniem takich metod.

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 Wprowadzenie

Poniższe opracowanie zawiera przegląd wybranych metod prognozowania szeregów czasowych na przykładzie konkretnych danych — tak zwane case study. Naszym celem jest zaprezentowanie nie tyle wszystkich możliwych metod, co tych najbardziej odpowiednich dla danego szeregu i pokazanie korzyści wynikających z odpowiedniego doboru narzędzi analitycznych.

Dane, z których będziemy korzystać, dotyczą liczby gości w obiektach hotelowych na terenie Dolnego Śląska. Są to dane miesięczne, zgromadzone od stycznia 2009 roku do marca 2014 roku; można je pobrać ze strony Banku Danych Lokalnych (stat.gov.pl/bdl). Wszystkie analizy zostały przeprowadzone w środowisku R.

1.1 Pytania biznesowe i analiza danych

Co może wynikać z analizy takich danych?

  1. Prognozowanie liczby gości w kolejnych okresach.
    Dla właścicieli obiektów hotelowych liczba gości jest z oczywistych względów bardzo ważną informacją. Prawdopodobnie wielu z nich intuicyjnie stosuje proste metody analityczne – na przykład szacują liczbę klientów w maju na podstawie danych z maja zeszłego roku. Do tego podejścia wrócimy w dalszej części artykułu. Pokażemy też, jak prognozować lepiej (bardziej dokładnie) niż standardowymi metodami.

  2. Analiza tendencji sezonowych i długoterminowych.
    Często dzięki analizom szeregów czasowych możemy zobaczyć, jaka jest kondycja w danym sektorze — czy możemy doświadczyć pewnej stagnacji, czy raczej jesteśmy w okresie dobrej koniunktury. Możliwe jest również zidentyfikowanie pewnych powtarzających się wzorców sezonowych – możemy na przykład spodziewać się większej liczby turystów w miesiącach letnich niż w zimie. Nie zawsze jednak jest to tak przewidywalne, czasami dopiero szczegółowa analiza szeregu ukazuje nam zaskakujący, sezonowy wzór, niewidoczny na pierwszy rzut oka.

  3. Podejmowanie decyzji dotyczących pokrewnych przedsięwzięć.
    Liczba osób nocujących w obiektach hotelowych oddaje po części liczbę turystów odwiedzających dany region — owszem, nie uwzględnia wszystkich, ale potrafi powiedzieć coś o wzroście bądź spadku liczby gości. Te informacje mogą być wykorzystane do przygotowania oferty turystycznej i dobrego rozplanowania czasowego wydarzeń kulturalnych.

1.2 Podstawowe własności danych

Pierwszym krokiem w analizie jest dokładne przyjrzenie się wykresowi szeregu czasowego. Wykres naszych danych dotyczących liczby gości hotelowych przedstawia rysunek 1.

Rysunek 1. Liczba gości w obiektach hotelowych w okresie 01.2009-03.2014

Rysunek 1. Liczba gości w obiektach hotelowych w okresie 01.2009-03.2014

Na wykresie szeregu bardzo wyraźnie widać wspomnianą wcześniej sezonowość — co roku w miesiącach zimowych jest zdecydowanie mniej turystów niż w miesiącach letnich; poza tym dla każdej z pór roku zachodzą pewne różnice pomiędzy poszczególnymi miesiącami, różnice, które są analogiczne we wszystkich latach.

Na rysunku 1 możemy też zauważyć, że dostrzegalny na wykresie wzorzec sezonowy nie jest identyczny dla kolejnych lat. Na rysunku 2 został przedstawiony wykres sezonowy, pokazujący jak w danym roku kształtowała się liczba gości na przestrzeni miesięcy. Widać na nim bardzo wyraźnie powtarzający się sezonowy układ; nie jest on jednak do końca taki, jak intuicyjnie byśmy przewidywali — rzeczywiście, następuje wzrost turystów w miesiącach letnich, jednak maksima w poszczególnych latach przypadają na maj bądź sierpień, ze spadkiem w czerwcu i lipcu.

Rysunek 2. Liczba gości na przestrzeni lat względem poszczególnych miesięcy

Rysunek 2. Liczba gości na przestrzeni lat względem poszczególnych miesięcy

Patrząc na wykres sezonowy przedstawiony na rysunku 2 możemy dostrzec jedno odstępstwo od corocznego, sezonowego wzorca — w czerwcu 2012 roku nastąpił niespodziewanie duży spadek liczby gości (w pozostałych latach w czerwcu następowało jedynie drobne zmniejszenie). Warto przypomnieć, że to właśnie wtedy we Wrocławiu były rozgrywane mecze w ramach rozgrywek Euro 2012, więc oczekiwalibyśmy wręcz wzrostu liczby osób nocujących w hotelach na terenie Dolnego Śląska. Z drugiej strony, przez zamieszanie związane z mistrzostwami być może wiele osób, które zazwyczaj przyjeżdżają do Wrocławia lub okolic w celach biznesowych, przełożyło swój wyjazd. Ponadto nie wszystkie mecze przyciągnęły tę samą liczbę nocujących turystów — na przykład Czesi, których drużyna narodowa grała we Wrocławiu, nie potrzebowali noclegu z powodu względnie niewielkiej odległości.

Kolejny wykres (przedstawiony na rysunku 3) pokazuje, jak zmieniała się liczba gości w danym miesiącu na przestrzeni lat. Dla większości miesięcy możemy dostrzec wzrost w pierwszych latach z drobnym spadkiem w ostatnim, 2013 roku. Jest to podstawa do wspomnianej wcześniej analizy tendencji długoterminowych — być może po 2012 roku nastąpiła pewna stagnacja liczby gości na terenie Dolnego Śląska (tak duża liczba gości w 2012 może też być długofalowym efektem Euro2012).

Rysunek 3. Liczba gości w poszczególnych miesiącach na przestrzeni lat

Rysunek 3. Liczba gości w poszczególnych miesiącach na przestrzeni lat

Wszystkie trzy omówione wykresy sugerowały, że w analizowanych danych obecny jest wyraźny trend wzrostowy. Możemy zastanawiać się teraz, jak w prosty sposób modelować tę tendencję. Najprostszym wyborem jest przyjęcie trendu liniowego, czyli założenie, że w każdym roku liczba gości wzrasta o pewną ustaloną liczbę — na przykład co roku z obiektów hotelowych korzysta o 10 000 osób więcej. Jednak to zmniejszenie liczby gości w 2012 roku mogłoby sugerować bardziej skomplikowaną postać trendu, na przykład trend kwadratowy. W dalszej części zobaczymy, czy nasze wstępne rozważania są słuszne.

1.3 Inne własności danych

Nasze dane są danymi miesięcznymi — jest to sumaryczna liczba gości w danym miesiącu. W tym wypadku może pojawić się efekt kalendarzowy związany ze zmienną długością poszczególnych miesięcy — liczba dni między miesiącami różni się nawet o 10% (na przykład między styczniem a lutym), co może skutkować błędną interpretacją wyników analizy (możemy na przykład odnaleźć fałszywy wzorzec sezonowy). W bardzo prosty sposób możemy dokonać korekty tak, by uwzględnić zmienną liczbę dni w poszczególnych miesiącach (tzw. month length adjustment).

Wykres 4 przedstawia porównanie wykresów szeregu czasowego przed i po dostosowaniu danych we wspomniany sposób. Widać na nim zdecydowaną zmianę kształtu szeregu w lutym każdego roku.

Rysunek 4. Wykres szeregu czasowego przed i po uwzględnieniu zmiennej długości miesięcy (*month length adjustment*)

Rysunek 4. Wykres szeregu czasowego przed i po uwzględnieniu zmiennej długości miesięcy (month length adjustment)

Rysunek 5. Wykres wzorca sezonowego przed i po uwzględnieniu zmiennej długości miesięcy (*month length adjustment*)

Rysunek 5. Wykres wzorca sezonowego przed i po uwzględnieniu zmiennej długości miesięcy (month length adjustment)

Warto również porównać wykresy sezonowe przed i po zastosowaniu wspomnianej transformacji (korekty kalendarzowej). Porównanie to zostało przedstawione na wykresie 5. Zauważmy, że po dostosowaniu danych pojawił się znaczny wzrost liczby gości w lutym — może być to efekt ferii zimowych, którego nie byliśmy w stanie wcześniej zauważyć.

1.4 Identyfikacja regularnych wzorców w danych

Jak wspomnieliśmy, w naszych danych możemy dostrzec pewien trend, świadczący o długotrwałej tendencji wzrostowej dotyczącej liczby gości, oraz wzorzec sezonowy, który wskazuje na roczne wahania. Możemy przyjąć, że szereg czasowy składa się z tych dwóch elementów oraz z losowych fluktuacji. Takie podejście, nazywane dekompozycją szeregu, jest bardzo intuicyjne.

Model dekompozycji możemy zapisać symbolicznie:
$$y_t = T_t + S_t + E_t \ \mbox{(dekompozycja addytywna)}$$ bądź $$y_t = T_t \times S_t \times E_t \ \mbox{(dekompozycja multiplikatywna)},$$ gdzie $T_t$ – trend, $S_t$ – sezonowość, $E_t$ – losowe fluktuacje.

Jak to rozumieć? Popatrzmy na nasze dane — stwierdziliśmy obecność pewnego trendu wzrostowego, który odpowiada długotrwałemu wzrostowi liczby gości na Dolnym Śląsku z roku na rok. Ponadto występuje wyraźny wzorzec sezonowy (sezonowość), który powtarza się co roku i odpowiada za to, że m.in. w maju jest więcej turystów niż w grudniu. Jednak nie może to być wszystko, ponieważ są też inne czynniki wpływające na to, czy turyści chcą odwiedzić dany region (chociażby pogoda) — te mieszczą się w części losowej, związanej z nieregularnymi fluktuacjami.

Podejście to jest często używane ze względu na prostą konstrukcję prognoz — po zidentyfikowaniu poszczególnych składowych możemy wyznaczyć ich prognozę i po zsumowaniu (bądź wymnożeniu) otrzymujemy prognozę dla wyjściowego szeregu.

Wybór wariantu dekompozycji — addytywna czy multiplikatywna — zależy od danych. Model addytywny jest bardziej odpowiedni, jeżeli wielkość sezonowych fluktuacji nie jest zależna od poziomu szeregu czasowego — czyli bez względu na to, ilu turystów rocznie odwiedzi Dolny Śląsk, różnica pomiędzy poszczególnymi miesiącami będzie (mniej więcej) taka sama. W przeciwnym razie model multiplikatywny może okazać się bardziej prawidłowy. Warto wspomnieć, iż alternatywą dla modelu multiplikatywnego jest zastosowanie modelu addytywnego do danych zlogarytmowanych.

Pierwszym krokiem w dekompozycji jest zidentyfikowanie trendu. Wspominaliśmy już wyżej, że możemy próbować dopasować na przykład trend liniowy czy kwadratowy, jednak była to tylko próba odgadnięcia jego postaci bezpośrednio z wykresu. Metoda, która może przyjść nam z pomocą w dokładniejszej identyfikacji tego elementu to średnia krocząca (moving average). Posłuży nam ona do wyodrębnienia z posiadanych danych składowej trendu. Średnią kroczącą w chwili $t$ określamy jako średnią arytmetyczną z pewnej liczby obserwacji poprzedzających i następujących (o ile oczywiście mamy odpowiednie poprzedzające i następujące obserwacje, co uściślimy poniżej). Średnią kroczącą rzędu $m$ określamy jako $$ \hat{T_t} = \frac{1}{m} \sum_{i=-k}^{k} y_{t+i},$$ gdzie $m=2k+1$ (tzw. rząd ruchomej średniej).

Oczywiście postać dopasowanego trendu zależy od wybranej wartości parametru $m$, który pełni rolę parametru wygładzającego. Warto porównać otrzymywane postaci trendu dla różnych rzędów $m$. Jeżeli zastosujemy średnią kroczącą do naszych danych, otrzymamy rezultaty pokazane na wykresie 6 (wybieramy $m=5$, $m=7$$m=12$).

Rysunek 6. Dopasowanie trendu przy użyciu średniej kroczącej

Rysunek 6. Dopasowanie trendu przy użyciu średniej kroczącej

Jak widać, średnie kroczące rzędu 5 i 7 oddają bardzo dobrze sezonowe wzrosty i spadki, ale nie to było naszym celem — chcieliśmy wyznaczyć ogólną tendencję rozwojową, co lepiej oddaje średnia ruchoma rzędu 12. W praktyce bardzo często stosuje się zasadę, że jeśli dane są sezonowe, to rząd ruchomej średniej jest równy pełnemu okresowi (u nas właśnie 12 miesięcy).

Warto zauważyć, że na wykresie 6 wszystkie dopasowane średnie kroczące są krótsze niż oryginalne dane. Dzieje się tak dlatego, iż były one dopasowane przy użyciu symetrycznej średniej kroczącej — uwzględniana była równa liczba obserwacji poprzedzających i następujących; dla skrajnych momentów czasowych (na początku i na końcu analizowanego okresu) nie mamy oczywiście poprzedzających lub następujących obserwacji, co nie pozwala na wyznaczenie wartości średniej kroczącej.

Metoda średniej kroczącej może być z powodzeniem stosowana do wstępnej identyfikacji prostego modelu trendu. Bardzo często rezultat, który otrzymamy przy zastosowaniu tej metody, może nam zasugerować jakiś prosty parametryczny model dla trendu, na przykład trend linowy czy kwadratowy.

Jeżeli mamy już wyznaczony trend, to po jego wyeliminowaniu (usunięciu) z danych możemy wyznaczyć składową sezonową (sezonowość), którą tworzą tzw. indeksy sezonowe — wartości przypisane każdemu z miesięcy, informujące o tym, o ile obserwacja z danego miesiąca różni się od średniej. W najprostszym przypadku, otrzymuje się je w następujący sposób: z danych z usuniętym już trendem wyznacza się średnią wartość dla konkretnego miesiąca (liczoną na bazie wszystkich lat) — czyli na przykład indeks sezonowy dla maja to średnia liczba gości (po usunięciu trendu) we wszystkich latach w maju.

Zobaczmy również, czy i jak bardzo zmieniłaby się składowa sezonowa, gdybyśmy wyznaczali indeksy sezonowe z danych po uwzględnieniu różnej liczby dni w poszczególnych miesiącach. Porównanie obu wykresów indeksów sezonowych przedstawione jest na wykresie 7. Zauważmy, że po dostosowaniu danych zmienił się indeks sezonowy na przykład w lutym — uwzględnienie mniejszej liczby dni w tym miesiącu sprawiło, iż średnia wartość w tym miesiącu nie odstaje tak bardzo od średniej całego szeregu.

Rysunek 7. Wykres indeksów sezonowych w omawianych danych

Rysunek 7. Wykres indeksów sezonowych w omawianych danych

Skoro zidentyfikowaliśmy już składową sezonową, to możemy wyodrębnić z naszych danych losowe fluktuacje, czyli reszty po usunięciu z nich trendu i sezonowości. Ta składowa przedstawiona jest na wykresie 8. Widać na nim jeden znaczny spadek, w czerwcu 2012 — są to wspomniane wcześniej rozgrywki Euro2012.

Rysunek 8. Losowe fluktuacje liczby gości -- dane po usunięciu trendu i sezonowości

Rysunek 8. Losowe fluktuacje liczby gości — dane po usunięciu trendu i sezonowości

Na zakończenie, wszystkie opisane powyżej składowe analizowanego szeregu (tzn.: trend, sezonowość i losowe fluktuacje) możemy przedstawić na jednym wykresie — rysunek 9.

Rysunek 9. Regularne wzorce w danych - dekompozycja

Rysunek 9. Regularne wzorce w danych – dekompozycja

2 Konstrukcja prognoz

Jednym z podstawowych zadań analizy szeregów czasowych jest konstrukcja prognoz dla kolejnych okresów, na podstawie zidentyfikowanych własności danego szeregu. Poniżej zaprezentujemy wyniki prognozowania, uzyskane na bazie różnych metod i modeli.

By móc ocenić dokładność otrzymanych prognoz, podzielimy nasze dane na podzbiór uczący i testowy, tj. część danych posłuży nam do zbudowania modelu, na podstawie których dokonamy prognoz na następny okres. Następnie przewidywane wartości zostaną porównane z prawdziwymi danymi z tego okresu — podzbiorem testowym.

Przyjmijmy, że podzbiór uczący to dane z okresu styczeń 2009 — marzec 2013, a podzbiór testowy to kwiecień 2013 – marzec 2014 (czyli prognozujemy wartości na jeden pełny rok).

2.1 Ocena dokładności prognoz

By dobrze porównać prognozy uzyskane różnymi metodami, należy przyjąć pewne kryteria oceniające popełniony błąd, czyli różnicę między prognozą a rzeczywistymi wartościami szeregu. W praktyce stosuje się różne kryteria (błędy względne, procentowe, skalowane, itp.); w tym opracowaniu będziemy korzystać z następujących dwóch:

  1. MAPE – mean absolute percentage error, czyli średni bezwzględny błąd procentowy. Wartość MAPE mówi nam, o ile procent średnio różniła się wartość prognozowana od rzeczywistej.
  2. MASE – mean absolute scaled error, czyli średni bezwzględny błąd skalowany. Jest to porównywanie otrzymanego błędu z błędem, który popełnilibyśmy konstruując prognozy za pomocą metody naiwnej (opisanej poniżej).

Będziemy wyliczać wartości obu rodzaju błędów dla każdej z przygotowywanych prognoz tak, by móc później wybrać najlepszą (obarczoną najmniejszym błędem) metodę konstrukcji prognoz.

2.2 Proste metody prognozowania

Zacznijmy od najprostszej metody prognozowania — metody naiwnej. W tym przypadku jako prognozę dla kolejnych (nowych) okresów przyjmujemy ostatnią z zaobserwowanych wartości. W naszym przypadku prognozowalibyśmy więc, że od kwietnia 2013 do marca 2014, w każdym z miesięcy Dolny Śląsk odwiedzi tyle samo gości hotelowych, ile zaobserwowano w marcu 2013.

Jednak prognozy wyznaczone w ten sposób nie biorą pod uwagę występującego trendu (tendencji długoterminowej), o którym wspominaliśmy wyżej. Pozwala na to metoda uwzględniająca dryf (random walk with drift), będąca modyfikacją metody naiwnej. W przypadku tej metody, zamiast przyjmować jako prognozę ostatnią zaobserwowaną wartość szeregu, przyjmujemy tę wartość zwiększoną o pewną stałą, tak zwaną stałą dryfu. Jest ona równa średniej wartości przyrostów szeregu (różnic pomiędzy sąsiadującymi obserwacjami). Graficznie predykcję tę możemy traktować jak przedłużenie (ekstrapolowanie) prostej poprowadzonej przez pierwszą i ostatnią wartość szeregu. Na wykresie 10 przedstawiono wynik prognozowania z wykorzystaniem metody uwzględniające dryf.

Rysunek 10. Prognoza otrzymana za pomocą dryfu

Rysunek 10. Prognoza otrzymana za pomocą dryfu

Wykres 10 sugeruje, że prognoza ta nie jest precyzyjna — są duże różnice pomiędzy rzeczywistymi danymi a wartościami przewidywanymi. Aby dokładnie ocenić jakość predykcji, wyznaczamy wspomniane już kryteria MAPE i MASE (tablica 1). Wartość MASE jest znacznie większa od 1, co oznacza, że wykorzystując sezonową metodę naiwną możemy spodziewać się lepszych rezultatów (przypomnijmy, że w przypadku danych sezonowych, to właśnie sezonowa metoda naiwna jest metodą referencyjną, wykorzystywaną do wyznaczenia błędu MASE).

Tabela 1. Błędy predykcji dla metody uwzględniającej dryf
MAPE MASE
Metoda uwzględniająca dryf 17.19 3.28

Zauważmy, że predykcja na bazie metod naiwnych jest tak niedokładna ze względu na sezonowość danych — gdyby we wszystkich miesiącach następował równomierny wzrost szeregu, przewidziane przez nas wartości byłyby zdecydowanie bliższe rzeczywistym. Jest to zatem przykład na to, że konstruując prognozy powinniśmy zwracać baczną uwagę na wszystkie regularne składowe występujące w szeregu.

Prostą metodą predykcji, która może być stosowana dla danych sezonowych, jest tzw. sezonowa metoda naiwna. W metodzie naiwnej jako prognozę w danej chwili przyjmuje się wartość szeregu z chwili poprzedniej; jeżeli chcemy uwzględnić sezonowość, jako wartości prognozowane w danym roku przyjmujemy wartości z roku poprzedniego, oczywiście z odpowiedniego miesiąca. W ten sposób prognozowana liczba gości w marcu 2013 będzie dokładnie taka, jaka była liczba gości w marcu 2012. Wykres 11 przedstawia prognozę wyznaczoną na bazie sezonowej metody naiwnej.

Rysunek 11. Prognoza otrzymana za pomocą sezonowej metody naiwnej

Rysunek 11. Prognoza otrzymana za pomocą sezonowej metody naiwnej

Wyznaczmy również wartości MAPE i MASE i porównajmy je z wartościami błędów dla metody uwzględniającej dryf. Porównanie zostało przedstawione w tabeli 2. Zwróćmy uwagę, iż MASE dla sezonowej metody naiwnej jest mniejsze od 1, choć to właśnie sezonowa metoda naiwna jest tą, do której się odnosimy wyliczając tę miarę błędu. Jednak dokładniej rzecz ujmując, MASE jest porównaniem błędu danej metody z błędem, który otrzymujemy dla sezonowej metody naiwnej na zbiorze uczącym, natomiast w tabeli jest przedstawiony błąd predykcji popełniony na zbiorze testowym.

Tabela 2. Błędy predykcji dla metody uwzględniającej dryf i metody sezonowej naiwnej
MAPE MASE
Metoda uwzględniająca dryf 17.19 3.28
Sezonowa metoda naiwna 4.28 0.70

Widać, że błąd procentowy (MAPE) jest czterokrotnie mniejszy. Nasuwa się jednak pytanie, czy predykcja z wykorzystaniem sezonowej metody naiwnej jest najlepszym rezultatem jaki możemy otrzymać. W kolejnych rozdziałach przedstawimy bardziej zaawansowane metody i porównamy je z wynikiem otrzymanym dla sezonowej metody naiwnej, by móc ocenić, czy (i w jakim stopniu) zastosowanie zaawansowanych metod poprawi dokładność skonstruowanych prognoz.

2.3 Konstrukcja prognoz na bazie dekompozycji

Dekompozycja szeregu (na regularne wzorce i losowe fluktuacje) pozwala nam również na łatwe wyznaczenie prognoz. Zobaczymy zatem, jak dokładne będą prognozy skonstruowane na podstawie różnych modeli dekompozycji.

Na początek załóżmy, że szereg składa się z trendu wielomianowego oraz czynnika losowego (pomińmy sezonowość). Rozpatrzmy trend liniowy oraz kwadratowy — wyznaczone prognozy są przedstawione na wykresie 12.

Rysunek 12. Prognoza na podstawie trendu liniowego i kwadratowego

Rysunek 12. Prognoza na podstawie trendu liniowego i kwadratowego

Jak widać, przyjęcie modelu dekompozycji, w którym występuje jedynie trend, bez sezonowości, nie prowadzi do dobrych wyników (błędy predykcji przedstawione są w tabeli 3). Prognozy na bazie modelu dekompozycji uwzględniającego sezonowość przedstawia rysunek 13.

Rysunek 13. Prognoza na podstawie modelu z trendem oraz sezonowością

Rysunek 13. Prognoza na podstawie modelu z trendem oraz sezonowością

Wartości błędów MAPE i MASE dla prognoz wyznaczonych na bazie dekompozycji uwzględniającej sezonowość są przedstawione w tablicy 3.

Tabela 3. Błędy predykcji dla modelu dekompozycji
MAPE MASE
Trend liniowy 18.19 2.95
Trend kwadratowy 14.61 2.65
Sezonowość + trend liniowy 4.90 0.75
Sezonowość + trend kwadratowy 4.78 0.73

Oba modele dekompozycji (tzn. trend liniowy z sezonowością oraz trend kwadratowy z sezonowością) prowadzą do niemal identycznych prognoz. Ponadto widzimy, jak istotne podczas konstrukcji prognoz jest uwzględnianie występujących w szeregu wahań sezonowych — wartości MASE dla prognoz z samym trendem są większe niż 1, czyli prognozy te są mniej dokładne niż uzyskane na bazie metody referencyjnej (czyli sezonowej metody naiwnej).

2.4 Wygładzanie wykładnicze

Kolejnym algorytmem, który możemy wykorzystać do przewidywania przyszłych wartości szeregu czasowego jest wygładzanie wykładnicze. Algorytm ten polega na przedstawieniu każdej obserwacji jako ważonej średniej poprzednich obserwacji. Skąd więc przymiotnik "wykładnicze"? Otóż wagi przypisywane poszczególnym przeszłym wartościom szeregu maleją wykładniczo — im starsza (bardziej oddalona w czasie) obserwacja, tym mniejsza odpowiadająca jej waga.

Najprostszą wersją wygładzania wykładniczego jest algorytm Simple Exponential Smoothing (SES) — stosowany do danych, w których nie ma trendu ani sezonowości. Zatem nie powinniśmy stosować tej wersji wygładzania do naszych danych związanych z liczbą turystów w hotelach dolnośląskich, jednakże jest to dobry przykład umożliwiający zaprezentowanie idei działania metody (bardziej precyzyjny opis algorytmu SES znajduje się w dodatku).

Jeżeli w analizowanych danych występuje trend, należy zastosować metodę Holta, będącą modyfikacją SES. Jeśli występuje również sezonowość, odpowiedni jest algorytm Holta-Wintersa. Niemniej jednak, możemy spróbować zastosować do naszych danych również najprostszy algorytm (SES), aby zilustrować konieczność zastosowania metody uwzględniającej zarówno trend, jak i sezonowość, czyli metody Holta-Wintersa.

Na wykresie 14 przedstawiamy wyniki prognozy przy użyciu metody SES, natomiast na wykresie 15 pokazane są prognozy wyznaczone na bazie algorytmu Holta-Wintersa.

Rysunek 14. Prognoza na podstawie prostego wygładzania wykładniczego

Rysunek 14. Prognoza na podstawie prostego wygładzania wykładniczego

Rysunek 15. Prognoza na podstawie metody Holta-Wintersa

Rysunek 15. Prognoza na podstawie metody Holta-Wintersa

Wartości MAPE i MASE są przedstawione w tabeli 4. Jak można zauważyć, wartość MASE dla metody Holta-Wintersa jest mniejsza niż 1 (dokładnie 3.4106908), czyli błędy popełnione przy użyciu tej metody są zdecydowanie mniejsze niż te popełnione przy użyciu sezonowej metody naiwnej (na podzbiorze uczącym). Z kolei MASE dla SES wynoszące 0.3929408 wskazuje, że ta metoda była gorsza niż referencyjna. Wartości błędów MAPE i MASE dla algorytmu Holta-Wintersa jednoznacznie pokazują więc, iż prognozy skonstruowane na bazie tego algorytmu są najlepsze z dotychczasowo otrzymanych.

Tabela 4. Błędy predykcji przy użyciu wygładzania wykładniczego
MAPE MASE
Algorytm Holta-Wintersa 2.50 0.39
Proste wygładzanie wykładnicze 17.85 3.41

Warto również wspomnieć, iż w przypadku algorytmu Holta-Wintersa możliwy jest wybór jednego z dwóch wariantów sezonowości, tzn. wariantu addytywnego lub multiplikatywnego. Przedstawione powyżej prognozy otrzymaliśmy na bazie addytywnej wersji algorytmu, aczkolwiek różnice w dokładności dla obu wariantów były bardzo małe.

2.5 Prognozowanie na bazie modelu ARIMA

Innym, klasycznym i często wykorzystywanym podejściem do prognozowania szeregów czasowych jest model ARIMA. Jego pełna nazwa to AutoRegressive Integrated Moving Average. W modelu tym zakładamy, że wartość szeregu (być może po zastosowaniu operacji różnicowania) w danej chwili $t$ jest liniową kombinacją wartości szeregu (część autoregresyjna) oraz zakłóceń (część związana ze średnią ruchomą), w pewnych chwilach poprzedzających $t$.

Wybór odpowiedniej postaci modelu nie jest prosty. W pakiecie R istnieje jednak funkcja pozwalająca na automatyczny wybór optymalnego modelu na podstawie odpowiednich kryteriów oceniających jakość dopasowania modelu do danych. Na rysunku 16 przedstawiony jest wynik predykcji za pomocą sezonowego modelu ARIMA. Poniżej, w tabeli 5 prezentujemy błędy dla otrzymanego modelu. Widzimy, że także ta prognoza jest lepsza od otrzymanej sezonową metodą naiwną (MASE mniejsze od 1).

Należy również dodać, że uzyskane prognozy otrzymaliśmy dopasowując model do zlogarytmowanych danych. W przypadku dopasowania modelu ARIMA dla oryginalnych danych dokładność skonstruowanych prognoz była mniejsza.

Rysunek 16. Prognoza na bazie modelu sARIMA

Rysunek 16. Prognoza na bazie modelu sARIMA

Tabela 5. Błędy predykcji przy użyciu modelu ARIMA
MAPE MASE
Model ARIMA 2.89 0.46

2.6 Która metoda jest najlepsza?

Spójrzmy, jak wyglądają wartości błędów MAPE i MASE w przypadku każdej z użytych wyżej metod (tablica 6).

Tabela 6. Porównanie błędów predykcji
MAPE MASE
Metoda uwzględniająca dryf 17.19 3.28
Sezonowa metoda naiwna 4.28 0.70
Trend liniowy 18.19 2.95
Trend kwadratowy 14.61 2.65
Sezonowość + trend liniowy 4.90 0.75
Sezonowość + trend kwadratowy 4.78 0.73
Proste wygładzanie wykładnicze 17.85 3.41
Algorytm Holta-Wintersa 2.50 0.39
Model ARIMA 2.89 0.46

W wypadku analizowanych przez nas danych, oba kryteria oceny dokładności prognoz wyraźnie wskazują na przewagę dwóch metod: Holta-Wintersa i sARIMA. Są to zaawansowane metody, więc fakt, że okazały się lepsze od pozostałych, nie dziwi. Wyłania się też grupa metod, o których z całkowitą pewnością można powiedzieć, że nie nadają się do naszych danych — są to modele uwzględniające wyłącznie trend (zarówno liniowy jak i kwadratowy), bez sezonowości, oraz proste wygładzanie wykładnicze (SES). Metodami konstrukcji prognoz, które uplasowały się pomiędzy, były metody prostej dekompozycji uwzględniające zarówno trend, jak i sezonowość, oraz sezonowa metoda naiwna. Nie były rzecz jasna tak dokładne, jak te najbardziej zaawansowane, jednak biorąc pod uwagę ich prostotę, można ich dokładność uznać za całkowicie zadowalającą. Ich niewątpliwym atutem jest łatwa i intuicyjna interpretacja wyników, co może być dużą zaletą w rozwiązaniach biznesowych, nawet za cenę odrobinę mniejszej dokładności.

W tabeli 6 wszystkie błędy wyliczone były dla 12-miesięcznego horyzontu prognozy. Natomiast w tabeli 7 pokazujemy również błędy dla 6-miesięcznego horyzontu. Widzimy, że ponownie metoda Holta-Wintersa i sARIMA dały najlepsze wyniki. Możemy zauważyć również, że część modeli lepiej sprawdza się dla krótszego horyzontu czasowego (h=6), a część dla dłuższego (h=12). Dla przykładu, porównując dokładność prognoz dla modeli wielomianowych zauważamy, że dla krótszego horyzontu ($h=6$) lepszy okazał się model liniowy, a już dla $h=12$ model kwadratowy. Ważny jest zatem wybór odpowiedniego modelu dla konkretnego zastosowania — warto korzystać z innego modelu do przygotowywania prognoz na najbliższy miesiąc czy pół roku, a innego do predykcji na cały następny rok.

Tabela 7. Porównanie błędów predykcji dla różnych horyzontów czasowych
h = 6 | MAPE h = 6 | MASE h = 12 | MAPE h = 12 | MASE
Metoda uwzględniająca dryf 27.67 3.69 17.19 3.28
Sezonowa metoda naiwna 3.20 0.41 4.28 0.70
Trend liniowy 14.99 1.96 18.19 2.95
Trend kwadratowy 19.67 2.63 14.61 2.65
Sezonowość + trend liniowy 2.90 0.36 4.90 0.75
Sezonowość + trend kwadratowy 2.84 0.35 4.78 0.73
Proste wygładzanie wykładnicze (SES) 28.48 3.80 17.85 3.41
Algorytm Holta-Wintersa 1.76 0.24 2.50 0.39
Model ARIMA 1.95 0.26 2.89 0.46

A Dodatki

A.1 Kryteria błędów zastosowane w analizie

  1. MAPE – mean absolute percentage error, czyli średni bezwzględny błąd procentowy. Można go wyliczyć według następującego wzoru: $$MAPE = \frac{1}{h} \sum_{i=1}^{h} \frac{e_i}{y_i} \cdot 100\%,$$ gdzie $y_i$ to $i$-ta obserwacja, a $e_i$ to błąd predykcji $i$-tej obserwacji (bezwzględna wartość różnicy wartości przewidywanej i wartości obserwowanej).
  2. MASE – mean absolute scaled error, czyli średni bezwzględny skalowany błąd. Można go wyliczyć według następującego wzoru: $$MASE = \frac{1}{h} \sum_{i=1}^{h} \frac{e_i}{\frac{1}{n-s}\sum_{j=2}^{n}|y_j-y_{j-s}|},$$ gdzie $y_i$, $e_i$ – jak poprzednio, $s$ – okres szeregu (dla szeregów niesezonowych równy 1). Jest to właściwie modyfikacja standardowego MASE, dostosowana do danych z sezonowością, gdzie $s$ jest długością okresu. W standardowej metodzie, w mianowniku mamy średni bezwzględny błąd jednokrokowej predykcji naiwnej, a nie sezonowej predykcji naiwnej, jak w tym przypadku. MASE będzie większe od 1 dla metod dających gorsze rezultaty niż prognozowanie metodą naiwną, a mniejsze od 1 dla lepszych.

A.2 Transformacja Boxa-Coxa

Transformacja Boxa-Coxa dla szeregu $y_t$ może być zapisana następująco: $$w_t=ln(y_t) \ \textrm{jezeli}\ \lambda=0 oraz $$w_t=(y_t^{\lambda}-1)/\lambda \ \textrm{w przeciwnym wypadku}.$$ Dla $\lambda=0$ transformacja Boxa-Coxa jest odpowiednikiem logarytmowania danych, dla pozostałych $\lambda=0$ otrzymujemy tzw. transformacje potęgowe (z dokładnością do stałej).

A.3 Dekompozycja

Model dekompozycji możemy zapisać symbolicznie: $$y_t = T_t + S_t + E_t \ \mbox{(dekompozycja addytywna)}$$ bądź $$y_t = T_t \times S_t \times E_t \ \mbox{(dekompozycja multiplikatywna)},$$ gdzie $T_t$ – trend, $S_t$ – sezonowość, $E_t$ – losowe fluktuacje.

Wybór wariantu dekompozycji zależy od danych. Model addytywny jest bardziej odpowiedni jeżeli amplituda sezonowych fluktuacji nie jest zależna od poziomu szeregu czasowego. W przeciwnym razie model multiplikatywny może okazać się bardziej prawidłowy. Alternatywą dla modelu multiplikatywnego jest zastosowanie modelu addytywnego do danych zlogarytmowanych.

A.4 Średnia krocząca (moving average)

Średnią kroczącą (średnią ruchomą) w chwili $t$ określamy jako średnią z pewnej liczby obserwacji poprzedzających i następujących (o ile oczywiście mamy odpowiednie obserwacje, co uściślimy poniżej). Średnią kroczącą rzędu $m$ wyznaczamy następująco: $$\hat{T_t} = \frac{1}{m} \sum_{i=-k}^{k} y_{t+i},$$ gdzie $m=2k+1$.

Postać dopasowanego trendu (stopień wygładzenia) zależy od wybranej wartości parametru $m$. Większy rząd $m$ prowadzi do większego wygładzenie.

A.5 Wygładzanie wykładnicze (Simple Exponential Smoothing (SES))

Niech $\alpha$ będzie parametrem wygładzającym, czyli wagą, którą przypisujemy obserwacji bezpośrednio poprzedzającej tę, której wartość chcemy wyznaczyć. Wtedy prognozowana wartość dla chwili $t$ wyraża się wzorem (zakładamy, że znamy wszystkie wartości wcześniejsze, tzn. nasza prognoza jest tylko na jeden krok naprzód): $$\hat{y_t} = \alpha y_{t-1} + \alpha (1-\alpha) y_{t-2} +\alpha (1-\alpha)^2 y_{t-3} + \cdots.$$ Każda kolejna przewidywana wartość jest oparta na poprzednich prognozach — innymi słowy, żeby w chwili $t-1$ przewidzieć, jaką wartość przyjmie szereg w chwili $t+1$, musimy najpierw wyznaczyć przewidywaną wartość $y_t$ i użyć jej do wyznaczenia $y_{t+1}$.

Dlaczego ta metoda to wygładzanie wykładnicze? Nazwa bierze się stąd, iż w jej efekcie otrzymujemy bardziej gładką kontynuację realizacji szeregu (bardziej gładką, to znaczy z mniejszymi wahaniami losowymi). Z kolei wykładnicze, bo w przedstawionym wyżej równaniu wagi przy poprzednich obserwacjach maleją wykładniczo.

A.6 Model ARIMA

Jak zostało wspomniane wcześniej, w modelu ARIMA można wyodrębnić trzy elementy:

  1. autoregresji — czyli zależności danej obserwacji $y_t$ od poprzednich obserwacji $y_{t-1}$, $y_{t-2},$ itd. To, co chcemy zrobić, to znaleźć takie współczynniki $\phi_1$, $\phi_2$, … , $\phi_p$, żeby $$y_t = \epsilon_t + \phi_1 y_{t-1} + \phi_2 y_{t-2} + \cdots + \phi_p y_{t-p}.$$ Jest to model autoregresyjny rzędu p (w skrócie: AR(p)).

  2. średniej ruchomej — czyli zależności danej obserwacji $y_t$ nie tylko od losowego błędu (zakłócenia) w danej chwili, czyli $\epsilon_t$, ale też od poprzednich zakłóceń: $\epsilon_{t-1}$, $\epsilon_{t-2},$ itd. Analogicznie do części autoregresyjnej, chcemy znaleźć takie współczynniki $\theta_1$, $\theta_2$,…,$\theta_q$, by $$y_t = \epsilon_t + \theta_1 \epsilon_{t-1} + \theta_2 \epsilon_{t-2} + \cdots + \theta_q \epsilon_{t-q}.$$ Jest to model średniej ruchomej rzędu q (w skrócie: MA(q)).

  3. związany z różnicowaniem szeregu — zakładamy, że nasz szereg być może trzeba będzie zróżnicować, to znaczy użyć do zbudowania modelu różnic między obserwacjami, zamiast oryginalnych obserwacji. Inaczej mówiąc, zamiast wyjściowego szeregu $y$ będziemy wykorzystywali szereg zróżnicowany ${y^\prime},$ stworzony w następujący sposób: $${y^\prime}_t = y_t – y_{t-k},$$ gdzie $k$ jest tzw. opóźnieniem (lag). Nie jest wykluczone, że ta operacja będzie musiała zostać wykonana więcej niż jeden raz (czasami jest potrzebne zróżnicowanie dwukrotne, rzadko większą liczbę razy). Ponieważ będziemy dopasowywać model do zróżnicowanych danych, to wyznaczając prognozy będziemy musieli wrócić do pierwotnej postaci szeregu, czyli w zastosować operację odwrotną do różnicowania.

Kompletny model ARIMA, zawierający wszystkie trzy elementy, możemy zapisać następująco: $${y^\prime}_t = c + \phi_1 {y^\prime}_{t-1} + \phi_2 {y^\prime}_{t-2} + \cdots + \phi_p {y^\prime}_{t-p} + \epsilon_t + \theta_1 \epsilon_{t-1} + \theta_2 \epsilon_{t-2} + \cdots + \theta_q \epsilon_{t-q},$$ przy czym $c$ oznacza pewną stałą.

Model ARIMA możemy uogólnić na przypadek szeregu w którym występuje sezonowość, otrzymując tzw. sezonowy model ARIMA, lub w skrócie sARIMA $(p,d,q)(P,D,Q)_s$. Model ma następującą postać: $$(1-\theta_1 B – \cdots – \theta_p B^p)(1- \Theta_1 B^s-\Theta_2 B^{2s}- \cdots – \Theta_P B^{Ps})y ^{\prime\prime} = $$ $$(1-\phi_1 B – \cdots – \phi_q B^q)(1- \Phi_1 B^s-\Phi_2 B^{2s}- \cdots -\Phi_Q B^{Qs})\epsilon_t,$$ gdzie $y^{\prime\prime}$ jest szeregiem różnicowanym najpierw $D$-krotnie z opóźnieniem $s$, a następnie $d$-krotnie z opóźnieniem 1, natomiast $B^ky_t = y_{t-k}$. Na przykład, po przekształceniu powyższego ogólnego równania, model sARIMA$(1,1,1)(1,1,1)_4$ wygląda następująco: $$y_t^{\prime\prime} – \theta_1 y_{t-1}^{\prime\prime} – \Theta_1 y_{t-4}^{\prime\prime} + \theta_1 \Theta_1 y_{t-5}^{\prime\prime} = \epsilon_t + \Theta_1\epsilon_{t-4} + \theta_1 \epsilon_{t-1} + \theta_1 \Theta_1 \epsilon_{t-5}$$

Wybór optymalnych rzędów modelu ARIMA zwykle nie jest proste. Wybór postaci modelu można przeprowadzić automatycznie (na przykład z wykorzystaniem odpowiednich narzędzi dostępnych w pakiecie R). Alternatywnie, identyfikacja optymalnego modelu może być wykonana przez analityka. Konieczne są do tego wykresy funkcji autokorelacji próbkowej (ACF) i funkcji cząstkowej autokorelacji próbkowej (PACF).

Funkcja autokorelacji (ACF) wskazuje nam jak, w zależności od odległości od siebie w czasie, związane są ze sobą obserwacje. ACF($h$), gdzie zwyczajowo $h=1,2,…,N/4$ (a $N$ to długość szeregu czasowego) jest współczynnikiem korelacji między $y_t$$y_{t-h} $, $t=h+1,…,N$. Może się jednak zdarzyć, że korelacja między $y_t$$y_{t-1}$ powoduje, że $y_t$$y_{t-2}$ są również skorelowane, ale nie jest to związane z postacią modelu, z którego pochodzi dany szereg. Aby uniknąć takiej sytuacji, rozpatruje się również funkcję cząstkowej autokorelacji (PACF). Mierzy ona korelację między $y_t$$y_{t-h}$ z pominięciem wpływu obserwacji pośrednich, tzn. $y_{t-1},…,y_{t-h+1}$. Dokładną informację, jak za pomocą takich wykresów wybrać rząd modelu, możemy znaleźć np. w książce Forecasting: principles and practice.

A.7 ETS (ExponenTial Smoothing, ErrorTrendSeasonality )

Model ETS, jest tzw. modelem przestrzeni stanów (state space), którego podstawą jest równanie różnicowe opisujące ewolucję układu. Przyjmowane są także odpowiednie założenia dotyczące rozkładu czynnika losowego (zakłócenia). Przedstawimy tutaj najprostszy przypadek modelu ETS odpowiadający metodzie prostego wygładzania wykładniczego (SES). Alternatywną formą zapisu SES jest $$l_t = l_{t-1} + \alpha(y_t – l_{t-1}) = l_{t-1} + \alpha \epsilon_t,$$ gdzie $\epsilon_t = y_t – l_{t-1},$$l_{t-1} = \hat y_t$. O błędach $\epsilon_t$ zakładamy, iż są niezależnymi zmiennymi losowymi z jednakowego symetrycznego rozkładu normalnego. To założenie, wraz z powyższym równaniem ewolucji, składa się na pełny model statystyczny. Dzięki niemu jesteśmy w stanie obliczyć funkcję wiarogodności oraz skonstruować przedziały ufności.

Inne warianty wygładzania wykładniczego (np. algorytmy Holta i Holta-Wintersa) również mają odpowiadające im modele ETS. Szczegółowe informacje można znaleźć np. w książce Forecasting: principles and practice.

A.8 Regresja harmoniczna

Model regresji harmonicznej jest oparty na założeniu, że możemy przybliżyć szereg czasowy pewną sumą funkcji okresowych, takich jak sinus i cosinus, o różnych okresach. Sumę taką możemy zapisać wzorem: $$y_t = a + \sum_{i=1}^{n/2} \left[ \alpha_i sin(2\pi i t/m) + \beta_i cos(2\pi i t/m) \right] + \epsilon_t,$$ gdzie $a$, $\alpha_i$$\beta_i$ są pewnymi współczynnikami, $\lambda_i$ to kolejne częstotliwości (odpowiadające kolejnym okresom), a $\epsilon_t$ oznacza szereg reszt, który możemy modelować przy użyciu ARIMA.

Dzięki użyciu funkcji okresowych możemy zobaczyć, czy wartość szeregu w danej chwili nie jest okresowo zależna od wcześniejszych obserwacji, nawet od kilku — na przykład, na dzisiejszą wartość szeregu wpływa wartość zarówno sprzed roku, jak i sprzed 3 miesięcy. Te zależności możemy wykryć analizując współczynniki $\alpha$$\beta$ dla poszczególnych składowych harmonicznych (im większy współczynnik, tym większy wpływ danej funkcji okresowej).

Spróbuj ponownie