파이썬에서 SARIMA 실행하기 - SARIMAX

파이썬의 statsmodels 패키지는 ARIMA, SARIMA 시계열 선형 모델을 지원합니다. 이번 포스트에서는 statsmodels의 SARIMAX 함수를 이용하여 시계열 데이터의 패턴을 학습하고 예측하는 걸 해보도록 하겠습니다. 먼저, SARIMAX 함수를 import 합니다. SARIMAX의 X는 외부 변수를 나타내는 eXogeneous의 줄임말로 자기 자신 (endogeneous) 뿐만 아니라 외부 변수까지 학습과 예측에 포함할 수 있다는 것입니다.

from statsmodels.tsa.statespace.sarimax import SARIMAX

 SARIMAX의 주요 파라미터는 다음과 같습니다.

파라미터	Description
endog	관측된 시계열 데이터 (endogeneous 데이터를 말합니다)
exog	관측된 시계열에 영향을 미치는 외부 변수 데이터 (exogeneous  데이터를 말합니다)
order	ARIMA의 $(p,d,q)$
seasonal_order	SARIMA의 seasonal component인 $(P,D,Q)_S$
enforce_stationary	AR 항이 stationary를 만족하게끔 강제하는 것으로 디폴트는 True
enforce_invertibility	MA 항이 stationar를 만족하게끔 강제하는 것으로 디폴트는 False

 

Train

SARIMAX 모델을 선언하기 위해서는 위의 파라미터를 설정해주어야 하며 특히 $(p,d,q)\times (P,D,Q)_S$ 파라미터를 설정해주어야 합니다. 이를 위해서는 각 파라미터가 가질 수 있는 값을 설정하고 모든 파라마티의 조합마다 SARIMAX를 훈련시켜 가장 낮은 AIC/BIC 값을 가지는 파라미터 조합을 선택합니다. 이를 grid search라고 합니다. 

Grid search는 파라미터의 개수가 많아지고 파라미터가 가질 수 있는 선택지가 많아질수록 탐색 공간이 지나치게 넓어져 탐색에 오랜 시간이 걸립니다. 여기서는 $S$는 12로 미리 안다고 가정하고 (ACF 등으로 미리 알 수 있습니다) $p,d,q,P,D,Q$는 0 에서 3 사이의 값을 가진다고 가정한 후 grid search를 진행하겠습니다. 먼저, itertools의 product 함수를 이용하여 가능한 모든 파라미터의 조합을 생성하고, 

import itertools

## Define the p, d and q parameters to take any value between 0 and 2
p = d = q = range(0, 3)

## Generate all different combinations of p, d and q triplets
pdq = list(itertools.product(p, d, q))

# generate all different combinations of seasonal p, q and q triplets
seasonal_pdq = [(x[0], x[1], x[2], 12) for x in list(itertools.product(p, d, q))]

AIC가 가장 낮은 파라미터 조합을 선택합니다.

best_aic = np.inf
best_pdq = None
best_seasonal_pdq = None
tmp_model = None
best_mdl = None

for param in pdq:
    for param_seasonal in seasonal_pdq:
        try:
            tmp_mdl = sm.tsa.statespace.SARIMAX(y_train, exog=None, order = param,
                                                seasonal_order = param_seasonal,
                                                enforce_stationarity=True,
                                                enforce_invertibility=True)
            
            res = tmp_mdl.fit()        
            print("SARIMAX{}x{}12 - AIC:{}".format(param, param_seasonal, res.aic))
            if res.aic < best_aic:
                best_aic = res.aic
                best_pdq = param
                best_seasonal_pdq = param_seasonal
                best_mdl = tmp_mdl
        except:
            #print("Unexpected error:", sys.exc_info()[0])
            continue
print("Best SARIMAX{}x{}12 model - AIC:{}".format(best_pdq, best_seasonal_pdq, best_aic))

Grid search를 통해 최상의 파라미터 조합을 선택하고 데이터에 대해 fitting 합니다.

## Define SARIMAX model and fit it to the data
mdl = sm.tsa.statespace.SARIMAX(endog=y_train, order=(1, 1, 0),
                                seasonal_order=(1, 2, 1, 12),
                                enforce_stationarity=True,
                                enforce_invertibility=True)
res = mdl.fit()

모델의 fit 결과에서 summary를 호출하면 모델이 데이터에 잘 fitting 되었는지 알 수 있습니다.

res.summary()

추정한 각 계수와 p-value, 지난 포스트에서 살펴봤던 모델의 검정 방법인 Ljung-Box, Jarque-Bera 테스트 결과를 볼 수 있습니다.

Ljung-Box 테스트는 fitting 이후의 잔차가 래그에 따라 correlated 되어 있는지 판단하는 것으로 잔차가 white noise를 따르는지 판단하는 테스트입니다. 따라서 Ljung-Box 테스트의 귀무 가설은 잔차가 white noise를 따른다 이고 p-value가 낮으면 귀무 가설을 기각하여 white noise를 따르지 않아 각 시간대 별 잔차가 correlated 되어있다고 판단할 수 있습니다.

Jarque-Bera 테스트는 잔차의 분포가 normal distribution을 따르는지에 대한 테스트로 normal distribution은 skewness가 0, kurtosis가 3인 것을 이용합니다. 이 테스트 또한 귀무 가설은 잔차가 normal distribution을 따른다 이고 p-value가 낮으면 귀무 가설을 기각합니다.

Heteroskedasticity는 각 시간대 별 잔차의 분산이 일정한지 (homoscedastic) 에 대한 테스트입니다. 이 테스트의 귀무 가설은 잔차의 분산이 일정하다 (homoscedastic)이고 p-value가 낮으면 귀무 가설을 기각합니다.

이를 종합해서 볼 때, fitting 모델의 결과의 p-value는 모든 테스트에 대해 0.05 이상을 가지므로 fitting이 잘 되었을 때의 잔차의 검정 테스트를 모두 잘 통과했음을 알 수 있습니다. plot_diagnostic을 이용하면 그림으로 확인할 수 있습니다.

res.plot_diagnostics(figsize=(16, 10))

 

Prediction

SARIMAX를 통해 예측한 이후에 get_prediction 혹은 get_forecast 메소드로 이후 시점에 대해 예측할 수 있습니다.

## Fit model to data
res = sm.tsa.statespace.SARIMAX(y_train,
                                order=(1, 1, 0),
                                seasonal_order=(1, 2, 1, 12),
                                enforce_stationarity=True,
                                enforce_invertibility=True).fit()

## In-sample-prediction and confidence bounds
pred = res.get_prediction(start=pd.to_datetime('1958-12-01'), 
                          end=pd.to_datetime('1960-12-01'),
                          dynamic=True)

get_prediction 함수는 예측할 기간의 start/end를 선택하고 dynamic 파라미터를 선택합니다. dynamic 파라미터가 True일 경우 in-sample prediction (학습한 데이터에 대한 예측)을 수행할 때, 예측값 대신 실제값을 이용하겠다는 것으로 여기서는 out-of-sample forecast 이어서 사용 가능한 실제값이 없으므로 True/False 차이가 없습니다. 또한, exogeneous 데이터를 훈련에 포함할 경우 exog 파라미터에 예측에 사용할 외부 수열을 넣을 수 있습니다.

get_prediction의 결과물로부터 prediction mean과 prediction의 95% 신뢰구간을 얻을 수 있습니다. 이를 plot하면 다음과 같습니다.

pred_ci = pred.conf_int()

## Plot in-sample-prediction
ax = y['1949':].plot(label='Observed',color='#006699');
pred.predicted_mean.plot(ax=ax, label='One-step Ahead Prediction', alpha=.7, color='#ff0066');

## Draw confidence bound (gray)
ax.fill_between(pred_ci.index, 
                pred_ci.iloc[:, 0], 
                pred_ci.iloc[:, 1], color='#ff0066', alpha=.25);

## style the plot
ax.fill_betweenx(ax.get_ylim(), pd.to_datetime('1958-12-01'), y.index[-1], alpha=.15, zorder=-1, color='grey');
ax.set_xlabel('Date');
ax.set_ylabel('Passengers');
plt.legend(loc='upper left');
plt.show()

또한, get_forecast(steps=n) 메소드로 out-of-sample forecast를 수행할 수 있습니다. get_prediction과 마찬가지로 예측값과 예측값의 신뢰 구간을 얻을 수 있습니다. steps 파라미터를 지정하면 훈련에 사용한 데이터 바로 다음부터 지정한 step만큼 예측합니다.

## Build model and fit
res = sm.tsa.statespace.SARIMAX(y,
                                order=(2, 1, 3),
                                seasonal_order=(1, 2, 1, 12),
                                enforce_stationarity=True,
                                enforce_invertibility=True).fit()

## Get forecast 120 steps ahead in future
pred_uc = res.get_forecast(steps=120)

## Get confidence intervals of forecasts
pred_ci = pred_uc.conf_int()

## Plot time series and long-term forecast
ax = y.plot(label='Observed', figsize=(16, 8), color='#006699');
pred_uc.predicted_mean.plot(ax=ax, label='Forecast', color='#ff0066');
ax.fill_between(pred_ci.index,
                pred_ci.iloc[:, 0],
                pred_ci.iloc[:, 1], color='#ff0066', alpha=.25);
ax.set_xlabel('Date')
ax.set_ylabel('Passengers')
plt.legend(loc='upper left')
plt.show()