W niniejszym artykule, postaram się przybliżyć podstawowe pojęcia w przypadku jednowymiarowej analizie szeregów czasowych. Dodatkowo zaimplementuję prosty model Autoregresji AR(p) w Pythonie wraz z oceną i interpretacją wyników. Ten cykl artykułów będzie skupiony bliżej na ekonometrii niż uczeniu maszynowemu czy głębokiemu, jednakże i takie modele prawdopodobnie pojawią się w dalszych artykułach. Zapraszam!
Małe duże różnice
Czym się różnią modele jednowymiarowej analizy szeregów czasowych od analizy przekrojowej. Wyróżnia się dwie przyczyny gdzie powstają problemy z analizą szeregów:
- Jeden szereg czasowy może wpłynąć na drugi z opóźnieniem
- Jeżeli mamy do czynienia z niestacjonarnymi zmiennymi powstaje problem regresji pozornej.
Co do pierwszego podpunktu spróbuje wyjaśnić na przykładzie, czemu analitycy są bardzo ostrożni w doborze zmiennych w modelach szeregów. Odnosi się to bezpośrednio do faktu, że niektóre zmienne działają na zmienną zależną z opóźnieniem co w przypadku predykcji szeregu czasowego jest silnie niepożądaną właściwością. Wyobraźmy sobie prosty przykład, badamy stopę inflacji i dodatkowo naszą zmienną objaśniającą jest wartość stóp procentowych. Bank centralny podnosi stopy żeby hamować inflację, zanim jednak skutki zaczną być zauważalne dla gospodarki oraz poziomu stopy inflacji, może minąć nawet rok. To jest przykład makroekonomiczny ale również w mniejszej skali wpływ niektórych zmiennych oddziałuje z opóźnieniem. Dlatego też skupimy się na badaniu jednowymiarowych szeregów czasowych.
Co do drugiego podpunktu pojęcie stacjonarności zostanie wyjaśnione w dalszej części artykułu, na razie możemy postrzegać stacjonarność jako coś pożądanego w modelowaniu szeregów czasowych. Regresja pozorna jest to przypadek w którym zmienna objaśniająca i zmienne objaśniane mają taki sam trend oznacza to, że model który wskazuje że zmienne (predyktory) objaśniają zmienną celu jest błędny a te wskazania to tylko i wyłącznie przypadek dlatego nazwa pozorna.
Notacja szeregów czasowych
Najważniejszym pojęciem dla zrozumienia ekonometrycznej analizy szeregów czasowych jest zrozumienie opóźnień zmiennej. Jednowymiarowy szereg czasowy może zostać przedstawiony w następujący sposób:
Opóźnienie jest to wskazanie o ile okresów czasu (t) cofamy się, żeby uwzględnić wartość zmiennej np. y. Przykładowo jeżeli chcemy dodać zmienną opóźnioną o jeden okres oznaczymy ją jako: Yt-1 a więc gdy nasza zmienna Y będzie wskazywała na wartość y4 to zmienna opóźniona Yt-1 wskaże wartość y5. Oczywiście można opóźniać zmienne o więcej okresów itd. Idea jest taka sama. Należy pamiętać https://becejprevoz.com/buchstabe-a/index.html , że uwzględniając kolejne opóźnienia zmniejsza nam się liczba obserwacji wchodzących do modelu. Jeżeli dodajemy zmienną opóźnioną o dwa okresy to do modelu wejdzie T-2 obserwacje.
Zły trend
W przypadku szeregów czasowych to czy zmienna jest niestacjonarna bądź stacjonarna jest bezpośrednio związane z tym czy występuje w szeregu trend czy nie. Bardziej formalne wyjaśnienie stacjonarności znajdzie się jeszcze niżej :).
Trend możemy zauważyć na poniższym przykładzie gdzie na podstawie zbioru danych 30 indeksów giełdowych z Kaggle, stworzyłem wykres indeksu akcji Google:
import pandas as pd
import matplotlib.pyplot as plt
%matplotlib inline
plt.style.use('fivethirtyeight')
#google
google = pd.read_csv(r'GOOGL_2006-01-01_to_2018-01-01.csv', index_col='Date', parse_dates=['Date'])
#wykres
google['High'].plot(figsize=(12,6));
Bardzo łatwo zauważyć trend wzrostowy. W wielu dziedzinach można zaobserwować tworzenie się trendów np. poziomy cen, produkcja, konsumpcja czy właśnie indeksy giełdowe. Skoro trendy są niepożądane w jaki sposób możemy się ich pozbyć?
Najłatwiej tego dokonać poprzez różnicowanie. To znaczy, że jeżeli Yt to nasz szereg czasowy to ∆Yt = Yt – Yt-1 jest to pierwsza różnica. Dokonajmy więc takiego przekształcenia:
#diff() od pandas - pierwsza roznica google['High'].diff().plot(figsize=(12,6));
Widać bardzo dużą różnicę pomiędzy pierwszym a drugim wykresem, na drugim nie obserwujemy już silnego trendu. Dodatkowo interpretacja takiego wykresu się zmienia. Otóż pierwsze różnice w naszym przypadku to będzie zmiana wartości ceny najwyższej w stosunku do dnia poprzedniego. Należy pamiętać o tej interpretacji.
Okej, okej… Wspomniałem, że trendy są niepożądane, ale za bardzo nie wyjaśniłem dlaczego. Chodzi tutaj o korelacje pomiędzy obserwacjami. W przypadku wielu szeregów gdy mamy do czynienia z trendem to ze względu na fakt, że wartość indeksu giełdowego z dzisiaj jest silnie skorelowana z wartością indeksu z zeszłego roku itd. Gdybyśmy policzyli korelację dla takich wartości otrzymalibyśmy wyniki bliskie 0.99 wsp. korelacji. Natomiast jeżeli policzymy korelacje pomiędzy obserwacjami w szeregu z pierwszych różnić ta wartość będzie bardzo niska. Zmiany nie są ze sobą skorelowane (zazwyczaj choć nie zawsze). To że składnik jest skorelowany sam ze sobą w przypadku szeregów określa się mianem autokorelacji. Trend oraz autokorelacja są bardzo ważne w kontekście modelowania szeregów czasowych.
Autokorelacja
Jak ustaliliśmy wcześniej autokorelacja jest to korelacja między zmienną i jej opóźnieniem. Funkcje korelacji możemy kalkulować dla wielu stopni opóźnień badanej zmiennej. Przykładowo możemy obliczyć jak wygląda korelacja dla poszczególnych opóźnień dla naszego zbioru danych google. Zarówno dla zwykłych opóźnień jak i dla różnicowanego szeregu. Postanowiłem zbadać opóźnienia od 1 do 12 stopnia:
opoznienie = np.arange(1,13)
google_normal = []
google_diff = []
for i in opoznienie:
google_normal.append(round(np.corrcoef(google['High'][:-i],google['High'][i:])[0][1],4))
for i in opoznienie:
google_diff.append(round(np.corrcoef(google['High'].diff().dropna()[:-i],google['High'].diff().dropna()[i:])[0][1],4))
#tabelka z opóźnieniami i autokorelacją
pd.DataFrame({'opoznienie': opoznienie,
'Indeks Korelacja': google_normal,
'Indeks Opozniony Korelacja': google_diff})
Przykład ten daje nam mnóstwo informacji. Przede wszystkim potwierdza fakt, że różnicowanie pozwala na usunięcie autokorelacji co jest dla nas pożądane. Po drugie to co zaobserwowaliśmy na wykresach również się zgadza. Funkcja autokorelacji jest przydatnym narzędziem i pozwala wyciągać pewne wnioski na temat szeregu:
- Y jest silnie skorelowane w czasie, wartości indeksu nawet do 12 opóźnienia wskazują na bardzo silną korelacje. Natomiast zmiany wartości indeksu praktycznie nie są skorelowane wcale.
- Na podstawie wartości indeksu Y moglibyśmy z dużym prawdopodobieństwem zgadywać przyszłe wartości indeksu. W przypadku zmian wartości nie możemy tego uczynić.
- Y niezróżnicowane jest niestacjonarnym szeregiem czasowym, natomiast Y zróżnicowane jest stacjonarnym.
W tym momencie warto zauważyć, że wysokie wartości funkcji autokorelacji są charakterystyczne dla niestacjonarnych szeregów czasowych. Dodatkowo przypomnę, że niestacjonarność jest niepożądana.
Model AutoRegresji – AR
Model autoregresji AR(i) jest najprostszym modelem szeregów czasowych. Nazwa odnosi się do tego, że zmienną objaśniającą jest również zmienna objaśniana tylko z opóźnieniem/opóźnieniami. W przypadku tego modelu AR(i) – i oznacza liczbę opóźnień, które wprowadzamy do modelu. I tak model AR(1) będzie prezentował się w sposób następujący:
Sam model prezentuje się niemal jak zwykła regresja prosta. Z tym, że jako zmienną objaśnianą traktujemy opóźnienia zmiennej celu. Bardzo ważnym składnikiem jest wartość \rho w tym modelu (AR(1)) jest ona bezpośrednio powiązana z funkcją autokorelacji i niestacjonarnością. Model AR(1) nazywany jest często procesem autoregresyjnym dlatego jeżeli taka terminologia wymieniana jest w literaturze będziemy wiedzieć do jakiego zjawiska się odnosi.
Aby zobaczyć na przykładzie w jaki sposób \rho wiąże się z autokorelacją. Wygeneruję 3 trzy różne szeregi w procesie autoregresji z 3 różnymi wartościami parametru rho.
Symulacja procesu AR(1):
import statsmodels.api as sm
#procesy AR
ar1 = np.array([1, -0.85]) #-0.85 odpowiada rho=0.85
ma1 = np.array([1]) #niwelowanie procesu MA
AR1 = sm.tsa.ArmaProcess(ar1, ma1)
sim1 = AR1.generate_sample(nsample=1000)
ar2 = np.array([1, 0])
AR2 = sm.tsa.ArmaProcess(ar2, ma1)
sim2 = AR2.generate_sample(nsample=1000)
ar3 = np.array([1, -1])
AR3 = sm.tsa.ArmaProcess(ar3, ma1)
sim3 = AR3.generate_sample(nsample=1000)
#wykresy
fig = plt.figure(figsize=(20,15))
(ax1, ax2, ax3) = fig.subplots(3,1)
ax1.plot(sim1)
ax1.set_title("Rho = 0.85", fontsize=26)
ax2.plot(sim2)
ax2.set_title("Rho = 0", fontsize=26)
ax3.plot(sim3)
ax3.set_title("Rho = 1", fontsize=26)
fig.tight_layout();Wykresy:
Wszystkie szeregi mają tą samą wartość alphy czyli 1. Co jest charakterystyczne dla tych procesów? Wykres procesu AR z parametrem rho = 0, pokazuje nam losowe fluktuacje wokół średniej w wysokości 1 czyli parametru alpha. Wykres ten przypomina nasz wcześniejszy obraz na którym przedstawiliśmy wykres zmian wartości indeksu (zróżnicowany). Wykres procesu AR z parametrem rho = 1 przedstawia silny trend, jest to szereg niestacjonarny. Natomiast wykres procesu AR rho = 0.85 prezentuje coś pomiędzy losowymi fluktuacjami a ukazanym trendem.
W przypadku szeregów czasowych dla procesu AR mówimy, że szereg jest stacjonarny jeżeli |\rho| < 1 i niestacjonarny jeżeli \rho = 0. Raczej nie bada się sytuacji gdy |\rho| > 1 ponieważ oznacza to sytuacje kiedy szereg czasowy eksploduje w czasie. Są to pewne nietypowe zdarzenia jak np. hiperinflacja itd.
W literaturze sytuacje w której \rho = 1 w procesie AR, określa się mianem występowania pierwiastka jednostkowego. Jeżeli pierwiastek jednostkowy występuje szereg jest niestacjonarny. Dodatkowo mówi nam to że korelacje są bliskie 1 i nie będą spadać wraz z wydłużeniem okresu opóźnienia. Wiemy również, że szereg z takim pierwiastkiem wykazuje trend. I jeżeli taki pierwiastek występuje, to zróżnicowany szereg tego pierwiastka posiadać nie będzie. Są to szeregi stacjonarne w pierwszych różnicach. Chociaż zdarzają się sytuacje, gdzie pierwsze różnice nadal są niestacjonarne ale ten temat poruszmy w dalszych artykułach.
Trochę bardziej formalnie. Model AR(1) który posiada pierwiastek jednostkowy jest również określany jako model błądzenia losowego (Random walk):
(bardziej formalnie powyższy model określa się mianem błądzenia losowego z dryfem jeżeli alpha jest różna od 0).
Najlepiej zrozumieć model błądzenia losowego na przykładzie. Uważa się, że taki model dobrze odwzorowuje zjawiska takie jak ceny na giełdzie. Cena akcji dzisiaj to cena akcji z wczoraj plus szum.
Stacjonarność formalnie
Zanim przejdę do omawiania diagnostyki modelu AR omówię formalnie wyczekiwaną niecierpliwie stacjonarność. Są to pewne warunki, które muszą zostać spełnione, żeby modele szeregów czasowych mogły działać.
Szereg czasowy Y_t jest stacjonarny jeżeli:
E(Y_t)– (wartość oczekiwana) – jest taka sama dla wszystkich wartości t,- wariancja
var(Y_t)jest skończona i taka sama dla wszystkich wartości t, - Kowariancja
cov(Y_t, Y_{t-s}), zależy tylko od s, a nie t.
Te warunki implikują to, że model jest niezmienny w czasie. Niestacjonarność określa po prostu fakt nie spełniania tych warunków.
Wykrywanie stacjonarności – dekompozycja i testy
Teraz, wiedząc dokładnie czym jest stacjonarność możemy zastosować kilka metod pozwalających na jej wykrywanie. Wśród takich metod najpopularniejszymi są dekompozycja szeregu oraz test Dickeya-Fullera.
Dekompozycja szeregu jak sama nazwa wskazuje jest rozłożeniem szeregu na 3 główne składniki, z których zbudowany jest szereg. Te czynniki to:
- Trend – autokorelacja szeregu czasowego, ciągłe nachylenia w czasie szeregu,
- Sezonowość – pewne powtarzające się wzorce w szeregu czasowym (np. szereg przypomina funkcje sinusoidalną),
- Szum/Nieregularność – jeżeli szum jest stacjonarny, taki szereg da się prognozować.
Dekompozycja służy wychwyceniu czy wszystkie te składniki istnieją w szeregu i pozwala na podjęcie pewnych kroków w celu eliminacji problemów z nimi związanych. Przeprowadźmy dekompozycję szeregu czasowego google oraz zróżnicowanego szeregu google:
plt.rcParams['figure.figsize'] = 11, 9 decomposed_google_volume = sm.tsa.seasonal_decompose(google["High"],period=360) # roczny okres, jak miesieczny to 30 itd decomposed_google_volume.plot() plt.show();
Otrzymujemy nasz oryginalny szereg, oraz 3 wymienione wyżej składniki. Warto zwrócić uwagę, że istotnie trend istnieje w takim szeregu widać ciągłą funkcję o dodatnim nachyleniu. Dodatkowo warto zauważyć, że występuje sezonowość widzimy pewien powtarzający się schemat mniej więcej rok do roku. Nasz szum również przyjmuje pewien schemat co nie jest pozytywne. Zobaczmy teraz dekompozycje szeregu pierwszych różnic:
plt.rcParams['figure.figsize'] = 11, 9 decomposed_google_volume = sm.tsa.seasonal_decompose(google["High"].diff()[1:],period=360) # roczny okres, jak miesieczny to 30 itd decomposed_google_volume.plot() plt.show();
Na powyższych wykresach jest zdecydowana różnica. W przypadku trendu nie jest już tak widoczny, jednakże w mojej prywatnej ocenie nie można stwierdzić ze stuprocentową pewnością, że trendu nie ma. Trzeba będzie posłużyć się testem Dickeya-Fullera. W przypadku sezonowości nadal widać powtarzające się schematy, oznacza to że pierwsze różnice nie były w stanie usunąć sezonowości z szeregu. Co do szumu, widać ogromną różnice. Nie ma tutaj prawie żadnego schematu. Przejdźmy dalej do weryfikacji stacjonarności szeregu testem Dickeya-Fullera.
Hipoteza zero testu Dickeya-Fullera zakłada, że w badanym szeregu występuje pierwiastek jednostkowy, czyli szereg jest niestacjonarny. Przeprowadźmy więc testowanie na naszych szeregach:
# Augmented Dickey-Fuller test
from statsmodels.tsa.stattools import adfuller
adf = adfuller(google["High"])
print("p-value Google: {}".format(float(adf[1])))
adf = adfuller(google["High"].diff()[1:])
print("p-value Google Pierwszych różnic: {}".format(float(adf[1])))
Jak widać na powyższym obrazku dla szeregu Googla nie ma podstaw do odrzucenia hipotezy zero poziomie istotności 0.05, stąd wniosek, że szereg jest niestacjonarny. W przypadku szeregu Googla pierwszych różnic, odrzucamy hipotezę zerową, szereg jest stacjonarny. Potwierdza to fakt, że różnicowanie jest często skuteczną metodą na stacjonarność szeregu.
Estymacja, diagnostyka, predykcja modelu AR(1)
Estymacja oraz diagnostyka modelu AR(1) nie różni się od estymacji i diagnostyki modelu regresji prostej. Otrzymujemy takie informacje jak istotność parametrów, kryteria informacyjne itd. To co najczęściej nas interesuje to predykcja, którą tutaj przeprowadzimy. Oczywiście pamiętamy, że do estymacji używamy szeregu stacjonarnego czyli szeregu Google pierwszych różnic. Przeprowadźmy estymację z wykorzystaniem modułu statmodels:
# Przewidywanie najwyższych wartości indeksu Google from statsmodels.tsa.arima_model import ARMA model = ARMA(google["High"].diff().iloc[1:].values, order=(1,0)) res = model.fit() res.summary()
Jak wspomniałem wcześniej, otrzymujemy p-wartość istotności oszacowań parametrów, oraz przedziały ufności. Ar.L1.y – oznacza pierwsze opóźnienie zmiennej y w procesie AR. Model ARMA(1,0) oznacza, że wykorzystaliśmy modele ARMA jednakże zdefiniowaliśmy tylko parametry dla modelu AR żeby nie estymować MA (szczegóły tego modelu pojawią się w dalszych artykułach). Otrzymaliśmy również kryteria informacyjne AIC, BIC, HQIC, wykorzystuje się je głównie do porównania modeli między sobą. Im niższe wartości tym model lepszy. Dokonajmy predykcji i policzmy błąd średniokwadratowy:
#metryka
from sklearn.metrics import mean_squared_error
import math
rmse = math.sqrt(mean_squared_error(google['High'].diff()[900:1000], res.predict(start=900,end=999)))
print("błąd średniokwadratowy {:.2f}".format(rmse))
Nie jest źle, pamiętajmy jednak, że gdy pracujemy na szeregu pierwszych różnic interpretacja parametru oraz wartości błędu się zmieniają. Zobaczmy jak prezentuje się predykcja na wykresie. Statmodels zapewnia bardzo proste funkcje dokonujące predykcji wraz z wykresami:
res.plot_predict(start=900, end=1010) plt.show()
Widać, tutaj różnicę pomiędzy predykcją a prawdziwymi wartościami. Co ciekawe kierunek predykcji najczęściej się zgadzał, jednakże wartości już nie. Dodatkowo możemy zobaczyć pewne opóźnienia w predykcji. Należy pamiętać, że używamy najprostszego modelu jaki istnieje do predykcji szeregów czasowych. Jak na takie warunki nasza predykcja była całkiem niezła.
W artykule opisywałem podstawy estymacji szeregów czasowych, wprowadzenie do podstawowych pojęć oraz jak wykorzystać Pythona do ich estymacji i pracy z nimi. Mam nadzieję, że nie pojawiło się wiele niejasności. Oczywiście to co opisałem, to jedynie wierzchołek góry lodowej dotyczącej modelowania szeregów czasowych. Mam nadzieję, że nie zbraknie mi motywacji na pisanie dalszych artykułów na ten temat i wdrażania nowych bardziej zaawansowanych modeli oraz technik modelowania. Stay tuned!
Artykuł napisano na podstawie:
- G. Koop, Wprowadzenie do ekonometrii
- https://el.us.edu.pl/ekonofizyka/index.php/Analiza_Szereg%C3%B3w_Czasowych/Dekompozycja_szeregu_czasowego#Klasyczna_dekompozycja_sygna.C5.82u
- https://www.kaggle.com/thebrownviking20/everything-you-can-do-with-a-time-series