W poprzednim poście było trochę wprowadzenia do szeregów czasowych, trochę o stacjonarności, sezonowości i trendzie. A teraz zajmiemy się kilkoma metodami takimi jak ARIMA czy SARIMA używanymi do prognozowania szeregów czasowych.
ARIMA oznacza Auto Regressive Integrated Moving Average. Ten model możemy użyć do niesezonowych danych. Model ARIMA składa się z dwóch części: model autoregresyjny (AR) oraz model średniej ruchomej (moving average MA) i ma trzy zmienne do uwzględnienia:
ARIMA(p, d, q)
p = Okresy zwłoki (jeśli P = 3, to użyjemy trzech poprzednich okresów z naszych szeregów czasowych w części autoregresyjnej obliczeń) P pomaga dostosować dopasowaną linię do prognozowania serii
d = W modelu ARIMA przekształcamy szereg czasowy w stacjonarny (seria bez trendu lub sezonowości) za pomocą różnicowania. D oznacza liczbę transformacji różnicowych wymaganych przez szereg czasowy, aby uzyskać stacjonarne.
q = Ta zmienna oznacza opóźnienie komponentu błędu, gdzie komponent błędu jest częścią szeregu czasowego nie wyjaśnionego przez trend lub sezonowość
ARIMA sprawdzi jeśli nie mamy elementu sezonowego, jeśli jednak istnieje mamy inny model – SARIMA. SARIMA jest podobny do modeli ARIMA, musimy tylko dodać kilka parametrów, aby uwzględnić sezonowość:
SARIMA (p, d, q) (P, D, Q) m
p – liczba autoregresji
d – stopień różnicowania
q – liczba średnich ruchomych terminów
m – odnosi się do liczby okresów w każdym sezonie
(P, D, Q) – reprezentuje (p, d, q) dla sezonowej części szeregu czasowego
Zapis dla modelu obejmuje określenie kolejności modeli AR (p), I (d) i MA (q) jako parametrów funkcji ARIMA i AR (P), I (D), MA (Q) i m parametry na poziomie sezonowym.
Przykład:
SARIMA (3,1,0) (1,1,0) 12
- Parametr m wpływa na parametry P, D i Q. Na przykład m=12 w przypadku danych miesięcznych sugeruje roczny cykl sezonowy.
- P = 1 użyta jest pierwsza sezonowa obserwacja w modelu, np. t- (m 1) lub t-12. P = 2, użyte są dwie ostatnie sezonowe obserwacje t- (m 1), t- (m * 2).
- D = 1 wyliczy różnicę sezonową pierwszego rzędu, a Q = 1 użyje błędy pierwszego rzędu w modelu (np. Średniej ruchomej).
import pandas as pd
import numpy as np
import matplotlib.pylab as plt
%matplotlib inline
from matplotlib.pylab import rcParams
rcParams['figure.figsize'] = 20, 6
data = pd.read_csv(r'C:\Users\VFW863\Documents\drugstore.csv')
data['date'] = pd.to_datetime(data['date'])
data = data.set_index('date'); data.tail(3)
plt.plot(data)
Czy dane wykazują tendencje sezonowe? Jest to ważne przy podejmowaniu decyzji, który typ modelu zastosować. Jeśli nie ma sezonowego trendu możemy użyć zwykłego modelu ARIMA. Nasz model jednak zawiera elementy sezonowości więc użyjemy model SARIMA.
Zanim przejdziemy do samego modelu podzielimy dane na treningowe i testowwe oraz przekształcimy dane aby były stacjonarne, czyli użyjemy różnicowania:
train = 0.75
split = round(len(data)*train)
training, testing = data[0:split], data[split:]
training_log = np.log(training)
training_roz = training_log.diff(periods=1).values[1:]
plt.plot(training_roz)
Aby określić parametry modelu możemy użyć dwóch wykresów: autokorelacji (ACF) i częściowej autokorelacji (PACF).
Autokorelacja odnosi się do korelacji szeregu czasowego z jego przeszłymi wartościami. W ACF współczynnik korelacji jest na osi X, podczas gdy liczba opóźnień (lag) jest pokazana na osi Y (-1:1). Wykres funkcji autokorelacji informuje, w jaki sposób dane szeregi czasowe są ze sobą skorelowane. Zwykle w modelu ARIMA używamy elementu AR, elementu MA lub oba. Wykres ACF pomaga nam zdecydować, który z tych terminów będziemy używać w naszych szeregach czasowych.
- Jeśli istnieje autokorelacja dodatnia w przypadku opóźnienia 1, wówczas używamy modelu AR
- Jeśli istnieje negatywna autokorelacja w fazie 1, wówczas używamy modelu MA
Różnice sezonowe uwzględniają np pory roku i różnice między bieżącą wartością a jej wartością w poprzednim sezonie, np .: Różnica za miesiąc może być wartością w maju 2018 r. – wartość w maju 2017 r.
W czysto sezonowym modelu AR, ACF rozpada się powoli, a PACF odcina się do zera
Modele AR są stosowane, gdy autokorelacja sezonowa jest dodatnia
W czysto sezonowym modelu MA, ACF odcina się do zera i na odwrót
Modele MA są stosowane, gdy autokorelacja sezonowa jest ujemna
Po wykreśleniu wykresu ACF przechodzimy do wykresu funkcji częściowej autokorelacji (PACF). Częściowa autokorelacja jest podsumowaniem związku między obserwacją w szeregu czasowym a obserwacjami w etapach poprzedzających czas, z usuniętymi relacjami obserwowanych interwencji. PACF opisuje tylko bezpośredni związek między obserwacją a jej opóźnieniem.
Użyjemy parametry AR w modelu, gdy:
- Wykresy ACF pokazują autokorelację rozpadającą się w kierunku zera
- Wykres PACF odcina się szybko w kierunku zera
- ACF z serii stacjonarnej pokazuje pozytywny w lag-1
Użyjemy parametry MA w modelu, gdy:
- Negatywnie Autokorelowane w Lag – 1
- ACF, który spada po kilku opóźnieniach
- PACF zmniejsza się stopniowo
Przedziały ufności są rysowane jako stożek. Możemy ograniczyć liczbę opóźnień na osi X do 50, aby działka była łatwiejsza do odczytania. Aby oszacować wartość p potrzebujemy wykres ACF a do wartości q wykres PACF
import numpy as np
import statsmodels.api as sm
from statsmodels.tsa.stattools import acf, pacf
lag_acf = acf(training_roz, nlags=50)
lag_pacf = pacf(training_roz, nlags=50, method='ols')
fig = plt.figure(figsize=(15,5))
ax1 = fig.add_subplot(121)
plt.axhline(y=-1.96/np.sqrt(len(training)), linestyle='--', color='red')
plt.axhline(y=1.96/np.sqrt(len(training)), linestyle='--', color='red')
plt.stem(lag_acf)
ax2 = fig.add_subplot(122)
plt.axhline(y=-1.96/np.sqrt(len(training)), linestyle='--', color='red')
plt.axhline(y=1.96/np.sqrt(len(training)), linestyle='--', color='red')
plt.stem(lag_pacf)
Z powyższych wykresów możemy otrzymać p (element AR) i q (element MA) do naszego modelu. Gdy prześledzimy wartości na obu wykresach możemy mniej więcej oszacować p i q: w naszym przykładzie dobrym początkiem będzie p=1 i q=1 (są to pierwsze wartości przekraczające czerwone linie). Parametr d zależy od różnicowania naszych danych, jeśli dane nie są stacjonarne – to d musi być co najmniej 1. Analizując wykresy dalej zauważymy, że mamy dużą wartość na 12 (oś X) co sugeruje sezon (miesiąc), m=12, jest to wartość dodatnia więc P=1 a Q=0. Dane nie są stałe więc przyjmiemy D=0.
Parametry do naszego modelu to:
SARIMA(1,1,1)(1,0,0)[12]
import warnings
warnings.filterwarnings("ignore")
from statsmodels.tsa.statespace.sarimax import SARIMAX
model = SARIMAX(training, order=(1,1,1), seasonal_order=(1,0,0,12), enforce_stationarity = False, enforce_invertibility = False)
model_fit = model.fit(disp=False)
Teraz możemy prognozować sprzedaż w aptece na kolejne lata:
K= len(testing)
prognoza = model_fit.forecast(K)
from sklearn.metrics import mean_squared_error
plt.figure(figsize=(15,5))
plt.plot(prognoza, 'red')
plt.plot(data, 'b')
plt.title ('RMSE: %.2f'% mean_squared_error(testing, prognoza))
plt.axvline(x=data.index[split], color = 'black')
Nasza prognoza pokrywa się z danymi przez pierwsze około 2 lata (2012-2014) a potem możemy zauważyć różnice a mianowicie prognozowana sprzedaż jest mniejsza od rzeczywistej. Ta różnica może być spowodowana m.in. doborem parametrów modelu SARIMAX.
Nasza kombinacja parametrów wcale nie musi być najbardziej optymalna, jeśli chcemy przetestować inne wartości możemy ręcznie zmieniać parametry lub użyć grid search.
Generalnie możemy podzielić modelowanie szeregów czasowych na etapy:
1 – Sprawdź stacjonarność: Jeśli szereg czasowy ma składową trendu lub sezonowość, musi on być stacjonarny, zanim będziemy mogli użyć ARIMA/SARIMA do prognozowania.
2 – Różnica: Jeśli szeregi czasowe nie są stacjonarne, muszą być stacjonarne poprzez różnicowanie. Weź pierwszą różnicę, a następnie sprawdź stacjonarność. Zrób tyle różnic, ile potrzeba. Upewnij się, że sprawdzasz także różnice sezonowe.
3 – Podziel dane na testowe i treningowe.
4 – Wybierz elementy AR i MA: Użyj ACF i PACF, aby zdecydować, czy uwzględnić AR, MA lub oba.
5 – Zbuduj model i ustaw liczbę okresów do prognozy na N
6 – Sprawdź poprawność modelu: porównaj przewidywane wartości z rzeczywistymi
