시계열 모델링 - ARIMA

"금융 파이썬 쿡북 Ch3. 시계열 모델링 "의 내용을 기반으로 작성하였습니다. (실습 깃헙)

---

 

1. ARIMA 상세

ARIMA 모델의 내용은 본 책과 고려대 김성범 교수님의 유튜브 강의를 참고하였다. (좋은 내용 공유 감사드립니다.)

https://www.youtube.com/watch?v=ma_L2YRWMHI 

 

 ARIMA는 데이터의 자기 상관을 설명하는 방법을 사용한다. 이는 ARMA(Autoregressive Moving Average) 모델을 확장한 것이다.

 

ARIMA의 구성요소에 대해 살펴보자.

 

• AR (Autoregressive) 모델
  • 관측값과 지연값 사이의 관계를 사용
  • 모멘텀과 평균 회귀 효과를 반영

DMQA 김성범 교수님 ARIMA 강의 슬라이드

 

• I (Integration)
  • 시계열의 차분을 의미, 이전 기간의 값을 현재 기간의 값에서 차감해 시계열을 정상성으로 만들 수 있음
    
    
• MA (Moving Average) 모델
  • 백색 노이즈 분포를 따르는 과거 시점의 error들을 설명변수로 하는 모델
  • 관측된 시계열에서 영향을 미치는 예측할 수 없는 충격(잔차에서의 관측)을 설명

 

 

 위 세가지 구성요소 중 AR과 MA가  합쳐져 ARMA를 구성하게 된다.

 

 이 때 I, 즉 차분을 통해 정상성이 확보된 이후에 시계열 모델링을 수행하는 것, 다시 말해 AR, I, MA 가 합쳐저 ARIMA가 되는 것이다.

이 세가지 요소는 한번에 적합화되어 특정 시계열 데이터 특성에 맞는 모델링을 수행한다.

 

 ARIMA는 모델링의 방식에 따라 ARIMA(p,d,q)로 표기되는데 각각의 의미는 위 수식들에서 확인할 수 있다.

 

• p : AR 모델을 구성하는 과거 시점 independent variable 갯수
• d : 차분의 차수 (1차, 2차 등등)
• q : MA 모델을 구성하는 과거 시점 independent variable 갯수

 즉, d차 차분한 데이터 위에 AR(p), MA(q)를 합친 모형인 것이다.

 

 

p,d,q의 매개변수 조합에 따라 일부 특수한 다음의 경우를 얻을 수 있다.

 

• ARIMA(0,0,0) : 백색 노이즈
• ARIMA(0,1,0) ,상수 없음 : 랜덤 워크
• ARIMA(p,0,q) : ARMA(p, q)
• ARIMA(p,0,0) : AR(p) 모델
• ARIMA(0,0,q) : MA(q) 모델
• ARIMA(0,1,2) : 댐핑된 홀트 모델
• ARIMA(0,1,1) , 상수없음 : SES 모델
• ARIMA(0,2,2) : 가산적(additive) 오차를 가진 홀트 선형 모델

 

2. ARIMA 기본 구현 과정

 

 구글의 주가 데이터로 구현해보자.

 

import yfinance as yf
import pandas as pd
import numpy as np
from statsmodels.tsa.arima_model import ARIMA
import statsmodels.api as sm
from statsmodels.graphics.tsaplots import plot_acf
from statsmodels.stats.diagnostic import acorr_ljungbox
import scipy.stats as scs

df = yf.download('GOOG',
                 start='2015-01-01',
                 end='2018-12-31',
                 progress=False)

goog = df.resample('W') \
         .last() \
         .rename(columns={'Adj Close': 'adj_close'}) \
         .adj_close

 

 우선 시계열 데이터 자체의 정상성부터 확인해보자. 아래 정상성 검정 내용은 이전 포스팅에서 확인할 수 있다.

 

# 정상성 확인
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from statsmodels.tsa.stattools import adfuller, kpss

# ADF (Augmented Dickey Fuller) 검정

def adf_test(x) :
  indices = ['Test Statistics', 'p-value', '# of Lags Used', '# of Observations Used']

  # AIC = 고려된 지연 수(lags)를 AIC(Akaike Information Criterion) 기반 자동 선택
  adf_test = adfuller(x, autolag = 'AIC')
  results = pd.Series(adf_test[0:4], index = indices)

  for key, value in adf_test[4].items() :
    results[f'Critical Value ({key})'] = value
  
  return results

# KPSS (Kwiatkowski-Phillips-Schmidt-Shin) 검정

# c = 시계열이 level 정상성이라는 귀무가설
# ct = 시게열이 trend 정상성이라는 귀무가설 (계열에서 trend를 제거하면 level 정상성이 됨)
def kpss_test(x, h0_type = 'c') :
  indices = ['Test Statistics', 'p-value', '# of Lags Used']
  kpss_test = kpss(x, regression = h0_type)
  results = pd.Series(kpss_test[0:3], index = indices)

  for key, value in kpss_test[3].items() :
    results[f'Critical Value ({key})'] = value
  
  return results


# 통합
def test_autocorrelation(x, n_lags=40, alpha=0.05, h0_type='c'):

    adf_results = adf_test(x)
    kpss_results = kpss_test(x, h0_type=h0_type)

    print('ADF test statistic: {:.2f} (p-val: {:.2f})'.format(adf_results['Test Statistics'],
                                                             adf_results['p-value']))
    print('KPSS test statistic: {:.2f} (p-val: {:.2f})'.format(kpss_results['Test Statistics'],
                                                              kpss_results['p-value']))

    fig, ax = plt.subplots(2, figsize=(16, 8))
    plot_acf(x, ax=ax[0], lags=n_lags, alpha=alpha)
    plot_pacf(x, ax=ax[1], lags=n_lags, alpha=alpha)

    return fig

test_autocorrelation(df['Adj Close'])

 

 

 통계적 수치와 ACF / PACF를 봤을 때 비정상성을 띔을 확인할 수 있다.

 

이 때 ARIMA는 정상성을 전제로 하기에, 차분을 통해 정상성을 달성할 수 있을지 확인해보자. (비정상성 교정은 이전 포스팅에서 확인할 수 있다.)

 

goog_diff = goog.diff().dropna()

fig, ax = plt.subplots(2, sharex=True)
goog.plot(title = "Google's stock price", ax=ax[0])
goog_diff.plot(ax=ax[1], title='First Differences')

plt.tight_layout()
plt.show()

 원래의 시계열 데이터는 비정상성을 띈다고 검정했던 것과 같이 일정 추세가 보이는 반면, 차분된 시계열은 정상성을 띄는 것 처럼 보인다.

 

이 또한 위 정상성 검정 함수를 통해 정확히 검정해보자.

 

fig = test_autocorrelation(goog_diff)

plt.tight_layout()
plt.show()

 

 차분된 결과는 정상성을 가짐을 알 수 있고, 이렇게 1차 차분한 결과를 바탕으로 ARIMA 모델을 구현한다.

 

앞서 ARIMA 모델에서는 ARIMA(p,d,q)라는 매개변수를 설정해야 하는데 이는 어떻게 설정할 수 있을까?

 

• d = 정상성을 띄게 된 차분의 차수 : 여기서는 1차 차분 이후 정상성으로 검정되었으므로 d = 1
• p = PACF에서 신뢰구간을 교차한 후 마지막 lag : 여기서는 p = 2 (PACF lag는 3부터 신뢰구간에 포함)
• q = ACF에서 신뢰구간을 교차한 후 마지막 lag : 여기서는 q = 1 (ACF lag는 2부터 신뢰구간에 포함)

 위 rule of thumb 대로라면 ARIMA(2,1,1) 모델을 사용해야 한다.

 

 위 방식을 달리 설명한 다음의 도표대로 모델을 설정할 수 있다.

 

# ARIMA 모델 설정

arima = ARIMA(goog, order = (2,1,1)).fit(disp = 0)
arima.summary()

 

 이렇게 적합화 모델이 잘 적합화 되었는지는 잔차를 기반으로 판단한다.

 

# Diagnosing the fit of the model based on its residual

def arima_diagnostics(resids, n_lags=40):
     
    # create placeholder subplots
    fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2)

    r = resids
    resids = (r - np.nanmean(r)) / np.nanstd(r)
    resids_nonmissing = resids[~(np.isnan(resids))]
    
    # residuals over time
    sns.lineplot(x=np.arange(len(resids)), y=resids, ax=ax1)
    ax1.set_title('Standardized residuals')

    # distribution of residuals
    x_lim = (-1.96 * 2, 1.96 * 2)
    r_range = np.linspace(x_lim[0], x_lim[1])
    norm_pdf = scs.norm.pdf(r_range)
    
    sns.distplot(resids_nonmissing, hist=True, kde=True, 
                 norm_hist=True, ax=ax2)
    ax2.plot(r_range, norm_pdf, 'g', lw=2, label='N(0,1)')
    ax2.set_title('Distribution of standardized residuals')
    ax2.set_xlim(x_lim)
    ax2.legend()
        
    # Q-Q plot
    qq = sm.qqplot(resids_nonmissing, line='s', ax=ax3)
    ax3.set_title('Q-Q plot')

    # ACF plot
    plot_acf(resids, ax=ax4, lags=n_lags, alpha=0.05)
    ax4.set_title('ACF plot')

    return fig

 

 

 잔차는 정상성을 띄며, 평균도 0에 가깝고 ACF 함수에서도 이상은 없어 보인다.

 

하지만 QQ plot 상에서 분포의 꼬리가 정규보다 무거워 보이기에 추가적으로 LJ 상자 (Ljung-Box)를 통해 잔차에 상관관계가 없는지 확인해 본다.

 

# Ljung-Box's test : 잔차에 상관관계가 없는지 

ljung_box_results = acorr_ljungbox(arima.resid)

fig, ax = plt.subplots(1, figsize=[16, 5])
sns.scatterplot(x=range(len(ljung_box_results[1])), 
                y=ljung_box_results[1], 
                ax=ax)
ax.axhline(0.05, ls='--', c='r')
ax.set(title="Ljung-Box test's results",
       xlabel='Lag',
       ylabel='p-value')

plt.tight_layout()
plt.show()

 

 

 이 검정에 대한 자세한 설명은 아래 링크에서 확인할 수 있다.

https://m.blog.naver.com/risk_girl/220834418182

Ljung-Box Q 통계량

 

 

3. AUTO ARIMA

 

 위 과정을 보면 모델을 설정하는데 있어서 정확하기 보다는 휴리스틱하게 결정하는 부분이 많았다.

 

이 때 시계열에 가장 적합한 모델을 자동으로 찾아주는 pmdarima의 auto_arima 기능이 있다.

 

여기서 가장 적합하다고 판단할 척도는 보통 AIC (Akaike Information Criteria)를 사용한다. 이는 모델의 적합성과 단순성 간의 균형을 제공하고 과적합과 과소적합을 핸들링한다. 자세한 내용은 아래 링크에서 확인할 수 있다.

 

https://rk1993.tistory.com/entry/AIC-BIC-Mallows-Cp-%EC%89%BD%EA%B2%8C-%EC%9D%B4%ED%95%B4%ED%95%98%EA%B8%B0

AIC, BIC, Mallow's Cp 쉽게 이해하기

 

 

# AUTO ARIMA : 지정된 범위에서 가능한 여러 매개 변수를 사용해 반복해 가며 그 중 AIC가 가장 낮은 모델을 선택

!pip install pmdarima

import pmdarima as pm

model = pm.auto_arima(goog, 
                      error_action = 'ignore', 
                      suppress_warnings = True,
                      seasonal = False)

model.summary()

 

 최적의 모델을 반환해주었다. 결과인 SARIMAX에 대해서는 이후 설명할 예정이다. ARIMA의 확장된 모델로 여기서는 최적의 모델로 반환했다 정도로만 알고 넘어가자.

 

조금 다르게 탐색을 진행하면 다음과 같다,

 

auto_arima = pm.auto_arima(goog,
                           error_action='ignore',
                           suppress_warnings=True,
                           seasonal=False,
                           stepwise=False,
                           approximation=False,
                           n_jobs=-1)
auto_arima.summary()

 

 다른 모델을 반환했음을 알 수 있다.

 

 

4. Prediction

 

 최적의 모델을 찾았으므로 실제 예측을 진행해보자.

최적의 모델을 찾기 전 임의 설정했던 ARIMA(2,1,1)과 비교해보자.

 

df = yf.download('GOOG',
                 start='2019-01-01',
                 end='2019-03-31',
                 progress=False)

print(f'Downloaded {df.shape[0]} rows of data.')

test = df.resample('W') \
         .last() \
         .rename(columns={'Adj Close': 'adj_close'}) \
         .adj_close

n_forecasts = len(test)

arima_pred = arima.forecast(n_forecasts)

arima_pred = [pd.DataFrame(arima_pred[0], columns=['prediction']),
              pd.DataFrame(arima_pred[2], columns=['ci_lower', 
                                                   'ci_upper'])]
arima_pred = pd.concat(arima_pred, axis=1).set_index(test.index)

auto_arima_pred = auto_arima.predict(n_periods=n_forecasts, 
                                     return_conf_int=True, 
                                     alpha=0.05)

aa = pd.DataFrame(auto_arima_pred[0], columns=['prediction'])

bb = pd.DataFrame(auto_arima_pred[1], columns=['ci_lower', 'ci_upper'])

aa['ci_lower'] = bb['ci_lower']
aa['ci_upper'] = bb['ci_upper']

auto_arima_pred = aa.set_index(test.index)

ig, ax = plt.subplots(1)

ax = sns.lineplot(data=test, color=COLORS[0], label='Actual')

ax.plot(arima_pred.prediction, c=COLORS[1], label='ARIMA(2,1,1)')
ax.fill_between(arima_pred.index,
                arima_pred.ci_lower,
                arima_pred.ci_upper,
                alpha=0.3, 
                facecolor=COLORS[1])

ax.plot(auto_arima_pred.prediction, c=COLORS[2], 
        label='ARIMA(3,1,2)')
ax.fill_between(auto_arima_pred.index,
                auto_arima_pred.ci_lower,
                auto_arima_pred.ci_upper,
                alpha=0.2, 
                facecolor=COLORS[2])

ax.set(title="Google's stock price  - actual vs. predicted", 
       xlabel='Date', 
       ylabel='Price ($)')
ax.legend(loc='upper left')

plt.tight_layout()

plt.show()

 

 여기서는 큰 차이는 보이지 않지만 auto_arima를 통해 찾은 최적 모델이 더 유연하게 예측함을 알 수 있다.

 

하지만 두 모델 모두 신뢰 구간 내에서 예측하고 있다.

 

---

(해당 게시물 학습을 위한 임의적 설정이므로, 수익을 대변하지 않습니다.)

 

Ref

• 금융 파이썬 쿡북, 에릭 르윈슨.
• https://www.youtube.com/watch?v=ma_L2YRWMHI 
• https://m.blog.naver.com/risk_girl/220834418182
• https://rk1993.tistory.com/entry/AIC-BIC-Mallows-Cp-%EC%89%BD%EA%B2%8C-%EC%9D%B4%ED%95%B4%ED%95%98%EA%B8%B0

Ljung-Box Q 통계량