[시계열분석] 시계열 알고리즘 - 선형확률과정의 분석사이클 자동화 (Auto ARIMA)(2) - 항공사 승객수요 Auro-ARIMA 모델링

항공사 승객수요 데이터 모델링

# 라이브러리 호출
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import statsmodels.api as sm
%reload_ext autoreload
%autoreload 2
from module import stationarity_adf_test, stationarity_kpss_test

# 데이터 준비
data = sm.datasets.get_rdataset("AirPassengers")
raw = data.data.copy()

# 데이터 전처리
## 시간 인덱싱
if 'time' in raw.columns:
    raw.index = pd.date_range(start='1/1/1949', periods=len(raw['time']), freq='M')
    del raw['time']

## 정상성 테스트
### 미변환
candidate_none = raw.copy()
display(stationarity_adf_test(candidate_none.values.flatten(), []))
display(stationarity_kpss_test(candidate_none.values.flatten(), []))
sm.graphics.tsa.plot_acf(candidate_none, lags=100, use_vlines=True)
plt.tight_layout()
plt.show()

### 로그 변환
candidate_trend = np.log(raw).copy()
display(stationarity_adf_test(candidate_trend.values.flatten(), []))
display(stationarity_kpss_test(candidate_trend.values.flatten(), []))
sm.graphics.tsa.plot_acf(candidate_trend, lags=100, use_vlines=True)
plt.tight_layout()
plt.show()

trend_diff_order_initial = 0
result = stationarity_adf_test(candidate_trend.values.flatten(), []).T
if result['p-value'].values.flatten() < 0.1:
    trend_diff_order = trend_diff_order_initial
else:
    trend_diff_order = trend_diff_order_initial + 1
print('Trend Difference: ', trend_diff_order)

### 로그+추세차분 변환
candidate_seasonal = candidate_trend.diff(trend_diff_order).dropna().copy()
display(stationarity_adf_test(candidate_seasonal.values.flatten(), []))
display(stationarity_kpss_test(candidate_seasonal.values.flatten(), []))
sm.graphics.tsa.plot_acf(candidate_seasonal, lags=100, use_vlines=True)
plt.tight_layout()
plt.show()

seasonal_diff_order = sm.tsa.acf(candidate_seasonal)[1:].argmax() + 1
print('Seasonal Difference: ', seasonal_diff_order)

### 로그+추세차분+계절차분 변환
candidate_final = candidate_seasonal.diff(seasonal_diff_order).dropna().copy()
display(stationarity_adf_test(candidate_final.values.flatten(), []))
display(stationarity_kpss_test(candidate_final.values.flatten(), []))
sm.graphics.tsa.plot_acf(candidate_final, lags=100, use_vlines=True)
plt.tight_layout()
plt.show()

## 데이터 로그 변환전

• ADF 테스트 결과 p-value 가 0.05보다 크므로 시간의존 구조 -> 비정상상태
• KPSS 테스트 결과 p-value가 0.05보다 작으므로 시간의존 구조 -> 비정상상태
• ACF 플롯도 추세와 계절성이 존재하고, 유의하지 않은 값들이 많다.

## 데이터 로그 변환 후

• ADF 테스트 결과 p-value 가 0.05보다 크므로 시간의존 구조 -> 비정상상태
• KPSS 테스트 결과 p-value가 0.05보다 작으므로 시간의존 구조 -> 비정상상태
• 하지만 p-value가 절반정도 감소하였다. -> 조금더 유의해졌음
• ACF 플롯도 추세와 계절성이 존재하고, 유의하지 않은 값들이 많다.

## 데이터 로그 변환 후 1회 추세 차분.

• ADF 테스트 결과 p-value 가 0.05보다 크므로 시간의존 구조 -> 비정상상태
• KPSS 테스트 결과 p-value가 0.05보다 크므로 시간의존 구조가 아님 -> 정상상태
• acf 플롯을 보면 여전히 특정한 차수가 반복되고 있음 -> 계절성이 제거가 안되었음
• 단순 로그변환 때보다 많이 개선되었다.

## 데이터 로그 변환 후 1회 추세 차분 후 계절 차분

• ADF 테스트 결과 p-value 가 0.05보다 작으므로 시간의존 구조가 아님 -> 비정상상태
• KPSS 테스트 결과 p-value가 0.05보다 크므로 시간의존 구조가 아님 -> 정상상태
• acf 플롯을 보면 값들이 대체로 유의하고, 추세성과 계절성 패턴이 나타나지 않는다.

 

항공사 승객수요 Auro-ARIMA 모델링

## 최종 타겟 선정 및 Train/Test 데이터 분리
candidate = candidate_trend.copy()
split_date = '1958-01-01'
Y_train = candidate[candidate.index < split_date]
Y_test = candidate[candidate.index >= split_date]

## 시각화 및 모수추론(p=1, q=0, d=1, P=1, Q=1, D(m)=12)
plt.figure(figsize=(14,4))
sm.tsa.graphics.plot_acf(Y_train, lags=50, alpha=0.05, use_vlines=True, ax=plt.subplot(121))
sm.tsa.graphics.plot_pacf(Y_train, lags=50, alpha=0.05, use_vlines=True, ax=plt.subplot(122))
plt.show()

# 모델링
## SARIMAX
logarithm, differencing = True, False

## 아래 부분을 돌려가면서 실행해본다.
# fit_ts_sarimax = sm.tsa.SARIMAX(Y_train, order=(1,trend_diff_order,0), trend='ct').fit()
# fit_ts_sarimax = sm.tsa.SARIMAX(Y_train, order=(1,trend_diff_order,0), 
#                                 seasonal_order=(0,0,1,seasonal_diff_order), trend='c').fit()
# fit_ts_sarimax = sm.tsa.SARIMAX(Y_train, order=(1,trend_diff_order,0), 
#                                 seasonal_order=(1,0,0,seasonal_diff_order), trend='c').fit()
fit_ts_sarimax = sm.tsa.SARIMAX(Y_train, order=(1,trend_diff_order,0), 
                                seasonal_order=(1,0,1,seasonal_diff_order), trend='c').fit()
display(fit_ts_sarimax.summary())
pred_tr_ts_sarimax = fit_ts_sarimax.predict() ## 점추정값
pred_te_ts_sarimax = fit_ts_sarimax.get_forecast(len(Y_test)).predicted_mean  ## 점추정 + 구간추정
pred_te_ts_sarimax_ci = fit_ts_sarimax.get_forecast(len(Y_test)).conf_int()
## 비정상성으로 변환
if logarithm: ##로그를 취한것을 원래 값으로 변환
    Y_train = np.exp(Y_train).copy()
    Y_test = np.exp(Y_test).copy()
    pred_tr_ts_sarimax = np.exp(pred_tr_ts_sarimax).copy()
    pred_te_ts_sarimax = np.exp(pred_te_ts_sarimax).copy()
    pred_te_ts_sarimax_ci = np.exp(pred_te_ts_sarimax_ci).copy()
if differencing: ## 차분한 것을 원래 값으로 변환
    pred_tr_ts_sarimax = np.cumsum(pred_tr_ts_sarimax).copy()
    
# 검증
%reload_ext autoreload
%autoreload 2
from module import *
Score_ts_sarimax, Resid_tr_ts_sarimax, Resid_te_ts_sarimax = evaluation_trte(Y_train, pred_tr_ts_sarimax, 
                                                                             Y_test, pred_te_ts_sarimax, graph_on=True)
display(Score_ts_sarimax)
ax = pd.DataFrame(Y_test).plot(figsize=(12,4))
## 점추정
pd.DataFrame(pred_te_ts_sarimax, index=Y_test.index, columns=['prediction']).plot(kind='line',
                                                                           xlim=(Y_test.index.min(),Y_test.index.max()),
                                                                           linewidth=3, fontsize=20, ax=ax)

## 구간추정
ax.fill_between(pd.DataFrame(pred_te_ts_sarimax_ci, index=Y_test.index).index,
                pd.DataFrame(pred_te_ts_sarimax_ci, index=Y_test.index).iloc[:,0],
                pd.DataFrame(pred_te_ts_sarimax_ci, index=Y_test.index).iloc[:,1], color='k', alpha=0.15)
plt.show()

• JB 테스트 결과 잔차는 정규분포에 가깝다고 본다.
• ljung-box 테스트 결과 자기상관이 없다고 한다
• 등분산성 테스트 결과 등분산이 아니라고 본다

- 예측과 실제값이 거의 비슷하게 나온다.

 

해당 포스팅은 패스트캠퍼스의 <파이썬을 활용한 시계열 데이터 분석 A-Z 올인원 패키지> 강의를 듣고 정리한 내용입니다