[시계열 머신러닝] 3. 단일 시계열 모델 분석 (ARIMA)

Kyung Jae, Cheong·2023년 3월 29일
데이터사이언스시계열 분석
5

시계열 데이터 분석 개념 정리

목록 보기
3/4
post-thumbnail

시계열 데이터 분석에 필요한 개념들을 요약하고 정리하는 시리즈 - 3. 단일 시계열 모델 분석 (ARIMA) 편

1. 시계열 데이터의 속성

  • 분석에 앞서서 시계열 데이터의 특징과 처리 방식을 간단하게 정리해봅니다.

1-1. 시계열 데이터의 정의

  • 데이터는 관측 시점과 변수에 따라 크게 2종류로 구분될 수 있습니다.
    • 횡단면(비시계열) 데이터 (Cross-Sectional)
      • 한 시점에 여러 변수에 대한 자료를 수집한 것
    • 시계열(종단면) 데이터 (Time-Series)
      • 시간의 흐름에 따라 한 변수에 대한 자료를 수집한 것
      • Longitudinal Data라고도 불립니다.
  • 횡단면 데이터와 시계열 데이터가 병합된 데이터를 Panel 데이터라고 부릅니다.
    • 관측치(변량) 2개 이상을 시점 2개 이상으로 누적
    • 일반적으로 패널 데이터들은 Truncation(누락)되거나 Censoring(검열)되어있는 경우가 많은데, 이를 Unbalanced Panel Data라고 부릅니다.

1-2. 시계열 데이터의 특징과 형태

  • 시계열 자료는 추세변동, 계절변동 등 크고 작은 다양한 변동요인들이 다중적으로 중첩되어 있습니다.
  • 시계열 데이터의 구성 요소
    • 계절성(Seasonal)
    • 추세성(Trend)
    • 순환성(Cycle)
    • 불규칙요소(Irregular, Residual)
  • 시계열 데이터의 형태는 크게는 다음과같이 분류됩니다.
    • 계절변동
    • 추세변동
    • 계절적추세변동
    • 순환변동
    • 우연변동 : 이를White Noise라고도 불립니다.

1-3. 백색잡읍(White Noise)

  • 백색잡음(White Noise)이란 패턴이 남아있지 않고 무작위로 야기되는 잡음을 의미합니다.
    • 시계열 데이터를 전처리한다는 것은 백색잡음으로 얻어내는 것을 의미한다 볼 수 있고, 이러한 백색잡음을 통해 모델에서 시계열 데이터 분석과 예측을 실시할 수 있게 됩니다.
  • 백색잡음은 2가지 속성을 만족해야하고, 하나라도 만족하지 않으면 모델에 대한 개선여지가 있음을 의미합니다.
    • 잔차들은 평균이 0인 정규분포이며, 일정한 분산을 가져야 한다.
    • 잔차들이 시간의 흐름에 따른 상관성이 없어야 한다.
      • 이러한 상관성을 자기상관성(Autocorrelation)이라 부릅니다.
      • 각 데이터 포인트가 각각 독립성을 지니는 것은 시계열 데이터의 기본적인 가정이기도 하며 이를 i.i.d. 라고도 부릅니다. (individually and independently distributed)
      • 이러한 상관성은 ACF, PACF등과 같은 함수들을 통해 확인해 볼 수 있으며 이에 대한 내용은 뒤에서 자세히 다룰 것입니다.
  • 시계열 예측 모델이 실제 현상의 트렌드 및 주기를 잘 반영할 수록 잔차의 변동이 작아지게 되고 이를 바탕으로 모델이 개선되었는지 여부를 파악할 수 있습니다.
    • 잔차 진단의 결과는 주로 시각화로 확인 할 수 있습니다.

1-4. 잔차 진단 (시계열 데이터 모델의 기본가정들)

  • 잔차에 대한 진단은 정상성, 정규분포, 자기상관성, 등분산성 총 4가지를 검정하게 되며, 이러한 특성은 시계열 데이터 모델의 기본가정들이기도 합니다.
  • 정상성(Stationarity)
    • 시간에 관계없이 일정한 상태를 유지하는 것
      • 자기상관성과 등분산성을 아우르는 개념입니다.
    • 내재적인 패턴의 존재를 의미하는 단위근(Unit Root) 존재여부를 검정합니다.
    • 대표적으로는 ADF 검정(Augemented Dickey-Fuller Test), KPSS 검정(Kwiatkowski-Phillips-Schmidt-Shin Test)이 있고, 이 둘의 귀무가설을 서로 반대임을 주의해야합니다.
      • ADF 검정의 귀무가설(H0H_0) : 단위근(Unit Root) 존재한다 (비정상)
      • KPSS 검정의 귀무가설(H0H_0) : 단위근(Unit Root)이 없다 (정상)
  • 정규분포(Normal distribution)
    • 귀무가설(H0H_0) : 정규분포를 따른다.
    • 대립가설(HaH_a) : 정규분포를 따르지 않는다.
    • Shapiro-Wilk test, Kolmogorov-Smirnov test 등 다양한 방법들이 있습니다.
  • 자기상관성(Autocorrelation)
    • 귀무가설(H0H_0) : 시간이 지나면 autocorrelation 은 0이다. 자기상관관계가 없다.
    • 대립가설(HaH_a) : 시간이 지나면 autocorrelation 은 0이 아니다. 자기상관관계가 있다.
    • Ljung-Box test, Portmanteau test 등 다양한 방법들이 있습니다.
  • 등분산성(Homoscedasticity)
    • 귀무가설(H0H_0) : 시간이 지나면 등분산이다.
    • 대립가설(HaH_a) : 시간이 지나면 등분산이 아니다. 발산하는 분산이다.
    • Goldfeld-Quandt test, Breusch-Pagan test 등 다양한 방법들이 있습니다.

1-5. 시계열 데이터의 처리

  • 시계열 데이터는 추세변동, 계절변동 등 크고 작은 다양한 변동요인들이 다중적으로 중첩되어 있습니다.
  • 이러한 시계열 데이터의 분석을 수행할 때에는 크게 2가지의 접근법을 이용할 수 있습니다.
    • 요소분해법(Decomposision) : 내포된 변동요인이 고정적인 패턴을 보이는 경우
      • 가법모형(Additive)
        • yt=St+Tt+Ct+Rty_t = S_t+T_t+C_t+R_t
        • 선형적이고, Trend, Seasonal, Cycle이 서로 독립적이라 가정할 수 있는 경우에 이용됩니다.
      • 승법모형(Multiplicative)
        • yt=St∗Tt∗Ct∗Rty_t = S_t*T_t*C_t*R_t
        • 비선형적이고, Trend, Seasonal, Cycle이 서로 영향을 주고받는 경우에 이용됩니다.
    • 평활법(Smooting) : 다양한 변동요인이 고정적인 패턴을 보이지 않는 경우
      • 이동평균 평활법(MA Smoothing)
        • 과거로부터 현재까지의 시계열 자료를 대상으로 일정 기간별 이동평균을 계산하고 이들의 추세를 통해 다음 기간을 예측
      • 지수평균 평활법(Exponential Smooting)
        • 모든 시계열 자료를 사용하여 평균을 구하고, 시간의 흐름에 따라서 최근 시계열에 더 많은 가중치를 부여하여 미래의 값을 예측
  • 시계열 데이터 처리 방법들
    • 변동폭이 일정하지 않은 경우 : log 변환
    • 추세, 계절성이 존재하는 경우 : 차분(differencing)

2. 시계열 확률모형(Stochastic Model)

  • 분석 모형에 따라 추정값과 관측값의 차이인 잔차(Residual)의 시계열 자료(Residual Error Series)는 확률 과정을 통해 나타낼 수 있습니다.
    • 이때 고려되는 모델이 확률모형(Stochastic Model)입니다.
    • 시계열 분석 모델이 실제 현상과 적합한 경우에 시계열 자료는 백색잡음(White Noise)이 된다는 가정을 바탕으로 합니다.
    • xt=yt−y^tx_t = y_t-\hat{y}_t = ~ N(0,σ2)N(0,\sigma^2)
    • 즉, 잔차는 서로 독립적인 경우에 시계열 분석 모델이 실제 현상을 잘 설명하고 있다고 해석할 수 있음을 의미합니다.
  • 확률모형의 분석은 다음의 4가지 단계를 걸쳐 수행됩니다.
    • 모델 설정(Model Identification)
      • ACF, PACF 같은 지표와 ARMA, ARIMA와 같은 프로그램을 이용해 적합한 모형을 선택합니다.
    • 모수 추정(Parameter Estimation)
      • 설정한 모델에 따라 모수(계수)를 추청하고 적합성을 검토합니다.
    • 통계적 검정(Statistical Testing)
      • 모델의 적합도에 따른 잔차의 정상성과 분석결과에 대한 검정을 수행합니다.
    • 예측(Forecasting)
      • 시계열 분석 모델을 통해 미래의 값을 예측합니다.

2-1. 단일 변량 시계열 모델

  • 하나의 변수에 대한 시계열 데이터를 분석하는 모델은 기본적으로 AR, MA, ARMA, ARIMA 모델이 있습니다.

AR 모델 (AutoRegressive Model)

  • 자기상관성을 시계열 모형으로 구성한 것으로, 변수의 과거 관측값의 선형 결합을 통해 변수의 미래값을 예측하는 모델입니다.
  • 이전의 관측값이 이후의 관측값에 영향을 준다는 아이디어에 기반하고 있습니다.
  • 파라미터 값으로는 p를 이용하며 AR(p)의 형태로 나타낼 수 있습니다.
  • yt=c+ϕ1yt−1+ϕ2yt−2+⋯+ϕpyt−p+ϵty_t = c + \phi_1y_{t-1} + \phi_2y_{t-2} + \cdots + \phi_py_{t-p} + \epsilon_t
  • yty_t는 t시점의 관측값, c는 상수, ϕ\phi는 가중치, ϵt\epsilon_t는 오차항을 의미합니다.
    • 오차항은 평균이 0, 분산이 1인 표준정규분포를 따르는 백색잡음입니다.

MA 모델 (Moving Average Model)

  • 예측 오차를 이용하여 미래값을 예측하는 모델입니다.
  • 파라미터 값으로는 q를 이용하며 MA(q)의 형태로 나타낼 수 있습니다.
  • yt=c+θ1ϵt−1+θ2ϵt−2+⋯+θqϵt−q+ϵty_t=c+\theta_1 \epsilon_{t-1} + \theta_2 \epsilon_{t-2} + \cdots + \theta_q \epsilon_{t-q} + \epsilon_t

ARMA 모델 (안정시계열 모델)

  • p개의 자기 자신의 과거값과 q개의 과거백색 잡음의 선형조합을 의미하고, AR(p) 모형과 MA(q)모형이 합쳐진 모델입니다.
  • AR(p)모형의 정상성, MA(q)모형의 가역성 조건을 만족해야 합니다.
  • 파라미터 값으로는 p,q를 이용하며, ARMA(p,q)의 형태로 나타낼 수 있습니다.
  • yt=AR(p)+MA(q)y_t=AR(p)+MA(q)
  • yt=c+ϕ1yt−1+ϕ2yt−2+⋯+ϕpyt−p+ϵt+θ1ϵt−1+θ2ϵt−2+⋯+θqϵt−qy_t = c + \phi_1y_{t-1} + \phi_2y_{t-2} + \cdots + \phi_py_{t-p} + \epsilon_t + \theta_1 \epsilon_{t-1} + \theta_2 \epsilon_{t-2} + \cdots + \theta_q \epsilon_{t-q}
  • ARMA는 값을 0을 주게 되면 AR, MA 모델과 동일해집니다.
    • ARMA(1, 0) = AR(1)
    • ARMA(0, 1) = MA(1)
  • ARMA(1,1)은 가장 흔한 형태의 정상시계열 모델입니다.
    • 성능 벤치마크로써 베이스 모델의 역할도 수행할 수 있습니다.

ARIMA 모델 (불안정 시계열 모델)

  • ARMA 모델에서 d차 차분을 추가적으로 적용시킨 모델로 비정상시계열에도 바로 적용할 수 있다는 장점이 있습니다.
  • ARIMA(p,d,q) : d차 차분한 데이터에 AR(p)모형과 MA(q)모형을 합친 모델
  • I는 Integrated를 뜻하고 이는 '누적'을 의미합니다.
  • ARIMA 모델 -> AR, MA
    • ARIMA(p,0,0) = AR(p)
    • ARIMA(0,0,q) = MA(q)
  • p,d,q 값은 AIC, BIC, ACF, PACF등 다양한 방법을 통해서 선택할 수 있습니다.
    • 자기상관(p,q)은 ACF, PACF 그래프를 그려서 확인합니다.

2-2. 자기상관함수(ACF, PACF)

ACF(AutoCorrelation Function)

  • 시차(lag)에 따른 연속적인 관측값 사이의 상호 연관관계를 의미하며, 시차가 커질 수록 0에 수렴하게 됩니다.
  • 동일한 변수를 시점을 달리하여 관측했을 때 시점에 따라 다른값 사이의 상호 연관관계를 척도로 나타낸 것입니다.
  • ACF를 통해서 시계열 자료의 정상성을 파악할 수 있습니다.
  • 정상시계열인 경우엔 상대적으로 빠르게 0에 수렴하고, 비정상시계열은 천천히 감소합니다.
  • yty_t와 yt+ky_{t+k}의 자기상관을 구하는 식
  • MA, ARMA, ARIMA에서 q값을 구하는 경우에 ACF 그래프를 통해서 구할 수 있으며, 0으로 수렴하게되는 시점을 q값으로 잡으면 됩니다.
    • 정확하게는 0을 중심으로하는 신뢰구간(95%에선 0±1.96n0\pm\frac{1.96}{\sqrt{n}})에 진입하는 최초의 시점을 q값으로 설정합니다.
  • 시계열 데이터가 MA의 특성을 띄는 경우 ACF는 급격하게 감소합니다.

PACF(Partial ACF)

  • 부분자기상관함수
  • 시차(lag)에 따른 편자기상관을 의미하며, 시차가 다른 두 시계열간의 순수한 상호연관관계를 의미합니다.
  • yty_t와 yt+ky_{t+k}의 순수한 상관관계로 두 시점 사이에 포함된 모든 값(yt−1,yt2,⋯ ,yt−k+1y_{t-1},y_{t_2},\cdots,y_{t-k+1})들의 영향을 제거하여 상관관계를 정의합니다.
  • yty_t와 yt+ky_{t+k}의 편자기상관을 구하는 식
  • AR, ARMA, ARIMA에서 p값을 구하는 경우에 PACF 그래프를 통해서 구할 수 있으며, 0으로 수혐하게 되는 시점을 p값으로 잡으면 됩니다.
    • 정확하게는 0을 중심으로하는 신뢰구간(95%에선 0±1.96n0\pm\frac{1.96}{\sqrt{n}})에 진입하는 최초의 시점을 p값으로 설정합니다.
  • 시계열 데이터가 AR의 특성을 띄는 경우 PACF는 급격하게 감소합니다.

3. Python 실습 (ARIMA)

  • Bindukuri 지역의 일간 평균기온 시계열 데이터를 통한 python 실습입니다.
# 라이브러리 불러오기
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns

데이터 불러오기

# 데이터 불러오기
df0 = pd.read_csv('Bindukuri_temperature.csv')
print(df0.shape)
display(df0.head(5))
display(df0.tail(5))
  • 데이터 확인하기 (결측치여부, column 타입)
df0.info()
'''
<class 'pandas.core.frame.DataFrame'>
RangeIndex: 751 entries, 0 to 750
Data columns (total 2 columns):
 #   Column    Non-Null Count  Dtype  
---  ------    --------------  -----  
 0   Date      751 non-null    object 
 1   MeanTemp  751 non-null    float64
dtypes: float64(1), object(1)
memory usage: 11.9+ KB
'''

EDA(탐색적 데이터 분석)

  • 'Date' 열을 datetime형으로 변환하기
ts0 = df0.copy()
ts0['Date'] = pd.to_datetime(ts0.Date, format='%Y-%m-%d')
ts0.info()
'''
<class 'pandas.core.frame.DataFrame'>
RangeIndex: 751 entries, 0 to 750
Data columns (total 2 columns):
 #   Column    Non-Null Count  Dtype         
---  ------    --------------  -----         
 0   Date      751 non-null    datetime64[ns]
 1   MeanTemp  751 non-null    float64       
dtypes: datetime64[ns](1), float64(1)
memory usage: 11.9 KB
'''
  • (번외) 연,월,일,요일 정보를 추출해보기
    • 필요에 따라서 분석을 진행하고자 할 때 유용하게 쓰일 수 있습니다.
ts1 = ts0.copy()
# 연도
ts1['Year'] = ts1.Date.dt.year
# 월
ts1['Month'] = ts1.Date.dt.month
# 일
ts1['Day'] = ts1.Date.dt.day
# 요일
ts1['Day_name'] = ts1.Date.dt.day_name()
ts1
  • 필요한 Column만 남기고 날짜를 인덱스로 변환하기
# 필요한 column만 남기기
ts = ts1.loc[:, ['Date', 'MeanTemp']]
# 인덱스 설정 후 drop
ts.index = ts.Date
ts = ts.drop(columns="Date")
ts
  • 시각화 함수 정의 및 분포 시각화
# 시각화 함수 정의
def plot_ts(data, color, alpha, label):
    
    plt.figure(figsize=(11,5))
    plt.plot(data, color=color, alpha=alpha, label=label)
    plt.title("Mean Temperature of Bindukuri")
    plt.ylabel('Mean temperature')
    plt.legend()
    plt.show()
plot_ts(ts, 'blue', 0.25, 'Original')

시계열 데이터의 정상성 검정: ADF

  • statsmodels의 adfuller 함수를 통해 검정 실시
    • 출력값 순서
      • ADF 통계량
      • p-value
      • usedlag (The number of lags used)
      • nobs (The number of observations used for the ADF regression and calculation of the critical values.)
      • critical values (기각역)
      • icbest (The maximized information criterion if autolag is not None.) 
from statsmodels.tsa.stattools import adfuller

adfuller(ts, autolag='AIC')

'''
(-1.4095966745887758,
 0.5776668028526356,
 11,
 739,
 {'1%': -3.439229783394421,
  '5%': -2.86545894814762,
  '10%': -2.5688568756191392},
 2474.854578321788)
'''
  • ADF 결과 중 통계량, p-value, 기각역을 계산하고 출력하는 함수 정의 및 실행
def ADF_test(data):
	# ADF 실시
    results = adfuller(data, autolag='AIC')
    
    # 통계량
    s = results[0]
    # p-value
    p = results[1]
    # 기각역
    cv = results[4]
    
    # 출력
    print('-'*30)
    print('Augemented Dickey-Fuller Test')
    print('H0 : 단위근이 존재한다 (비정상 시계열)')
    print('Ha : 단위근이 없다 (정상 시계열)')
    print('Critical Values : {}'.format(cv))
    print('-'*30)
    print('Test Statistics : {:.4f}'.format(s))
    print('p-value : {:.4f}'.format(p))
    print('-'*30)

# ts 데이터로 ADF 실행해보기
ADF_test(ts)

'''
------------------------------
Augemented Dickey-Fuller Test
H0 : 단위근이 존재한다 (비정상 시계열)
Ha : 단위근이 없다 (정상 시계열)
Critical Values : {'1%': -3.439229783394421, '5%': -2.86545894814762, '10%': -2.5688568756191392}
------------------------------
Test Statistics : -1.4096
p-value : 0.5777
------------------------------
'''
# 영가설 기각 불가 : 단위근이 존재함. 비정상 시계열이다.
  • 이동 평균 함수를 이용한 평균 및 표준편차 분포 시각화
def plot_rolling(data, roll_size):
    # 이동평균함수(rolling) - 평균, 표준편차
    roll_mean = data.rolling(window=roll_size).mean()
    roll_std = data.rolling(window=roll_size).std()
    
    # 시각화
    plt.figure(figsize=(11,5))
    plt.plot(data, color='blue', alpha=0.25, label='Original')
    plt.plot(roll_mean, color='black', label='Rolling mean')
    plt.plot(roll_std, color='red', label='Rolling std')
    plt.title('Rolling Mean & Standard Deviation')
    plt.ylabel("Mean Temperature")
    plt.legend()
    plt.show()

# 함수 실행
plot_rolling(ts, 6)

차분(Differencing)으로 정상시계열 만들기

  • 방법1 : .shift()함수를 이용
# ts에서 ts.shift(1)을 빼기 (1차 차분)
ts_diff = ts - ts.shift() # default=1임

# 시각화
plot_ts(ts_diff, color='blue', alpha=0.25, label='differencing')
  • 이동평균 시각화 및 ADF 테스트 실시
plot_rolling(ts_diff, 6)
ADF_test(ts_diff.dropna())
  • '''
------------------------------
Augemented Dickey-Fuller Test
H0 : 단위근이 존재한다 (비정상 시계열)
Ha : 단위근이 없다 (정상 시계열)
Critical Values : {'1%': -3.439229783394421, '5%': -2.86545894814762, '10%': -2.5688568756191392}
------------------------------
Test Statistics : -11.6790
p-value : 0.0000
------------------------------
'''
  • 방법2 : .diff()함수 이용
    • 이게 더 간단하고 결측치 처리도 한번에 실시할 수도 있음
# 차분
ts_diff2 = ts.diff().dropna()

# ADF 테스트
ADF_test(ts_diff2)
'''
------------------------------
Augemented Dickey-Fuller Test
H0 : 단위근이 존재한다 (비정상 시계열)
Ha : 단위근이 없다 (정상 시계열)
Critical Values : {'1%': -3.439229783394421, '5%': -2.86545894814762, '10%': -2.5688568756191392}
------------------------------
Test Statistics : -11.6790
p-value : 0.0000
------------------------------
'''
# 영가설 기각 : 단위근이 없다. 정상시계열이다.

ARIMA로 평년 기온 예측하기

  • ACF, PACF 실시 및 시각화
# ACF and PACF 
from statsmodels.tsa.stattools import acf, pacf

# ACF
acf_20 = acf(x=ts_diff2, nlags=20)
# PACF
pacf_20 = pacf(x=ts_diff2, nlags=20, method='ols')

# 95% 신뢰구간 계산하기
confidence = 1.96/np.sqrt(len(ts_diff2))

# 시각화
plt.figure(figsize=(8,4))
# ACF
plt.subplot(1,2,1)
plt.plot(acf_20)
plt.axhline(y=0, linestyle='--', alpha=0.5)
plt.axhline(y=-confidence, linestyle='--', alpha=0.5)
plt.axhline(y=confidence, linestyle='--', alpha=0.5)
plt.title('AutoCorrelation Function')
# PACF
plt.subplot(1,2,2)
plt.plot(pacf_20)
plt.axhline(y=0, linestyle='--', alpha=0.5)
plt.axhline(y=-confidence, linestyle='--', alpha=0.5)
plt.axhline(y=confidence, linestyle='--', alpha=0.5)
plt.title('Partial AutoCorrelation Function')

plt.tight_layout()
  • lag를 20에서 10정도로 줄여서 봐도 무방할 듯
# ACF
acf_10 = acf(x=ts_diff2, nlags=10)
# PACF
pacf_10 = pacf(x=ts_diff2, nlags=10, method='ols')

# 95% 신뢰구간
confidence = 1.96/np.sqrt(len(ts_diff2))

# 시각화
plt.figure(figsize=(8,4))

plt.subplot(1,2,1)
plt.plot(acf_10)
plt.axhline(y=0, linestyle='--', alpha=0.5)
plt.axhline(y=-confidence, linestyle='--', alpha=0.5)
plt.axhline(y=confidence, linestyle='--', alpha=0.5)
plt.title('AutoCorrelation Function')

plt.subplot(1,2,2)
plt.plot(pacf_10)
plt.axhline(y=0, linestyle='--', alpha=0.5)
plt.axhline(y=-confidence, linestyle='--', alpha=0.5)
plt.axhline(y=confidence, linestyle='--', alpha=0.5)
plt.title('Partial AutoCorrelation Function')

plt.tight_layout()
# ARIMA(p, d, q)
# 최적 q값 (신뢰구간 최초 진입 시점)
	# ACF (q, MA) = 3
# 최적 p값 (신뢰구간 최초 진입 시점)
	# PACF (p, AR) = 6 

# 예상되는 최적의 ARIMA 모델
	# ARIMA(6, 1, 3)
  • ARMA(1,1) = ARIMA(1,0,1) 실시해보기 (벤치마크)
# ARIMA
from statsmodels.tsa.arima.model import ARIMA

# index를 period로 변환해주어야 warning이 뜨지 않음
ts_copy = ts.copy()
ts_copy.index = pd.DatetimeIndex(ts_copy.index).to_period('D')

# 예측을 시작할 위치(이후 차분을 적용하기 때문에 맞추어주었음
start_idx = ts_copy.index[1]

# ARIMA(1,0,1)
model1 = ARIMA(ts_copy, order=(1,0,1))
# fit model
model1_fit = model1.fit()

# 전체에 대한 예측 실시
forecast1 = model1_fit.predict(start=start_idx)
  • 예측 시각화 및 오차함수(MSE) 함수 정의 및 실시
from sklearn.metrics import mean_squared_error

def plot_and_error(data, forecast):
    # MSE 계산
    mse = mean_squared_error(data, forecast)
    # 시각화
    plt.figure(figsize=(11,5))
    plt.plot(data, color='blue', alpha=0.25 , label='Original')
    plt.plot(forecast, color='red', label='predicted')
    plt.title("Time Series Forecast")
    plt.ylabel("Mean Temperature")
    plt.legend()
    plt.show()
    # MSE 출력
    print('Mean Squared Error : {:.4f}'.format(mse))

plot_and_error(ts[1:], forecast1)
  • 최적화 파라미터로 ARIMA실시해보기
# 모델 파라미터 최적화 (p=6, d=1, q=3)
model2 = ARIMA(ts_copy, order=(6,1,3))
# fit
model2_fit = model2.fit()
# 예측
forecast2 = model2_fit.predict(start=start_idx)
# 시각화 및 MSE 연산
plot_and_error(ts[1:], forecast2)
  • MSE가 상대적으로 줄어들었으므로(1.8389 -> 1.7304)성능의 개선이 이루어 졌음을 확인 가능!

4. Python 실습 2 - Auto ARIMA

  • 현재까지는 ARIMA의 p,d,q값을 직접적으로 설정하는 방식을 통해 모델링과 최적화를 진행하였지만, ACF나 PACF의 경우엔 다소 주관적으로 판단해야하는 경우가 있었습니다.
  • Auto ARIMA는 이러한 파라미터의 탐색과정을 자동화하여 보다 쉽고 정확하게 가장 높은 성능을 이끌어 낼 수 있도록 최적화할 수 있도록 하는 Tool이며, Auto ARIMA를 지원하는 다양한 라이브러리들이 있습니다.
    • 보통 머신러닝에서는 Grid Search나 Bayesian Search를 통해 하이퍼파라미터를 최적화하는데, Auto ARIMA는 Bayesian Search의 방식에 좀 더 가까운 방식으로 동작합니다.
    • stepwise algorithm : Hyndman and Khandakar (2008)의 알고리즘으로 가능한 모든 조합을 찾는 것이 아닌 주어진 상황에서 확률적으로 파라미터들을 바꾸어보면서 최적화를 진행하는 방식
  • 본 실습에서는 pmdarima라는 파이썬 라이브러리를 통해 Auto ARIMA를 실시하였습니다.
    • 2017년의 주식정보를 통해 2018년도의 주식가격을 예측해보고, 실제값과 얼만큼의 차이가 나는지를 확인해보는 것이 목표!
# 기본 라이브러리 불러오기
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns

데이터 불러오기

  • 어느정도 가공을 거친 데이터를 이용하였고, 기업명은 익명으로 처리하였습니다.
# 데이터 불러오기 (첫 열을 datetime형태로 인덱스로 지정하여 불러오기)
df0 = pd.read_csv('stock_data.csv', index_col=0, parse_dates=True)
df0.info()
'''
<class 'pandas.core.frame.DataFrame'>
DatetimeIndex: 312 entries, 2017-01-03 to 2018-03-29
Data columns (total 2 columns):
 #   Column     Non-Null Count  Dtype  
---  ------     --------------  -----  
 0   company_A  312 non-null    float64
 1   company_M  312 non-null    float64
dtypes: float64(2)
memory usage: 7.3 KB
'''
  • 불러온 데이터 확인
print(df0.shape)
display(df0.head())
display(df0.tail())
  • 결측치 확인
df0.isnull().sum()
'''
company_A    0
company_M    0
dtype: int64
'''
  • 필요한 column만 남기기 (A회사 데이터만 남기기)
stock_a = df0[['company_A']]
stock_a

데이터 분리

  • 앞의 80% 데이터는 train 데이터로, 뒤의 20% 데이터는 test 데이터로 나누기
    • 시계열 데이터는 시간 변수에 따른 자기상관성이 있기 때문에 함부로 shuffle을 통한 분리를 실시할 수 없는 것이 특징입니다.
train = stock_a['company_A'][:int(0.8*len(stock_a))]
test = stock_a['company_A'][int(0.8*len(stock_a)):]
  • 시각화로 확인하기
train.plot(label='Train')
test.plot(label='Test')
plt.legend()
plt.show()
  • 트렌드와 주기성이 있는 것으로 보이므로 차분을 통해 정상성 데이터로 만들어주어야 할 듯 합니다.

최적 차분 수 구하기

  • pmdarima라이브러리를 통해 최적의 차분 수를 구함
  • KPSS와 ADF모두 지원하므로 둘 모두 시행하여 최대값으로 ARIMA의 d값을 결정
import pmdarima as pm

kpss_diffs = pm.arima.ndiffs(train, alpha=0.05, test='kpss', max_d=5)
adf_diffs = pm.arima.ndiffs(train, alpha=0.05, test='adf', max_d=5)
n_diffs = max(kpss_diffs, adf_diffs)

print(f"Optimized 'd' = {n_diffs}")

'''
Optimized 'd' = 1
'''
  • 1차 차분만으로도 정상성을 만족할 수 있으므로 d값은 1로 설정

Auto ARIMA를 통한 모델링

# 직접 지정해주어 모델링 실시
model = pm.auto_arima(y=train,		# 데이터
                      d=n_diffs,	# 차분 (d), 기본값 = None
                      start_p= 0,	# 시작 p값, 기본값 = 2
                      max_p = 5,	# p 최대값, 기본값 = 5
                      start_q= 0,	# 시작 q값, 기본값 = 2
                      max_q = 5,	# q 최대값, 기본값 = 5
                      m=1,			# season의 주기, 기본값 = 1
                      seasonal=False,	# sARIMA를 실시, 기본값 = True
                      stepwise=True,	# stepwise algorithm, 기본값 = True
                      trace=True)		# 각 step을 출력할지, 기본값 = False

'''
Performing stepwise search to minimize aic
 ARIMA(0,1,0)(0,0,0)[0] intercept   : AIC=967.793, Time=0.01 sec
 ARIMA(1,1,0)(0,0,0)[0] intercept   : AIC=969.462, Time=0.03 sec
 ARIMA(0,1,1)(0,0,0)[0] intercept   : AIC=969.432, Time=0.07 sec
 ARIMA(0,1,0)(0,0,0)[0]             : AIC=969.948, Time=0.01 sec
 ARIMA(1,1,1)(0,0,0)[0] intercept   : AIC=971.366, Time=0.14 sec

Best model:  ARIMA(0,1,0)(0,0,0)[0] intercept
Total fit time: 0.266 seconds
'''
  • 위에서는 0부터 지정하여 실시하였지만, 필요한 파라미터를 제외하고 모든 값을 기본값으로 설정하여 모델링을 다시 진행해보았습니다.
model2 = pm.auto_arima(train, d=1, seasonal=False, trace=True)

'''
Performing stepwise search to minimize aic
 ARIMA(2,1,2)(0,0,0)[0] intercept   : AIC=966.160, Time=0.33 sec
 ARIMA(0,1,0)(0,0,0)[0] intercept   : AIC=967.793, Time=0.01 sec
 ARIMA(1,1,0)(0,0,0)[0] intercept   : AIC=969.462, Time=0.03 sec
 ARIMA(0,1,1)(0,0,0)[0] intercept   : AIC=969.432, Time=0.04 sec
 ARIMA(0,1,0)(0,0,0)[0]             : AIC=969.948, Time=0.01 sec
 ARIMA(1,1,2)(0,0,0)[0] intercept   : AIC=972.919, Time=0.10 sec
 ARIMA(2,1,1)(0,0,0)[0] intercept   : AIC=972.738, Time=0.13 sec
 ARIMA(3,1,2)(0,0,0)[0] intercept   : AIC=969.483, Time=0.52 sec
 ARIMA(2,1,3)(0,0,0)[0] intercept   : AIC=971.491, Time=0.22 sec
 ARIMA(1,1,1)(0,0,0)[0] intercept   : AIC=971.366, Time=0.14 sec
 ARIMA(1,1,3)(0,0,0)[0] intercept   : AIC=972.667, Time=0.13 sec
 ARIMA(3,1,1)(0,0,0)[0] intercept   : AIC=973.267, Time=0.13 sec
 ARIMA(3,1,3)(0,0,0)[0] intercept   : AIC=966.431, Time=0.57 sec
 ARIMA(2,1,2)(0,0,0)[0]             : AIC=968.841, Time=0.20 sec

Best model:  ARIMA(2,1,2)(0,0,0)[0] intercept
Total fit time: 2.567 seconds
'''
  • 초기값이 달라지는 경우엔 최적화된 값을 찾아가는 과정과 결과가 달라짐을 알 수 있고, 이는 stepwise algorithm을 따르기 때문에 발생하는 차이입니다.
    • 두번째 모델의 AIC 점수가 더 낮게 나왔으므로 ARIMA(2,1,2)를 최종적인 모델로 선정하고 예측을 실시해보도록 하겠음

최적화된 ARIMA 모델 분석

model2.summary()
  • 시각화를 통한 정상성 확인
    • pmdarima에서 제공하는 plot_diagnostics 기능을 이용
    • Correlogram은 ACF 그래프를 의미
model2.plot_diagnostics(figsize=(16,8))
plt.show()

예측 및 평가

  • 훈련에는 적용되지 않은 테스트데이터를 통해 예측을 실시해봅니다.
test.shape
'''
(63,)
'''
  • 한번에 예측을 진행해봅니다 (트렌드만을 고려한 예측이 실시됩니다)
# 예측 -> 리스트로 변환
pred = model2.predict(n_periods=len(test)).to_list()
# 데이터프레임 생성
test_pred = pd.DataFrame({'test':test, 'pred':pred}, index=test.index)
test_pred
  • 시각화로 확인해보기
plt.plot(train, label='Train')
plt.plot(test, label='Test')
plt.plot(test_pred.pred, label='Predicted')
plt.legend()
plt.show()
  • 이번에는 한 지점에 대한 예측을 진행하고 모델을 업데이트 하는 방식으로 예측을 진행해보도록 하겠습니다. (트렌드 이외에도 변동요인들이 모두 반영됩니다)
# one point forcast 함수 정의, 신뢰구간도 함께 담아보기
def forcast_one_step():
    fc, conf = model2.predict(n_periods=1, return_conf_int=True)
    return fc.tolist()[0], np.asarray(conf).tolist()[0]

# 값들을 담을 빈 리스트를 생성
y_pred = []
pred_upper = []
pred_lower = []

# for문을 통한 예측 및 모델 업데이트를 반복함
for new_ob in test:
    fc, conf = forcast_one_step()
    y_pred.append(fc)
    pred_upper.append(conf[1])
    pred_lower.append(conf[0])
    model2.update(new_ob)
  • 예측 결과를 데이터프레임으로 만들기
test_pred2 = pd.DataFrame({'test':test, 'pred':y_pred})
y_pred_df = test_pred2['pred']	# Series로 반환
y_pred_df
'''
Date
2017-12-28    170.843675
2017-12-29    171.682863
2018-01-02    169.646252
2018-01-03    172.295839
2018-01-04    172.064889
                 ...    
2018-03-23    169.268631
2018-03-26    165.184800
2018-03-27    172.805556
2018-03-28    168.247672
2018-03-29    166.840596
Name: pred, Length: 63, dtype: float64
'''
  • 시각화로 확인
plt.plot(train, label='Train')
plt.plot(test, label='Test')
plt.plot(y_pred_df, label='Predicted')
plt.legend()
plt.show()
  • 업데이트된 모델 분석
model2.summary()
  • 예측 모델 오차 계산 (MAPE)
# sklearn으로 MAPE 계산
from sklearn.metrics import mean_absolute_percentage_error
print(f"MAPE : {mean_absolute_percentage_error(test, y_pred):.3f}")
'''
MAPE : 0.013
'''

# numpy로 직접 계산
def MAPE(y_test, y_pred): 
	return np.mean(np.abs((test - y_pred) / y_test)) 
print(f"MAPE : {MAPE(test, y_pred):.3f}")
'''
MAPE : 0.013
'''

마무리

  • 이번 글에서는 시계열 데이터의 속성과 처리, 시계열 확률모형을 살펴보고 python 실습을 진행해 보았습니다.

Keyword 요약

    1. 시계열 데이터의 속성
    • 1-1. 시계열 데이터의 정의
      • 횡단면(비시계열) 데이터 (Cross-Sectional)
      • 시계열(종단면) 데이터 (Time-Series)
    • 1-2. 시계열 데이터의 특징과 형태
      • 계절성(Seasonal)
      • 추세성(Trend)
      • 순환성(Cycle)
      • 불규칙요소(Irregular, Residual)
    • 1-3. 백색잡읍(White Noise)
      • 우연변동 시계열
      • White Noise 조건
        • 잔차들은 평균이 0인 정규분포이며, 일정한 분산을 가져야 한다.
        • 잔차들이 시간의 흐름에 따른 상관성이 없어야 한다.
      • 자기상관성(Autocorrelation)
      • i.i.d.
    • 1-4. 잔차 진단 (시계열 데이터 모델의 기본가정들)
      • 정상성(Stationarity)
        • 단위근(Unit Root)
        • ADF 검정(Augemented Dickey-Fuller Test)
        • KPSS 검정(Kwiatkowski-Phillips-Schmidt-Shin Test)
      • 정규분포(Normal distribution)
      • 자기상관성(Autocorrelation)
      • 등분산성(Homoscedasticity)
    • 1-5. 시계열 데이터의 처리
      • 요소분해법(Decomposision)
        • 가법모형(Additive)
        • 승법모형(Multiplicative)
      • 평활법(Smooting)
        • 이동평균 평활법(MA Smoothing)
        • 지수평균 평활법(Exponential Smooting)
      • 변동폭이 일정하지 않은 경우 : log 변환
      • 추세, 계절성이 존재하는 경우 : 차분(differencing)
    1. 시계열 확률모형(Stochastic Model)
    • 2-1. 단일 변량 시계열 모델
      • AR 모델 (AutoRegressive Model)
        • AR(p)
      • MA 모델 (Moving Average Model)
        • MA(q)
      • ARMA 모델 (안정시계열 모델)
        • ARMA(p,q)
      • ARIMA 모델 (불안정 시계열 모델)
        • ARIMA(p,d,q)
    • 2-2. 자기상관함수(ACF, PACF)
      • ACF(AutoCorrelation Function)
        • MA, ARMA, ARIMA에서 q값을 구하는 경우
      • PACF(Partial ACF)
        • AR, ARMA, ARIMA에서 p값을 구하는 경우

  • 다음 글에서는 다변량 시계열 분석 기법에 대해 다루어보고, 이후엔 딥러닝 기법들도 다루어볼 예정입니다.

profile
Kyung Jae, Cheong
Machine Learning (AI) Engineer & BackEnd Engineer (Entry)
이전 포스트

[시계열 머신러닝] 2. 머신러닝과 지도학습 모델

다음 포스트

[시계열 머신러닝] 4. 다변량 시계열 모델 (ARCH, VAR, VECM)

4개의 댓글

comment-user-thumbnail
김수희
2024년 3월 21일

안녕하세요!
글 잘 보았습니다.

관련해서 질문이 있어 댓글 남깁니다.
test 데이터는 우리가 현재 예측 시점에서 알지 못하는 정보라고 생각하는데요

Ex) 향후 16주를 예측한다고 했을 때, 다음 주 ~ 16주 이후 모두 모르는 값.

아래의 모델 업데이트는 실제 test 데이터 값인 new_ob를 사용해서 이루어지는 것이
맞는 방법인지 의문이 들더군요.
오히려 모델의 예측 값으로 업데이트를 해야하는 것이 아닌가 라는 생각도 들구요

혹시 관련해서 제가 잘못 이해하고 있는 부분이 있거나, 작성자님의 의견 말씀해주시면 감사하겠습니다!

for문을 통한 예측 및 모델 업데이트를 반복함

for new_ob in test:
fc, conf = forcast_one_step()
y_pred.append(fc)
pred_upper.append(conf[1])
pred_lower.append(conf[0])
model2.update(new_ob)

1개의 답글