🔥알림🔥
① 테디노트 유튜브 - 구경하러 가기!
② LangChain 한국어 튜토리얼 바로가기 👀
③ 랭체인 노트 무료 전자책(wikidocs) 바로가기 🙌
④ RAG 비법노트 LangChain 강의오픈 바로가기 🙌
⑤ 서울대 PyTorch 딥러닝 강의 바로가기 🙌

18 분 소요

시계열 데이터 분석이라는 다소 넓은 범위의 주제이지만, 그 중에서도 전통적인 시계열 분석에서 언급되는 주요 분석 기법에 대해서 간략히 알아보고, 파이썬(Python) 코드로 구현하는 방법에 대해 알아보겠습니다.

필요한 모듈 import

시계열 데이터 분석에 필요한 패키지를 설치합니다.

# !pip install pmdarima -q
# !pip install statsmodels -q

import os
import random
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
import seaborn as sns
import warnings

# SEED 설정
SEED = 123

#####################    Pandas   #####################

# DataFrame 행과 열의 출력 개수 지정. None은 전체 출력.
pd.set_option('display.max_columns', None)

# 1,000 단위 표기
pd.options.display.float_format = '{:,.3f}'.format

##################### Matplotlib #####################

# Unicode warning 제거 (폰트 관련 경고메시지)
plt.rcParams['axes.unicode_minus']=False

### 한글 폰트 설정 ###
# 나눔고딕 폰트
plt.rcParams['font.family'] = "NanumGothic"
# 맑은고딕 폰트
# plt.rcParams['font.family'] = "Malgun Gothic"
# 애플고딕 (Mac OS)
# plt.rcParams['font.family'] = "AppleGothic"

# 그래프 출력 사이즈 설정
plt.rcParams["figure.figsize"] = (10, 6)

# Seaborn
# sns.set_style("white")

COLORS = {
    'green': '#278F45',
    'red': '#DB4A25', 
    'blue': '#23D3DB',
    'black': '#000000', 
    'yellow': '#DBBE2E',
}

%matplotlib inline
######################## Others  #########################

# 노트북 출력에 나타나는 경고(warning) 메시지 무시
warnings.filterwarnings('ignore')

데이터 로드

# Finance DataReader 설치
# !pip install -U finance-datareader

import FinanceDataReader as fdr

df = fdr.DataReader('BTC/USD')
df.head()
Open High Low Close Adj Close Volume
Date
2014-09-17 465.864 468.174 452.422 457.334 457.334 21056800
2014-09-18 456.860 456.860 413.104 424.440 424.440 34483200
2014-09-19 424.103 427.835 384.532 394.796 394.796 37919700
2014-09-20 394.673 423.296 389.883 408.904 408.904 36863600
2014-09-21 408.085 412.426 393.181 398.821 398.821 26580100

freq에 원하는 시계열 데이터의 주기를 입력할 수 있습니다.

  • 참고링크: https://pandas.pydata.org/docs/user_guide/timeseries.html#timeseries-offset-aliases
# 주간(Weekly) 데이터만 추출하기 위하여 date_range()를 주 단위로 생성합니다.
target_index = pd.date_range('20130107', '20230701', freq='W-MON')
target_index

DatetimeIndex(['2013-01-07', '2013-01-14', '2013-01-21', '2013-01-28',
               '2013-02-04', '2013-02-11', '2013-02-18', '2013-02-25',
               '2013-03-04', '2013-03-11',
               ...
               '2023-04-24', '2023-05-01', '2023-05-08', '2023-05-15',
               '2023-05-22', '2023-05-29', '2023-06-05', '2023-06-12',
               '2023-06-19', '2023-06-26'],
              dtype='datetime64[ns]', length=547, freq='W-MON')

# 추출한 시계열 인덱스 데이터로 필터링
df.loc[df.index.isin(target_index)].head(10)
Open High Low Close Adj Close Volume
Date
2014-09-22 399.100 406.916 397.130 402.152 402.152 24127600
2014-09-29 376.928 385.211 372.240 375.467 375.467 32497700
2014-10-06 320.389 345.134 302.560 330.079 330.079 79011800
2014-10-13 377.921 397.226 368.897 390.414 390.414 35221400
2014-10-20 389.231 390.084 378.252 382.845 382.845 16419000
2014-10-27 354.777 358.632 349.809 352.989 352.989 13033000
2014-11-03 325.569 334.002 325.481 327.554 327.554 12948500
2014-11-10 362.265 374.816 357.561 366.924 366.924 30450100
2014-11-17 388.349 410.199 377.502 387.408 387.408 41518800
2014-11-24 366.948 387.209 366.669 376.901 376.901 30930100 
# 2020년 ~ 2022년 3년치의 시세 데이터를 조회합니다. (Close)는 종가를 의미
data = df.loc[(df.index >= '2020-01-01') & (df.index <= '2022-12-31'), 'Close']
data

Date
2020-01-01    7,200.174
2020-01-02    6,985.470
2020-01-03    7,344.884
2020-01-04    7,410.657
2020-01-05    7,411.317
                ...    
2022-12-27   16,717.174
2022-12-28   16,552.572
2022-12-29   16,642.342
2022-12-30   16,602.586
2022-12-31   16,547.496
Name: Close, Length: 1096, dtype: float64

시계열 데이터 결측치 처리

Simple Imputer

  • 단순하게 평균, 중앙값 등으로 채웁니다.

  • 하지만, 시계열 데이터의 특성상 앞,뒤 데이터의 추세를 반영하지 못합니다.

# 결측치 확인
data.isnull().sum()

0

# 평균 값으로 채우는 경우
data.fillna(data.mean())

Date
2020-01-01    7,200.174
2020-01-02    6,985.470
2020-01-03    7,344.884
2020-01-04    7,410.657
2020-01-05    7,411.317
                ...    
2022-12-27   16,717.174
2022-12-28   16,552.572
2022-12-29   16,642.342
2022-12-30   16,602.586
2022-12-31   16,547.496
Name: Close, Length: 1096, dtype: float64

Interpolate

  • 앞, 뒤 데이터의 추세를 보고 interpolation 값을 계산하여 채웁니다.
# interpolation 방식으로 결측치를 채웁니다.
data = data.interpolate()
data

Date
2020-01-01    7,200.174
2020-01-02    6,985.470
2020-01-03    7,344.884
2020-01-04    7,410.657
2020-01-05    7,411.317
                ...    
2022-12-27   16,717.174
2022-12-28   16,552.572
2022-12-29   16,642.342
2022-12-30   16,602.586
2022-12-31   16,547.496
Name: Close, Length: 1096, dtype: float64

# 결측치 확인
data.isnull().sum()

0

지수평활 (Exponential Smoothing)

지수평활(Exponential Smoothing) 은 시계열 데이터 분석에서 일반적으로 사용되는 기법 중 하나입니다. 이는 데이터의 패턴을 예측하는 데 도움이 되는, 데이터의 노이즈를 줄이고 추세를 더 잘 파악 할 수 있는 방법입니다.


지수평활의 기본적인 아이디어는 최근의 관찰값에 더 많은 가중치 를 주고, 과거의 관찰값에는 적은 가중치 를 주는 것입니다. ‘지수’라는 용어는 각 관찰값에 곱해지는 가중치가 지수적으로 감소하기 때문에 사용됩니다.


이 방법은 다음과 같은 이유로 사용됩니다:

  • 간단함: 지수평활은 계산이 간단하고 이해하기 쉬운 방법입니다. 이는 데이터가 복잡한 패턴이나 추세를 보이지 않을 때 특히 유용합니다.

  • 노이즈 감소: 지수평활은 원래의 시계열 데이터에서 노이즈를 줄이고 주요한 패턴을 강조함으로써 데이터를 더 잘 이해하는 데 도움이 됩니다.

  • 예측: 지수평활은 미래 값의 예측에 사용할 수 있습니다. 이 방법은 과거의 데이터를 사용하여 미래의 데이터를 예측하는 방법 중 하나로, 데이터의 패턴이 시간이 지남에 따라 일정하게 유지될 것으로 예상되는 경우 특히 유용합니다.

  • 빠른 업데이트: 지수평활은 새로운 데이터가 추가될 때마다 쉽게 업데이트될 수 있습니다. 이는 특히 실시간 분석이 필요한 경우에 유용합니다.


지수평활에는 다양한 변형이 있습니다.

  • 가장 간단한 형태의 지수평활은 단순 지수평활(Simple Exponential Smoothing, SES)이며, 이는 데이터에 선형적인 패턴이 없을 때 사용됩니다.
  • 더 복잡한 형태의 지수평활은 홀트의 선형 지수평활(Holt's Linear Exponential Smoothing) 이나 홀트-윈터스의 계절적 지수평활(Holt-Winters' Seasonal Exponential Smoothing) 과 같이 트렌드와 계절성을 포함하는 데이터를 다룰 수 있습니다.
fig, ax = plt.subplots(1, 1)
fig.set_size_inches(10, 6)
fig.set_dpi(300)
ax.plot(data, label='Close', color='black', linestyle='solid', linewidth=0.5)
ax.plot(data.ewm(alpha=0.1).mean(), label='alpha=0.1', color=COLORS['red'], linestyle='solid', linewidth=0.8)
ax.plot(data.ewm(alpha=0.5).mean(), label='alpha=0.5', color=COLORS['blue'], linestyle='solid', linewidth=0.7)
ax.plot(data.ewm(alpha=0.8).mean(), label='alpha=0.8', color=COLORS['green'], linestyle='solid', linewidth=0.5)

plt.xticks(rotation=30)
plt.legend()
plt.show()

이동평균(Moving Average)

이동 평균은 주어진 시점에서 이전 n 개의 데이터 포인트의 평균을 계산 하는 것을 의미합니다.

이는 시계열 데이터의 잡음을 줄이고 데이터의 기본적인 패턴을 확인하는 데 도움이 됩니다.


이동평균(moving average) 를 구하는 방식은 다음과 같습니다.


(예시)

  • 10일 이동평균을 계산

data.rolling(window=10).mean()

fig, ax = plt.subplots(1, 1)
fig.set_size_inches(10, 6)
fig.set_dpi(300)
ax.plot(data, label='Close', color='black', linestyle='solid', linewidth=0.5)
ax.plot(data.rolling(window=10).mean().bfill(), label='ma10', color=COLORS['red'], linestyle='solid', linewidth=0.8)
ax.plot(data.rolling(window=20).mean().bfill(), label='ma20', color=COLORS['blue'], linestyle='solid', linewidth=0.7)
ax.plot(data.rolling(window=60).mean().bfill(), label='ma60', color=COLORS['green'], linestyle='solid', linewidth=0.5)

plt.xticks(rotation=30)
plt.legend()
plt.show()

정상성(stationary) / 비정상성(non-stationary)

이미지 출처: https://t1.daumcdn.net/cfile/tistory/99E7CC3E5D13859D30

정상 시계열 (stationary)

정상 시계열은 시간의 흐름에 따라 그 통계적 특성이 변하지 않는 시계열 을 의미합니다.

다시 말해, 정상 시계열의 평균과 분산(또는 표준편차), 그리고 공분산은 시간이 변해도 일정 합니다.


이러한 특성 때문에 정상 시계열은 분석이 용이하며, 많은 통계적 시계열 모델들(AR, MA, ARMA 등)은 데이터가 정상성을 가정 하고 있습니다. (여기서, 정상성을 가정한다는 점이 매우 중요합니다!!)


(예시)

  • 시간에 따른 일정한 주기를 가진 신호(일정한 주기로 진동하는 심장 박동)

  • 백색 잡음(랜덤한 값들의 연속) 등이 있습니다.

비정상 시계열 (non-stationary)

비정상 시계열은 시간의 흐름에 따라 그 통계적 특성이 변하는 시계열 을 의미합니다. 즉, 평균, 분산, 공분산 등이 시간에 따라 변화합니다.


비정상 시계열 데이터는 시간의 흐름에 따라 통계적 속성이 변하는 데이터를 의미합니다. 이는 추세(trend) 또는 계절성(seasonality) 때문에 발생할 수 있습니다.

  • 추세(trend): 시계열 데이터가 시간에 따라 증가하거나 감소하는 경향 을 보일 때, 이를 추세 라고 합니다.

  • 계절성(seasonality): 시계열 데이터가 일정한 주기로 반복되는 패턴 을 보일 때, 이를 계절성 이라고 합니다.

대부분의 실제 시계열 데이터는 비정상적인 성질을 가지고 있습니다. 이 경우에는 데이터를 정상적으로 만드는 변환(예: 차분, 로그 변환 등) 이 필요할 수 있습니다.


(예시)

  • 시간에 따라 트렌드가 있는 주가 데이터

  • 인구 증가 데이터 등

이러한 데이터는 시간에 따라 평균과 분산이 변하기 때문에 비정상적(non-stationary) 입니다.

정상성(stationary) / 비정상성(non-stationary) 확인

시각화를 통한 확인

ACF(Auto Correlation Function) Plot - 시각화

  • 링크: https://www.statsmodels.org/stable/generated/statsmodels.graphics.tsaplots.plot_acf.html

분석 tip 1. 통계적 유의성 확인

  • 자기 상관이 통계적으로 유의미한지 여부를 나타내는 “신뢰구간”을 회색 음영 영역(shading area) 으로 표시

  • 기본적으로 95% 신뢰구간 을 사용하여 shading area를 그립니다.

  • 각 Lag 에서 Shading 영역을 벗어나 자기상관계수(Autocorrelation) 추정치는 95% 신뢰구간을 벗어났기 때문에, 통계적 유의성을 찾을 수 없다 는 의미로 해석할 수 있습니다.

분석 tip 2. 정상성 확인

  • 보통은 시차가 커질수록 상관계수는 점차 0에 가까워지며, 시차가 커질 수록 자기상관성이 떨어지는 양상을 보여주는 것이 일반적입니다. 하지만, Seasonal 한 데이터의 경우, 값이 다시 튀어오르는 현상은 있을 수 있습니다.

  • ACF는 시계열 데이터의 정상성을 판단할 때 유용 하게 활용될 수 있습니다.

    • 빠르게 0으로 수렴하는 경우 -> 정상(stationary) 시계열, 천천히 감소합니다.

    • 큰 상관계수를 가지는 경우 비정상(non-stationary) 시계열 데이터로 평가할 수 있습니다.

수식

  • ACF는 시계열 데이터의 lagged 복사본과 원본 사이의 상관 관계를 측정합니다.


$\rho_k = \frac{Cov(X_t, X_{t+k})}{\sqrt{Var[X_t]*Var[X_{t+k}]}}$


PACF(Partial Autocorrelation Function) Plot - 시각화

  • 시간차(lag)가 늘어나면서 관측치 간의 상관 관계가 어떻게 변하는지 설명 합니다.

  • 원리적으로, PACF는 선형 회귀분석을 통해 계산됩니다. 각 시차의 관측치를 이전 시차들의 관측치로 선형 회귀분석 하고, 이때의 회귀 계수를 PACF로 간주합니다. 이렇게 하면, PACF는 이전 시차들의 영향을 배제하고 특정 시차만의 자기상관을 측정할 수 있습니다.

분석 tip 1. AR 모델의 차수 결정

  • PACF를 사용하면 AR 모델의 차수 를 결정하는 데 도움이 될 수 있습니다.

수식

  • PACF는 다른 모든 시점의 영향을 배제하고 두 시점 사이의 상관 관계를 측정합니다.


$\phi_{kk} = Cor r(X_t, X_{t+k} | X_{t+1}, …, X_{t+k-1})$


from statsmodels.graphics.tsaplots import plot_acf, plot_pacf


fig, axes = plt.subplots(1, 2)
fig.set_size_inches(12, 4)

# ACF Plot
plot_acf(data, lags=30, ax=axes[0])
# PACF Plot
plot_pacf(data, lags=10, zero=False, ax=axes[1])

for ax in axes:
    ax.set_ylim(-0.5, 1.25)
plt.show()

정상성 검증

kpss test

  • 귀무가설: 해당 시계열은 정상(stationary) 시계열이다. <-> 대립가설: 해당 시계열은 비정상(non-stationary) 시계열이다.

    • p-value <= 0.05 : 귀무가설 기각/ 대립가설 채택 -> 비정상(non-stationary) 시계열.

    • p-value > 0.05 : 귀무가설 채택/ 대립가설 기각 -> 정상(stationary) 시계열.

from statsmodels.tsa.stattools import kpss

def kpss_test(series, **kw):    
    stats, p_value, nlags, critical_values = kpss(series, **kw)
    
    print(f'KPSS Stat: {stats:.5f}')
    print(f'p-value: {p_value:.2f}')
    print(f'Lags: {nlags}')
    
    print(f'검증결과: {"비정상(non-stationary)" if p_value <= 0.05 else "정상(stationary)"} 시계열 데이터입니다.')

kpss_test(data, nlags=30)

KPSS Stat: 1.18454
p-value: 0.01
Lags: 30
검증결과: 비정상(non-stationary) 시계열 데이터입니다.

# 1차 차분에 대한 결과
diff_1 = data.diff(periods=1).iloc[1:]

kpss_test(diff_1, nlags=30)

KPSS Stat: 0.29060
p-value: 0.10
Lags: 30
검증결과: 정상(stationary) 시계열 데이터입니다.

# 1차 차분 결과에 대한 시각화
fig, axes = plt.subplots(1, 1)
fig.set_size_inches(8, 4)

# ACF Plot
plot_acf(diff_1, lags=30, ax=axes)

axes.set_ylim(-0.5, 1.25)
plt.show()

Adfuller 테스트

  • kpss test 와 귀무가설이 반대입니다.

  • 귀무가설: 해당 시계열은 비정상(non-stationary) 시계열이다. <-> 대립가설: 해당 시계열은 정상(stationary) 시계열이다.

    • p-value <= 0.05 : 귀무가설 기각/ 대립가설 채택 -> 정상(stationary) 시계열.

    • p-value > 0.05 : 귀무가설 채택/ 대립가설 기각 -> 비정상(non-stationary) 시계열.

# 데이터 정상 상태 수치 확인 
from statsmodels.tsa.stattools import adfuller #ADF Test를 위한 함수 호출 

st_result = adfuller(data)

def adfuller_test(series, **kw):    
    adf, p_value, nlags, number_of_observations, critical_values, _ = adfuller(series, **kw)
    
    print(f'ADF: {adf:.5f}')
    print(f'p-value: {p_value:.2f}')
    print(f'Lags: {nlags}')
    print(f'Number of Observations: {number_of_observations}')
    
    print(f'검증결과: {"비정상(non-stationary)" if p_value > 0.05 else "정상(stationary)"} 시계열 데이터입니다.')

# Adfuller 테스트
adfuller_test(data, maxlag=30)

ADF: -1.44857
p-value: 0.56
Lags: 0
Number of Observations: 1095
검증결과: 비정상(non-stationary) 시계열 데이터입니다.

# 1차 차분에 대한 Adfuller 테스트
adfuller_test(diff_1, maxlag=30)

ADF: -10.34401
p-value: 0.00
Lags: 8
Number of Observations: 1086
검증결과: 정상(stationary) 시계열 데이터입니다.

ARIMA Model

시계열 데이터 분석에서는 AR(AutoRegressive), MA(Moving Average), ARMA(AutoRegressive Moving Average), ARIMA(AutoRegressive Integrated Moving Average) 등의 모델이 있습니다.

  • AR (AutoRegressive): 이전 관측치의 값들이 현재 관측치에 영향을 주는 모델입니다. 예를 들어, AR(1) 모델은 현재 관측치가 바로 이전 관측치에 의해 영향을 받는다는 것을 나타냅니다. AR 모델은 파라미터 p로 몇 개의 이전 관측치를 고려할지를 결정합니다.

  • MA (Moving Average): 이전 오차항들이 현재 관측치에 영향을 주는 모델입니다. MA(1) 모델은 현재 관측치가 바로 이전의 오차항에 의해 영향을 받는다는 것을 나타냅니다. MA 모델은 파라미터 q로 몇 개의 이전 오차항을 고려할지를 결정합니다.

AR & MA 모델의 차이

AR(AutoRegressive) 모델과 MA(Moving Average) 모델은 둘 다 시계열 분석에서 널리 사용되는 모델이지만, 그 방식이 다릅니다.


AR 모델

  • AR 모델은 ‘자기회귀 모델’ 로, 과거의 자신의 값에 의존하여 현재 값을 예측하는 방법을 사용합니다. 즉, 이전 시점의 관측값이 현재 시점의 관측값에 영향을 줍니다.
  • 예를 들어, 주식 시장에서 어제 주가가 상승했다면 오늘 주가도 상승할 가능성이 있다는 것 이 AR 모델의 관점입니다.


AR 모델은 파라미터 p로 이전 몇 개의 관측치를 고려할지 결정 하게 됩니다. AR(1) 모델은 바로 이전의 한 단계 관측치를 고려하고, AR(2) 모델은 바로 이전의 두 단계 관측치를 고려하게 됩니다.


(예시)

예를 들어, 전력 소비를 예측하는 문제를 생각해 보겠습니다. 사람들은 특정 시간대에 일관되게 전력을 사용하는 경향이 있습니다. 아침에 사람들이 일어나서 활동을 시작하면 전력 사용량이 증가하고, 밤에 사람들이 잠에 들면 전력 사용량이 감소합니다. 이러한 경우, AR 모델은 이전 시간대의 전력 사용 패턴을 학습하여 다음 시간대의 전력 사용량을 예측하는 데 사용 될 수 있습니다.


MA 모델

  • MA 모델은 ‘이동 평균 모델’ 로, 과거의 오차항에 의존하여 현재 값을 예측합니다.
  • 예를 들어, 우리가 날씨를 예측한다고 가정해 봅시다. 오늘의 기온 예측치가 실제 기온보다 더 높았다면, 이는 오차를 의미합니다. MA 모델에 따르면, 이 오차는 내일의 기온 예측에 반영되어야 합니다. 즉, 오늘 예측이 너무 높았다면, 내일의 예측치는 그 오차를 보정하기 위해 조금 더 낮아질 것입니다.


MA 모델은 파라미터 q로 이전 몇 개의 오차항을 고려할지 결정하게 됩니다. MA(1) 모델은 바로 이전의 한 단계 오차항 을 고려하고, MA(2) 모델은 바로 이전의 두 단계 오차항 을 고려하게 됩니다.

즉, AR 모델은 자신의 과거 값을, MA 모델은 과거의 예측 오차를 사용하여 현재 값을 예측하는 방법론입니다. 어떤 모델을 사용할지는 주로 데이터의 특성과 목표에 따라 결정됩니다.


(예시)

매장의 판매량 예측을 예로 들어 보겠습니다. 이전의 판매 예측이 과소평가되었다면 (즉, 실제 판매량이 예측보다 높았다면), 그 다음날의 판매량 예측을 증가시키는 것이 합리적일 수 있습니다. 반대로, 이전의 판매 예측이 과대평가 되었다면 (즉, 실제 판매량이 예측보다 낮았다면), 그 다음날의 판매량 예측을 감소시키는 것이 합리적 일 수 있습니다. 이런 상황에서는 MA 모델이 유용하게 사용될 수 있습니다.

ARMA, ARIMA

  • ARMA (AutoRegressive Moving Average): AR 모델과 MA 모델의 결합으로, 이전 관측치와 이전 오차항 모두를 고려하는 모델입니다. ARMA 모델은 두 파라미터 p와 q를 모두 사용합니다.

  • ARIMA (AutoRegressive Integrated Moving Average): ARMA 모델에 비정상 시계열을 정상 시계열로 변환하는 과정을 포함한 모델입니다. 이 변환 과정은 차분(differencing)이라고 하며, 몇 차 차분을 사용할지를 결정하는 파라미터가 d입니다. 따라서 ARIMA 모델은 세 파라미터 p, d, q를 모두 사용합니다.

ARIMA 적용시 고려사항

  • AR(p) = ARIMA(p, 0, 0)

  • MA(q) = ARIMA(0, 0, q)

  • ARMA(p, q) = ARIMA(p, 0, q)

from statsmodels.tsa.arima.model import ARIMA
from statsmodels.tsa.statespace.sarimax import SARIMAX
import itertools
from tqdm import tqdm
import warnings

warnings.filterwarnings('ignore')

p = range(0, 4)
d = range(0, 3)
q = range(0, 5)

pdq = list(itertools.product(p,d,q))

aic = []
params = []

with tqdm(total=len(pdq)) as pg:
    for i in pdq:
        pg.update(1)
        try:
            model = ARIMA(data, order=(i)).fit()
            aic.append(round(model.aic, 2))
            params.append((i))
        except:
            continue 

100%|███████████████████████████████████████████| 60/60 [00:22<00:00,  2.66it/s]

가장 Optimal 한 (p, d, q) 하이퍼 파라미터를 찾습니다.

optimal = [(params[i],j) for i,j in enumerate(aic) if j == min(aic)]
arima = ARIMA(data, order = optimal[0][0], freq='D').fit()
arima.summary()
SARIMAX Results
Dep. Variable: Close No. Observations: 1096
Model: ARIMA(2, 2, 3) Log Likelihood -9369.370
Date: Mon, 03 Jul 2023 AIC 18750.740
Time: 01:46:18 BIC 18780.726
Sample: 01-01-2020 HQIC 18762.087
- 12-31-2022
Covariance Type: opg
coef std err z P>|z| [0.025 0.975]
ar.L1 -1.9073 0.014 -137.751 0.000 -1.934 -1.880
ar.L2 -0.9732 0.014 -71.464 0.000 -1.000 -0.946
ma.L1 0.8959 0.030 29.584 0.000 0.837 0.955
ma.L2 -0.9402 0.046 -20.326 0.000 -1.031 -0.850
ma.L3 -0.9557 0.028 -34.322 0.000 -1.010 -0.901
sigma2 1.637e+06 8.1e-08 2.02e+13 0.000 1.64e+06 1.64e+06
Ljung-Box (L1) (Q): 0.03 Jarque-Bera (JB): 1219.10
Prob(Q): 0.86 Prob(JB): 0.00
Heteroskedasticity (H): 5.39 Skew: -0.19
Prob(H) (two-sided): 0.00 Kurtosis: 8.16 


arima_model = ARIMA(data, order=optimal[0][0], freq='D').fit()
forecast = arima_model.forecast(steps=10)

fig, ax = plt.subplots(1, 1)
fig.set_size_inches(10, 5)
fig.set_dpi(300)
ax.plot(data.index, data, color=COLORS['green'])
ax.plot(forecast, color=COLORS['red'])
plt.show()

결과를 시각화 합니다. 다음의 시각화 결과는 향후 15일 치 데이터에 대한 예측 결과입니다.

fig, axes = plt.subplots(1, 2)
fig.set_size_inches(12, 6)
fig.set_dpi(300)
 
# number of prediction days
n_futures = 15

freq = 'D'
future_idx = pd.date_range(data.index[-1], periods=n_futures+1, freq=freq)[-n_futures:]
 
forecast = arima_model.fittedvalues
forecast_futures = arima_model.forecast(steps=n_futures)

axes[0].plot(data.index, data, label='data', color=COLORS['black'], linestyle='dotted', alpha=0.5)
axes[0].plot(data.index, forecast, label='fitted', color=COLORS['green'], linestyle='solid', alpha=0.5)
axes[0].set_title('Fitted Result')
axes[0].set_xticklabels(mdates.num2date(axes[0].get_xticks()), rotation = 45)
axes[0].xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m'))

axes[1].plot(future_idx, forecast_futures, color=COLORS['red'], label='prediction')
axes[1].plot(data.index[-30:], data[-30:], color=COLORS['black'], linestyle='dotted', label='real')
axes[1].set_title(f'Predict Next {n_futures} days')
axes[1].set_xticklabels(mdates.num2date(axes[1].get_xticks()), rotation = 45)
axes[1].xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m-%d'))

# plt.xticks(rotation=45)
plt.legend()
plt.tight_layout()
plt.show()

SARIMAX

SARIMAX: Seasonal AutoRegressive Integrated Moving Average with eXogenous regressors


주요 하이퍼파라미터

  • AR(AutoRegressive): 이전 관측치의 값 을 사용하여 현재 값을 예측하는데 사용되는 모델입니다.

  • I(Integrated): 시계열의 비정상성을 처리하기 위해 사용되는 차분(Differencing) 단계를 나타냅니다.

  • MA(Moving Average): 이전 오차항의 평균 을 사용하여 현재 값을 예측하는데 사용되는 모델입니다.

  • Seasonality: 계절성 요인 을 고려한 모델입니다. 예를 들어, 연간 패턴, 월간 패턴, 주간 패턴 등을 모델링 할 수 있습니다.

  • Exogenous: 시계열에 영향을 주는 외부 변수 를 고려하는 모델입니다.

ARIMA 모델은 AR, I, MA의 세 가지 구성요소만을 포함합니다. 따라서 계절성과 외생변수는 ARIMA 모델로는 처리하기 어렵습니다. 이 때문에 이러한 요소를 고려해야 하는 복잡한 시계열 데이터를 분석할 때는 ARIMA 보다는 SARIMAX와 같은 확장된 모델 을 사용하는 것이 적합할 수 있습니다.


즉, ARIMA와 SARIMAX의 주요 차이점은 SARIMAX가 ARIMA에 비해 계절성 과 외생변수를 모델링할 수 있다는 점입니다. 이로 인해 SARIMAX는 ARIMA보다 더 다양한 시계열 데이터를 분석할 수 있으며, 이에 따라 예측의 정확성을 높일 수 있습니다.


공식문서: https://www.statsmodels.org/dev/generated/statsmodels.tsa.statespace.sarimax.SARIMAX.html

sarimax = SARIMAX(data, order=(2, 2, 3), season_order = (1, 1, 2, 14)).fit(disp=False)
forecast = sarimax.forecast(steps=30)

plt.figure(figsize=(10, 5))
plt.plot(data.index, data, label='real', color=COLORS['green'])
plt.plot(forecast, label='predict', color=COLORS['red'])
plt.legend()
plt.show()

Finance DataReader 라이브러리를 활용하여 금 선물가격 과 코스피 지수 데이터를 가져와서 병합하도록 하겠습니다.

# 금 선물
gold = fdr.DataReader('ZG')
# 금 선물 데이터 컬럼 변환
gold = gold['Close'].reset_index().rename(columns={'Close': 'Gold'})
gold
Date Gold
0 2011-07-20 11.154
1 2011-07-21 10.352
2 2011-07-22 10.686
3 2011-07-25 11.101
4 2011-07-26 11.360
... ... ...
3002 2023-06-26 48.190
3003 2023-06-27 49.980
3004 2023-06-28 51.390
3005 2023-06-29 49.010
3006 2023-06-30 49.200

3007 rows × 2 columns

# 코스피 데이터
kospi = fdr.DataReader('KS11')
# 코스피 데이터 병합
kospi = kospi['Close'].reset_index().rename(columns={'Close': 'Kospi'})
kospi
Date Kospi
0 1996-12-11 704.680
1 1996-12-12 689.380
2 1996-12-13 689.070
3 1996-12-16 673.920
4 1996-12-17 663.350
... ... ...
6689 2023-06-26 2,582.200
6690 2023-06-27 2,581.390
6691 2023-06-28 2,564.190
6692 2023-06-29 2,550.020
6693 2023-06-30 2,564.280

6694 rows × 2 columns

# Kospi 데이터와 금값 시세 데이터 병합
kospi_gold = pd.merge(kospi, gold, how='left')
kospi_gold
Date Kospi Gold
0 1996-12-11 704.680 NaN
1 1996-12-12 689.380 NaN
2 1996-12-13 689.070 NaN
3 1996-12-16 673.920 NaN
4 1996-12-17 663.350 NaN
... ... ... ...
6689 2023-06-26 2,582.200 48.190
6690 2023-06-27 2,581.390 49.980
6691 2023-06-28 2,564.190 51.390
6692 2023-06-29 2,550.020 49.010
6693 2023-06-30 2,564.280 49.200

6694 rows × 3 columns

병합한 데이터에서 결측치가 관측됩니다. 결측치를 interpolate() 으로 결측치를 채우도록 하겠습니다.

kospi_gold['Gold'] = kospi_gold['Gold'].interpolate().bfill()
kospi_gold['Kospi'] = kospi_gold['Kospi'].interpolate().bfill()
kospi_gold
Date Kospi Gold
0 1996-12-11 704.680 11.154
1 1996-12-12 689.380 11.154
2 1996-12-13 689.070 11.154
3 1996-12-16 673.920 11.154
4 1996-12-17 663.350 11.154
... ... ... ...
6689 2023-06-26 2,582.200 48.190
6690 2023-06-27 2,581.390 49.980
6691 2023-06-28 2,564.190 51.390
6692 2023-06-29 2,550.020 49.010
6693 2023-06-30 2,564.280 49.200

6694 rows × 3 columns

data 와 kospi_gold 데이터를 병합합니다.

data_merged = pd.merge(data, kospi_gold, left_on=data.index, right_on=kospi_gold['Date'], how='left').drop('Date', axis=1)
data_merged = data_merged.rename(columns={'key_0': 'Date'})
data_merged
Date Close Kospi Gold
0 2020-01-01 7,200.174 NaN NaN
1 2020-01-02 6,985.470 2,175.170 45.000
2 2020-01-03 7,344.884 2,176.460 44.560
3 2020-01-04 7,410.657 NaN NaN
4 2020-01-05 7,411.317 NaN NaN
... ... ... ... ...
1091 2022-12-27 16,717.174 2,332.790 31.540
1092 2022-12-28 16,552.572 2,280.450 30.390
1093 2022-12-29 16,642.342 2,236.400 31.100
1094 2022-12-30 16,602.586 NaN NaN
1095 2022-12-31 16,547.496 NaN NaN

1096 rows × 4 columns

# Gold, Kospi 의 결측치를 채웁니다.
data_merged['Gold'] = data_merged['Gold'].interpolate().bfill()
data_merged['Kospi'] = data_merged['Kospi'].interpolate().bfill()
data_merged
Date Close Kospi Gold
0 2020-01-01 7,200.174 2,175.170 45.000
1 2020-01-02 6,985.470 2,175.170 45.000
2 2020-01-03 7,344.884 2,176.460 44.560
3 2020-01-04 7,410.657 2,169.330 44.517
4 2020-01-05 7,411.317 2,162.200 44.473
... ... ... ... ...
1091 2022-12-27 16,717.174 2,332.790 31.540
1092 2022-12-28 16,552.572 2,280.450 30.390
1093 2022-12-29 16,642.342 2,236.400 31.100
1094 2022-12-30 16,602.586 2,236.400 31.100
1095 2022-12-31 16,547.496 2,236.400 31.100

1096 rows × 4 columns

SARIMAX 모델로 학습하고, 예측 후 검증결과를 확인하기 위하여 데이터셋을 분할합니다.

  • 30일 기준으로 분할하여, train 데이터로 학습된 모델로 test 데이터를 예측하도록 하겠습니다.
# train / test 분할
train_data = data_merged[:-30]
test_data = data_merged[-30:]

# 향후 30일의 데이터에 대한 예측
n_predictions = 30

sarimax = SARIMAX(train_data['Close'], 
                exog=train_data[['Gold', 'Kospi']], 
                order=(2, 2, 3), 
                season_order = (1, 1, 2, 14)).fit(disp=False)

forecast = sarimax.forecast(steps=30, exog=test_data[['Gold', 'Kospi']])

fig, axes = plt.subplots(1, 2)
fig.set_size_inches(12, 5)
fig.set_dpi(300)

axes[0].plot(train_data['Date'], train_data['Close'], label='real', color=COLORS['green'])
axes[0].plot(test_data['Date'], forecast, label='predict', color=COLORS['red'])
axes[0].set_xlim(pd.Timestamp('2020-09-01'), )
# Tick을 날짜형식으로 변환
axes[0].set_xticklabels(mdates.num2date(axes[0].get_xticks()), rotation = 45)
axes[0].xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m'))
axes[0].set_title('Overall Time Series')
axes[0].legend()

axes[1].plot(train_data['Date'], train_data['Close'], label='real', color=COLORS['green'])
axes[1].plot(test_data['Date'], forecast, label='predict', color=COLORS['red'])
axes[1].set_xlim(pd.Timestamp('2022-11-01'), pd.Timestamp('2023-01-15'))
# Tick을 날짜형식으로 변환
axes[1].set_xticklabels(mdates.num2date(axes[1].get_xticks()), rotation = 45)
axes[1].xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m'))
axes[1].set_title(f'Predict Next {n_predictions} days')
axes[1].legend()

plt.tight_layout()
plt.legend()
plt.show()

SARIMAX 하이퍼파라미터 튜닝

p = range(0, 5)  # AR
d = range(1, 3)  # 차분
q = range(0, 4)  # MA
m = 30           # 계절성 주기(Seasonal Trends)

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

aic = []
params = []

with tqdm(total = len(pdq) * len(seasonal_pdq)) as pg:
    for i in pdq:
        for j in seasonal_pdq:
            pg.update(1)
            try:
                model = SARIMAX(train_data['Close'], 
                                exog=train_data[['Gold', 'Kospi']], 
                                order=(i), 
                                season_order = (j)).fit(disp=False)
                aic.append(round(model.aic, 2))
                params.append((i,j))
            except:
                print('Error Occured..')
                continue

100%|███████████████████████████████████████| 1600/1600 [17:53<00:00,  1.49it/s]

# 향후 30일의 데이터에 대한 예측
n_predictions = 30

optimal = [(params[i], j) for i, j in enumerate(aic) if j == min(aic)]

sarimax_tunned = SARIMAX(train_data['Close'], 
                         exog=train_data[['Gold', 'Kospi']], 
                         order = optimal[0][0][0], 
                         seasonal_order = optimal[0][0][1]).fit(disp=False)

sarimax_tunned.summary()
SARIMAX Results
Dep. Variable: Close No. Observations: 1066
Model: SARIMAX(2, 1, 2)x(0, 1, [], 30) Log Likelihood -9211.826
Date: Mon, 03 Jul 2023 AIC 18437.652
Time: 02:04:26 BIC 18472.247
Sample: 0 HQIC 18450.778
- 1066
Covariance Type: opg
coef std err z P>|z| [0.025 0.975]
Gold 80.1426 12.380 6.474 0.000 55.879 104.406
Kospi 0.4043 1.656 0.244 0.807 -2.842 3.651
ar.L1 -1.8651 0.012 -150.177 0.000 -1.889 -1.841
ar.L2 -0.9719 0.012 -81.006 0.000 -0.995 -0.948
ma.L1 1.8475 0.013 140.016 0.000 1.822 1.873
ma.L2 0.9670 0.013 75.694 0.000 0.942 0.992
sigma2 3.265e+06 4.59e-05 7.11e+10 0.000 3.27e+06 3.27e+06
Ljung-Box (L1) (Q): 1.17 Jarque-Bera (JB): 242.89
Prob(Q): 0.28 Prob(JB): 0.00
Heteroskedasticity (H): 5.96 Skew: 0.02
Prob(H) (two-sided): 0.00 Kurtosis: 5.37 


sarimax_tunned.forecast(steps=n_predictions, exog=test_data[['Gold', 'Kospi']])

1066   16,750.097
1067   16,527.811
1068   17,467.031
1069   17,493.765
1070   17,086.020
1071   16,559.103
1072   14,505.218
1073   11,872.131
1074   13,345.895
1075   12,583.305
1076   12,389.935
1077   12,193.529
1078   12,474.401
1079   12,560.788
1080   12,330.134
1081   12,341.078
1082   12,425.906
1083   12,470.121
1084   11,983.324
1085   11,553.372
1086   11,896.263
1087   12,174.399
1088   12,190.526
1089   12,071.561
1090   12,055.353
1091   11,979.975
1092   11,682.838
1093   11,908.014
1094   12,426.192
1095   12,175.555
Name: predicted_mean, dtype: float64

예측 결과 시각화

forecast = sarimax_tunned.forecast(steps=n_predictions, exog=test_data[['Gold', 'Kospi']])

plt.figure(figsize=(10, 5))
plt.plot(data.index, data, label='real', color='#1AA341', alpha=0.5)
plt.plot(test_data['Date'], forecast, label='predict', color='#F02828')
plt.axvline(x = pd.Timestamp('2022-12-02'), color = 'blue', label = 'After Prediction', linestyle='dotted')
plt.xlim(pd.Timestamp('2022-11-01'), )
plt.ylim(10000, 25000)
plt.xticks(rotation=45)
plt.legend()
plt.show()

시계열 성분 분해 (Seasonal Decomposition)

시계열 데이터 분석에서 성분 분해(Decomposition)는 관측된 시계열 데이터를 여러 성분으로 나누는 것을 의미합니다. 이는 시계열 데이터에 있는 복잡한 패턴을 이해하고 분석하는 데 도움이 됩니다. 주로 데이터를 추세(Trend), 계절성(Seasonality), 그리고 잔차(Redisual)/불규칙성(Irregularity) 으로 분해합니다.

  • 추세 성분(Trend Component) : 이는 시계열 데이터의 장기적인 패턴 을 의미합니다. 예를 들어, 회사의 매출이 시간이 지남에 따라 전반적으로 증가하거나 감소하는 경향이 있다면, 이를 추세라고 합니다. 추세는 일반적으로 선형적일 수도 있고, 비선형적일 수도 있습니다.

  • 계절 성분(Seasonal Component) : 이는 시계열 데이터에서 반복적으로 나타나는 패턴 을 의미합니다. 예를 들어, 매년 여름에 아이스크림 판매가 증가하는 것이 계절성의 예입니다. 계절성은 일반적으로 연간, 분기별, 월별, 주별, 일별 등 주기적으로 나타납니다.

  • 불규칙 성분(Irregular Component) : 이는 추세나 계절성에 의해 설명되지 않는 잔차(Residual) 를 의미합니다. 불규칙성은 랜덤 노이즈일 수도 있고, 데이터에 있는 예측 불가능한 변동성일 수도 있습니다. 예를 들어, 예측할 수 없는 이벤트(자연재해, 금융위기 등)에 의한 판매량의 변동이 이에 해당됩니다.

additive 방식을 활용한 분해


$Y_t = T_t + S_t + I_t$


  • 시계열 데이터가 선형적인 형태 라는 가정하에 데이터를 분해합니다.

  • 즉, 데이터에서의 추세, 계절성, 그리고 잔차가 더해져 원래의 시계열 데이터를 구성한다는 가정입니다.

이는 시계열 데이터의 계절적 효과가 시간이 흐름에 따라 일정하게 유지되는 경우에 적합합니다.

import statsmodels.api as sm

# additive 방식을 활용한 분해
decomposition = sm.tsa.seasonal_decompose(data, model='additive', period=365)
fig = decomposition.plot()
fig.set_dpi(300)
fig.set_tight_layout('constrained')
fig.set_size_inches(10, 6)
plt.show()

Multiplicative Model


$Y_t = T_t \times S_t \times I_t$


이 모델은 세 가지 성분이 곱해져서 시계열 데이터를 구성한다고 가정합니다

  • 시계열 데이터가 비선형적인 형태, 특히 시간이 지남에 따라 계절적 효과가 증가하거나 감소하는 경우에 적합합니다.

  • 즉, 데이터에서의 추세, 계절성, 그리고 잔차가 곱해져 원래의 시계열 데이터를 구성한다는 가정입니다.

(주의) multiplicative 모델을 사용하기 위해서는 데이터에 0이 존재해서는 안됩니다.

decomposition = sm.tsa.seasonal_decompose(data, model='multiplicative', period=365)
fig = decomposition.plot()
fig.set_dpi(300)
fig.set_tight_layout('constrained')
fig.set_size_inches(10, 6)
plt.show()

additive 와 multiplicative의 차이

additive 모델은 시간에 따라 빈도와 진폭이 변하지 않는다는 점에 유의하세요.

multiplicative 모델은 이 두 번째 모델에서 행동이 증가하는 깔때기 역할을 합니다(감소할 수도 있음).

  • additive 는 Trend 와 Seasonality 를 별개로 보고,

  • multiplicative 는 Trend 에 따라 Seasonality 도 변화한다고 보면 됩니다.

따라서 데이터의 계절적 패턴의 크기가 데이터의 크기에 따라 달라지는 경우 multiplicative 모델을 사용합니다. 반면, additive 모델에서는 계절성의 크기가 시간에 따라 변하지 않습니다.

출처: https://sigmundojr.medium.com/seasonality-in-python-additive-or-multiplicative-model-d4b9cf1f48a7

댓글남기기

참고

poetry 의 거의 모든것 (튜토리얼)

5 분 소요

Python 개발에 있어서 poetry는 매우 강력한 도구로, 프로젝트의 의존성 관리와 패키지 배포를 간소화하는 데 큰 도움을 줍니다. 지금부터 poetry 활용 튜토리얼을 살펴 보겠습니다.

[Assistants API] Code Interpreter, Retrieval, Functions 활용법

35 분 소요

OpenAI의 새로운 Assistants API는 대화와 더불어 강력한 도구 접근성을 제공합니다. 본 튜토리얼은 OpenAI Assistants API를 활용하는 내용을 다룹니다. 특히, Assistant API 가 제공하는 도구인 Code Interpreter, Retrieval...