[데이터분석] 시계열 분석 1 - ARIMA

기업 내 데이터분석이 많아지면서 자연스럽게 시계열분석에 대한 Needs가 많아졌다. 마케팅/구매/SCM 등 가치 사슬 내 거의 모든 부분이 시계열 데이터에 의존하고 있는데, 기존의 통계적 접근에서 빅데이터 분석, 딥러닝까지 발전(?)해 오고 있다. Kaggle의 예제 데이터 기반으로 그 방법들을 살펴 볼 예정이다. 

 

데이터는 아래 Kaggle 에서 가져올 수 있다. 

Store Item Demand Forecasting Challenge

구조는 간단해서, 2013년부터 17년말까지 가게(store)별로 각 아이템들이 얼마나 팔렸는지를 나타내고 있다. 월별로 살펴보면 판매는 여름에 높아지고 겨울로 갈 수록 줄어드는 경향을 보인다. 

 

 

데이터 EDA는 위와 같이 간단히만 해 두고, 판매 예측을 위해 우선 전통적인 ARIMA(Auto Regressive Integrated Moving Average)를 적용해 본다. ARIMA는 AR과 MA가 합쳐진 것으로, AR은 p시점 이전의 데이터가, MA는 최근의 추세(평균) 변화가 현재의 결과에 영향을 준다는 것이다. 코드를 통해 단계별 결과를 확인해 보자.

 

1) 전처리

 

아래는 전처리 하는 부분으로 우선 월별로 판매를 합산하고 훈련과 테스트 데이터로 나눈다.

 

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from datetime import datetime

from pylab import rcParams
import statsmodels.api as sm
import warnings
import itertools
from statsmodels.tsa.ar_model import AR
from statsmodels.tsa.arima_model import ARIMA

# Kaggle에서 데이터를 받아오면 train.csv와 test.csv가 있는데 train.csv를 
df = pd.read_csv("data/train.csv")

# string 형태의 date 필드를 datetime 형태로 변환
df['date'] =  pd.to_datetime(df['date'])

# date 필드를 index로 설정
df = df.set_index('date')

# 빠르게 테스트 해 보기 위해 월별로 아이템 판매 예측을 해 보기로 함
salesbymonth = df.sales.resample('M').sum()

#2013-2016 데이터를 train으로 2017 데이터를 test로 분리
split = "2017-01-01"
salesbymonth_train= salesbymonth[:split]
salesbymonth_test= salesbymonth[split:]
salesbymonth_test_final=salesbymonth_test.copy()

# 데이터를 시즌별로 분해해서 살펴 봄
rcParams['figure.figsize'] = 18, 8
decomposition = sm.tsa.seasonal_decompose(salesbymonth_train, model='additive')
fig = decomposition.plot()
plt.show()

 

statsmodels 패키지의 seasonal_decompose함수는 seasonal 데이터를 트렌드와 주기성을 살펴보는데 사용하는 데 결과는 아래와 같다.

 

 

점진적으로 증가하는 트렌를 보이고 있으며 계절별로 주기성이 뚜렷한 편이다. 네번째 항목은 잔차를 나타낸다.

 

2) 파라미터 찾기

 

이제 본격적으로 ARIMA 모델을 만들어야 한다. ARIMA 모델을 만들기 위해 p, d, q 값을 지정해 주어야 하는데 일반적으로 p가 0이면 MA 모델을 따르고, q가 0이면 AR 모델을 따른다고 한다. d가 0이면 정상성을 보유한 모델이라고 보는데, 발생확률이 변하지 않는 ARMA 모델이라고 간주하면 된다고 한다. 우린 아직 이 데이터가 어떤 모형을 따르는 지 모르기 때문에 최적의 결과를 가져오는 p,d,q 값을 찾아야한다. 그 방법은 크게 두가지로 임의의 값들을 다 넣어 테스트 해보거나, pmdarima 패키지의 auto_arima 함수로 찾아내면 된다. 

 

#방법 1. p,d,q의 조합을 만들어 하나하나 ARIMA 모델을 돌려봄
p = d = q = range(0, 2)

import itertools
pdqa = list(itertools.product(p, d, q))

seasonal_pdq = [(x[0], x[1], x[2], 12) for x in list(itertools.product(p, d, q))]

for param in pdqa:
    for param_seasonal in seasonal_pdq:
        try:
            mod = sm.tsa.statespace.SARIMAX(salesbymonth_train, order=param, seasonal_order=param_seasonal,enforce_stationarity=False,enforce_invertibility=False)                                
            results = mod.fit()
            print('ARIMA{}x{}12 - AIC:{}'.format(param, param_seasonal, results.aic))
        except:
            continue
            
#방법 2. auto_arima 함수로 자동 추출
from pmdarima import auto_arima
stepwise_model = auto_arima(salesbymonth_train, start_p=1, start_q=1,
                           max_p=3, max_q=3, m=12,
                           start_P=0, seasonal=True,
                           d=1, D=1, trace=True,
                           error_action='ignore',  
                           suppress_warnings=True, 
                           stepwise=True)

 

후자의 결과는 아래와 같이 나타나는데 전자의 결과에서도 AIC(Akaike Information Criterion) 값이 제일 작은 조합을 선택하면 된다.

 

Performing stepwise search to minimize aic
 ARIMA(1,1,1)(0,1,1)[12]             : AIC=775.474, Time=0.07 sec
 ARIMA(0,1,0)(0,1,0)[12]             : AIC=784.706, Time=0.01 sec
 ARIMA(1,1,0)(1,1,0)[12]             : AIC=774.061, Time=0.04 sec
 ARIMA(0,1,1)(0,1,1)[12]             : AIC=773.779, Time=0.05 sec
 ARIMA(0,1,1)(0,1,0)[12]             : AIC=786.609, Time=0.01 sec
 ARIMA(0,1,1)(1,1,1)[12]             : AIC=775.686, Time=0.20 sec
 ARIMA(0,1,1)(0,1,2)[12]             : AIC=775.734, Time=0.11 sec
 ARIMA(0,1,1)(1,1,0)[12]             : AIC=774.156, Time=0.04 sec
 ARIMA(0,1,1)(1,1,2)[12]             : AIC=inf, Time=0.55 sec
 ARIMA(0,1,0)(0,1,1)[12]             : AIC=770.820, Time=0.03 sec
 ARIMA(0,1,0)(1,1,1)[12]             : AIC=772.724, Time=0.07 sec
 ARIMA(0,1,0)(0,1,2)[12]             : AIC=772.374, Time=0.11 sec
 ARIMA(0,1,0)(1,1,0)[12]             : AIC=767.799, Time=0.05 sec
 ARIMA(0,1,0)(2,1,0)[12]             : AIC=772.845, Time=0.09 sec
 ARIMA(0,1,0)(2,1,1)[12]             : AIC=inf, Time=0.41 sec
 ARIMA(1,1,1)(1,1,0)[12]             : AIC=775.833, Time=0.23 sec
 ARIMA(0,1,0)(1,1,0)[12] intercept   : AIC=773.761, Time=0.09 sec

Best model:  ARIMA(0,1,0)(1,1,0)[12]          
Total fit time: 2.150 seconds

 

3) 학습 및 테스트

 

이제 ARIMA 모델을 학습해 본다. 학습 코드는 간단하다.

 

SARIMAMonth = sm.tsa.statespace.SARIMAX(salesbymonth, order=(0, 1, 0), seasonal_order=(1, 1, 0, 12) ,enforce_stationarity=False,enforce_invertibility=False)

SARIMA_results_month = SARIMAMonth.fit()

SARIMA_results_month.plot_diagnostics(figsize=(16, 8))
plt.show()

 

학습 결과는 아래와 같다. 상관도가 높게 잘 나타났다.

 

 

만들어 두었던 예측 데이터를 넣어 예측하고 결과를 비교해 본다.

 

# 2017년 12개월 데이터로 예측
SARIMA_predict_month_1 = SARIMA_results_month.predict(start=48,end=60)

# 결과 비교를 위해 기존에 마련해둔 test데이터에 결과를 붙임
salesbymonth_test_final['SeasonalARIMA'] = SARIMA_predict_month_1

# RMSE를 살펴 봄
RMSE_Month_Seasonal_ARIMA  = np.mean(np.sqrt((salesbymonth_test_final['SeasonalARIMA'] - salesbymonth_test_final['sales']) ** 2)) 
print(RMSE_Month_Seasonal_ARIMA)
-> 12190.886296802802

# test 데이터와 예측 결과치를 비교
salesbymonth_test_final[1:].plot()

 

비교적 결과를 제대로 예측하고 있음을 확인할 수 있다.