Auto ARIMA Pythonban

A mai bejegyzésben egy idősorokon gyakran használt modellt fogunk megismerni: az AutoRegressive Integrated Moving Average-t (ARIMA)

ARIMA

Az ARIMA három részből épül fel:

AR: Autoregression — ebben a részben autoregressziót végzünk az idősorunkon.
I: Integrated — deriválás ellentéte ebben az esetben
MA: Moving Average — ebben a részben egy mozgó átlaggal számolt modell hibáját használjuk a következő megfigyelés előrejelzésére.

Most nézzük ezeket a részeket egyesével.

AR

Ebben a részben azt becsüljük, hogy a megfigyelés mostani, $Y_t$ értéke hogyan függ az előző időszakok $Y_{t-1},Y_{t-2}, \dots Y_{t-p}$ értékeitől. Ez így írható fel:

$Y_t = c + \sum_{i=1}^p \alpha_i Y_{t-i} + \epsilon_t$

Ahol:

$Y_t$ — az időpillana, amit becsülni szeretnénk
$c$ — az eltolás, ami egy konstans
$p$ — milyen mesze megyünk vissza az időben
$Y_i$ — i időpillanat megfigyelése
$\alpha_i$ — az i időpillanat súlya
$\epsilon$ — zaj

Érdemes megjegyezni, hogy az AR modell nem feltétlen időfüggetlen (angolul: stationary).

MA

A Mozgó átlag első lépésében veszünk egy Mozgó átlag modellt, és megnézzük mekkora a modell ennek a modellnek a hibája. Vagyis a múltban elkövetett becslési hibák segítségével próbáljuk becsülni a következő megfigyelést. Ez így írható fel:

$Y_t = \mu + \epsilon_t + \sum_{i=1}^q \beta_i \epsilon_{t-i}$

Ahol:

$Y_t$ — az időpillanat, amit becsülni szeretnénk
$q$ — milyen mesze megyünk vissza az időben
$\epsilon_t$ –egy standard normál eloszlás.
$\epsilon_i$ — i időpillanatra adott becslésünk hibája
$\beta_i$ — az i időpillanat súlya
$\mu$ — az idősor átlaga, általában 0-nak tekintik

Talán ez így egy kicsit zavaros, tehát nézzünk két példát egy idősor MA(1) és MA(2) modelljét:

Az MA(1) modell a fentiek alapján: $Y_t = 20 + \epsilon_t - 0,8 \epsilon_{t-1}$ , az MA(2) pedig: $Y_t = 0 + \epsilon_t - \epsilon_{t-1} + 0,8 \epsilon_{t-2}$

Vagyis becsült értékét a t időpontban befolyásolja a jelenlegi hibatag és a múltbeli hibatagok súlyozott kombinációja.

AR és MA kapcsolata

Ami érdekes, hogy egy AR(p) modellt felírhatunk úgy, mint egy MA(∞) modellt. Tegyük fel, hogy az AR model eltolása ( $c$ -je) 0. Ebben az esetben az AR(1) modell:

$Y_t = \alpha_1 Y_{t-1} + \epsilon_t$

Viszont a $Y_{t-1}$ felírható a t_1 időpillanatra készített AR modellel: $Y_{t-1} = \alpha_2 Y_{t-2} + \epsilon_{t-1}$ . Helyettesítsük ezt be:

$Y_t = \alpha_1 ( \alpha_2 Y_{t-2} + \epsilon_{t-1}) + \epsilon_t$

Most végezzük el a szorzást:

$Y_t = \alpha_1 \alpha_2 Y_{t-2} + \alpha_1 \epsilon_{t-1} + \epsilon_t$

Viszont vegyük észre, hogy az előzőekhez hasonlóan a $Y_{t-2}$ felírható a t_2 időpillanatra készített AR modellel. És ezt végtelenségig lehet folytatni. Vegyük azt is észre, hogy a súly kisebb mint 1, tehát ahogy az $\alpha$ értékeket hatványozzuk egyre kisebb értékeket fogunk kapni. Az egész vége a következő kifejezés lesz:

$Y_t = \epsilon_t + \alpha_1 \epsilon_{t-1} + \alpha_1 \alpha_2 \epsilon_{t-2} + \alpha_1 \alpha_2 \alpha_3 \epsilon_{t-3} + \dots$

Ami pedig egy MA(∞) modell, aminél minél távolabb megyünk vissza az időben annál kisebb a megfigyelés súlya.

Meg lehet ezt fordítva is csinálni? Igen. MA(1) modellt felírhatjuk úgy is mint AR(∞).

Ugye mivel mind a két modell ugyanazt a $Y_t$ -t jelzi előre, ezért ha feltesszük, hogy a $\mu=0$ , és a $\sum \alpha + \sum \beta = 1$ , akkor:

$Y_t = c + \sum_{i=1}^p \alpha_i Y_{t-i} + \epsilon_t + \sum_{i=1}^q \beta_i \epsilon_{t-i}$

Nem meglepő módon a fenti modellt ARMA(p,q)-nak nevezik. Már csak egy betűre vagyunk az ARIMA modelltől. A következő részben megnézzük mi a hiányzó I. De még előtte beszéljünk az idősorok stationary-ról:

Idősor Stationary

Definíció szerint akkor nevezünk stationary-nak egy idősort, ha a megfigyelések értékei nem függnek a megfigyelés időpontjától. Személy szerint én ezt a meghatározást kicsit viccesnek érzem, mert lényegében azt jelenti, hogy innentől fogva nem is kell idősorról beszélnünk. A megfigyeléseket kezelhetjük egymástól független mintáknak az idő paraméter elhagyásával.

Miért szeretnénk stationary idősorokat? Pont azért mert nem kell különleges, időfüggő adatként kezelni őket. Hogy néz ki egy stationary idősor? Nem tudunk semmiféle mintát felfedezni, kb. mintha fehér zaj lenne. Értelemszerűen egy idősor, aminek trendje vagy szezonalitása van, nem fog megfelelni ennek az elvárásnak. Nézzük például a következő 9 ábrát:

Mi lesz ezek közül stationary? Az (a), (c), (e), (f) és (i) egyértelműen nem, mert trendjük van. A (d), (h) és (i)-nek szezonanitása van. A (b) stationarynak tűnik néhány kiugró értékkel. De mi van a (g)-vel? Elsőre annak is mintha lenne szezonanitása, de ha jobban megvizsgáljuk, akkor ki fog derülni, hogy nincs neki. Mégpedig azért, mert elsőre talán nem tűnik fel, de a hullámok hossza nem megjósolható. A (g) egy hiúz populáció alakulását mutatja. A hullámvölgyek a túlszaporodást követő visszaesés időszakai. Ahhoz, hogy a túlszaporodás bekövetkezzen több tényezőnek együtt kell állnia, és ezért nem megjósolható az egyes ciklusok hossza.

Az ARIMA szempontjából számunkra az (a) és (b) ábra igazán érdekes. Mint látható ez az (a) adatsor, Google részvényárfolyam, deriváltja. Az, hogy az egyik stationary, a másik nem, rá is mutat számunkra egy eljárásra, ami segítségével időfüggetlené tehetünk egy idősort.

I

Ha a fenti ARMA modellt nem az eredeti idősoron, hanem annak, akár többszörösen, derivált változatán készítjük el, akkor beszélünk ARIMA(p,d,q) modellről. Nem meglepő módon a deriválást azért szeretnék, mert így stationary lesz az idősorunk. Az, hogy hányszor deriváljuk az idősorunk fejezi ki a d paraméter. Ennek van egy “lag operatorral” leírt változata:

$(1-\alpha_1 B - \cdots - \alpha_p B^p) (1-B)^d Y_{t} = c + (1 + \beta_1 B + \cdots + \beta_q B^q) \epsilon_t$

SARIMA

Eddig olyan adatokon dolgoztunk, amikben nincs szezonalitás. Ha ezt is hozzáadjuk a modellhez akkor beszélünk $SARIMA(p,d,q)(P,D,Q)_m$ modellről. Ahol a $(P,D,Q)$ a szezonális komponense a ARIMA modellben megismert módon, az m pedig, hogy mennyi megfigyelés van egy ciklusban.

SARMAX

A fenti SARIMA szép és jó, de gyakorlati alkalmazásánál szinte egyből felmerül a kérdés: hogyan integráljunk más tulajdonságokat is a modellbe. Az ARIMA nagy problémája, hogy csak egyetlen egy tulajdonságot, idősort tud modellezni. De a világunk általában ennél sokkal összetettebb, pl: a jövőbeli GDP értékek nem csak a múltbeli GDP értékektől függnek, hanem sok másik tényezőtől is. Ennek a problémának a felismerése vezetett el az ARMAX-hoz. A misztikus X pedig a “exogenous inputs”-t rövidíti, amik a többi tulajdonság.

SARIMAX Pythonban

Az elméleti áttekintés után nézzük meg, hogy is lehet ARIMA modelletek építeni Python-ban. Gondolom nem árulok el nagy titkot, ha elmondom, hogy a legnagyobb probléma a helyes paraméterek beállítása. Egy SARIMA modell esetén összesen 7 paraméterünk van, ami kézi optimalizációja fájdalmas tud lenni. Ezek közül az m, a ciklus hossza talán a legkönnyebb. Ha biztosak vagyunk abban, hogy van szezonalitás, akkor végezzünk egy autokorrelációt és kész. Tegyük is ezt:

from scipy.signal import argrelextrema
import statsmodels.api as sm
from statsmodels.graphics.tsaplots import plot_acf
def autocorr(x):
    # autokorreláció 95% konfidencia intervallummal 
    acorr,confidence_interval = sm.tsa.acf(x, alpha=0.05, fft=True, nlags=x.shape[0])
    confidence_interval[:, 0] = confidence_interval[:, 0]-acorr
    confidence_interval[:, 1] = confidence_interval[:, 1]-acorr
    lag = 1
    r = np.NaN
    pottential_lag = []

    # az optimális lag ahol acorr a konfidencia intervallumon kívül van
    # és a lokális min vagy max a legtávolabb van a szomszédjaitól
    maxs = argrelextrema(acorr, np.greater)
    maxs = maxs[0]
    for maxi in maxs:
        # a lokális maximum nagyobb mint 95%  konfidencia intervallumon
        if acorr[maxi] > confidence_interval[maxi, 1]:
            # távolság a két szomszédtól
            dist = np.abs(acorr[maxi]-acorr[maxi-1])+np.abs(acorr[maxi]-acorr[maxi+1])
            pottential_lag.append( 
                (dist, int(maxi))
            )
    mins = argrelextrema(acorr, np.less)
    mins = mins[0]
    for mini in mins:
        # a lokális minimum kisebb mint 95% konfidencia intervallumon
        if acorr[mini] < confidence_interval[mini, 0]:
            dist = np.abs(acorr[maxi-1]-acorr[maxi])+np.abs(acorr[maxi]-acorr[maxi+1])
            pottential_lag.append( 
                (dist, int(mini))
            )
    
    if len(pottential_lag) > 0:
        lag = sorted(pottential_lag)[-1][1]
        r = acorr[lag]        
    if lag > 1:
        print('Van autokorreláció: r = {}, lag = {}'. format(r, lag))
    else: 
        print('Nincs autokorreláció.')
    plot_acf(x)
    plt.show()
    return r, lag
_, m = autocorr(train)

# ábrázoljuk
import matplotlib.pyplot as plt
plot_acf(train, lags=50)
plt.show()

Az autokorreláció ábráján jól látszik, hogy a példaadatoknál 12-es ciklusok vannak. És a számított m érték is ezt támasztja alá.

A többi paraméter optimalizációja sajnos ennél több munkát igényel, ezért egy kicsit csalni fogunk és automatizálni fogjuk. A legkézenfekvőbb megoldás az automatizálásra egy Grid search implementálása. De sebesség tekintetében ettől ne várjunk csodát. Még ha az m értékét ismerjük is, akkor is 6 paramétert kell végigpróbálni. Most nézzünk egy példa megvalósítást:

import itertools
import warnings
warnings.filterwarnings("ignore")

def sarimax_grid_search(y,seasonal_period):
    p = d = q = range(0, 3)
    pdq = list(itertools.product(p, d, q))
    seasonal_pdq = [(x[0], x[1], x[2],seasonal_period) for x in list(itertools.product(p, d, q))]
    
    mini = float('+inf')
    
    
    for param in pdq:
        for param_seasonal in seasonal_pdq:
            try:
                mod = sm.tsa.statespace.SARIMAX(y,
                                                order=param,
                                                seasonal_order=param_seasonal,
                                                enforce_stationarity=False,
                                                enforce_invertibility=False)

                results = mod.fit()
                
                if results.aic < mini:
                    mini = results.aic
                    param_mini = param
                    param_seasonal_mini = param_seasonal
                    
                print('SARIMA{}x{} - AIC:{}'.format(param, param_seasonal, results.aic))
            except Exception as e:
                print(e)
                continue
    print(f"Legjobb modell: SARIMA{param_mini}x{param_seasonal_mini}, Minimum AIC: {mini}")
    
sarimax_grid_search(train, m)

Hosszú-hosszú várakozás után a legjobb modell a SARIMA(1,1,2)x(0,2,2,12). Most végezzünk egy kis előrejelzést ezzel a modellel:

# model
model = sm.tsa.statespace.SARIMAX(train,
                          order=(1, 1, 2),
                          seasonal_order=(0, 2, 2, 12),
                          enforce_stationarity=False,
                          enforce_invertibility=False)
model_fit = model.fit()
# előrejelzés, az utolsó megfigyelés nem kell, mert annak valós értékét még nem tudjuk
yhat = model_fit.predict(start=len(train), end=len(y)-1)
# ábra
x = np.arange(y.shape[0])
plt.plot(x[:150], train, c='blue', label="tanuló adatok")
plt.plot(x[150:], yhat, c='green', label="model")
plt.plot(x[150:], test, c="red", label="valóság")
plt.legend()
plt.xlabel("Idő")
plt.ylabel("Megfigyelés")
plt.title("Auto-SARIMAX példa")
plt.show()

Az eredmény ránézésre elég jónak tűnik, de gondolom feltűnt, hogy ez SARIMA modell és nem SARIMAX. A helyzet az, hogy a statsmodel SARIMAX függvényt alkalmazzuk a fenti eljárás során, ami nem meglepő módon alkalmas SARIMAX modellek építésére. Sajnos én hirtelen nem találtam megfelelő példa adathalmazt, így csak SARIMA modellt építettem. De ha valakit érdekel: SARIMAX modell esetén annyi a dolgunk, hogy megadjuk őket az exog paraméter értékében. Majd a fit függvénynél a start_modell értékét kiterjesztjük k-1 értékkel. Ahol k az exogenous tulajdonságunk száma. Pl: ha 2 exogenous tulajdonságot használunk akkor a start_params valami ilyesmi lesz: [0, 0, 0, 0, 1]. Összességében valami ilyesmi egy SARIMAX modell:

mod = sm.tsa.statespace.SARIMAX(endog=y,
                                exorg=valami,
                                order=param,
                                seasonal_order=param_seasonal,
                                enforce_stationarity=False,
                                enforce_invertibility=False
      )
mod.fit(start_params=[0, 0, 0, 0, 1])

Auto-SARIMA

Abban az esetben, ha a SARIMA modell elegendő számunkra van egy a fenti példánál is egyszerűbb lehetőségünk: egy előre megírt könyvtár használata. R-ben az Auto-ARIMA csomag elég nagy népszerűségnek örvend, olyan nagynak, hogy Taylor G Smith eldöntötte, hogy átírja Python-ra. Ez a pmdarima projekt. Személyes véleményem szerint ez egy elég gyors és könnyen használható megvalósítás, így mindenkinek ajánlom. Az egyetlen problémája, hogy az m értéket ismerni kell. De szerencsére ezt már tudjuk.

Most készítsük el a Auto-SARIMA modellt:

import pmdarima

# Model illesztése
#   m -- ezt tudni kell 
modell = pmdarima.auto_arima(train, error_action="ignore", suppress_warnings=True, trace=True, seasonal=True, m=12)
print(modell.summary())

# Előrejelzés
elorejezes, konv_int  = modell.predict(test.shape[0], return_conf_int= True, alpha=0.01)  

# ábra
x = np.arange(y.shape[0])
plt.plot(x[:150], train, c='blue', label="tanuló adatok")
plt.plot(x[150:], elorejezes, c='green', label="legvalószínűbb modell")
plt.fill_between(
    x[150:],
    konv_int[:,0],
    konv_int[:,1],
    alpha=0.5,
    color="green",
    label="99% konfidenciaintervallum"
                )
plt.plot(x[150:], test, c="black", label="valóság")
plt.legend()
plt.xlabel("Idő")
plt.ylabel("Megfigyelés")
plt.title("Auto-SARIMA példa")
plt.show()

Ami eredménye:

Ez elég jónak tűnik. És nagyságrendekkel gyorsabb mint az előző Auto-SARIMAX módszer.

Irodalom

Dr. Aric LaBarr: The Bayesians are Coming to Time Series
Sushmitha Pulagam: Time Series forecasting using Auto ARIMA in python
Jason Brownlee: How to Create an ARIMA Model for Time Series Forecasting in Python
Selva Prabhakaran: ARIMA Model – Complete Guide to Time Series Forecasting in Python
Rob J Hyndman és George Athanasopoulos: Forecasting: Principles & Practice
Shay Palachy: Stationarity in time series analysis
Jose Marcial Portilla: Using Python and Auto ARIMA to Forecast Seasonal Time Series
Anthony Skolozdrzyk-Ardouin: A simple grid search to find (p,d,q)*(p,d,q,s) in SARIMAX forecasting algorithm