베이지안을 활용한 비즈니스 시계열 예측

Prophet Model

Bayesian Concept을 활용해 비즈니스 시계열 추정을 하는 가장 대표적인 모형이 페이스북에서 릴리즈한 Prophet Model이다.

https://peerj.com/preprints/3190/

$$y(t) = g(t)+s(t)+h(t)+\varepsilon_t$$

 

- g(t): 추세 함수(비주기적 변화 포착)

- s(t): 주기적 변화(예) 주간 혹은 연간 계절성)

- h(t): 비규칙적인 이벤트

 

세 가지의 성분을 덧셈으로 연결하였다.(GAM: 일반화가법모형) GAM은 반응평균의 함수를 독립변수들의 함수꼴로 모형화시킨것을 말한다. 컨셉은 결국 시계열 그래프를 가장 그럴듯한 함수 개형을 찾아서 피팅시키는 건데, 그 함수 개형을 크게 3개로 나눠 전반적인 트렌드, 반복하는 시즈널리티, 불규칙적인 이벤트 이렇게 세가지 함수를 각각 찾아서 더해주자는 것이다.

$$g(E(Y|X)) = \beta_0 + f_1(X_1) + f_2(X_2) + \cdots + f_p(X_p)$$

 

https://github.com/facebook/prophet/tree/main/python/prophet

(implementation code)

 

먼저 함수에서 굵직한 트렌드를 어떻게 베이지안으로 잘 찾아낼 수 있는건지, 다음을 보자.

Trend Model (logistic growth model)

$g(t) = \dfrac{C}{1+exp(-k(t-m))}$

g(t) 함수 form을 이와 같이 가정하자. $k$는 성장률(growth rate), $m$은 위치조정 파라미터(offset parameter), $C$는 자연 성장 시 최대값(carrying capacity)을 의미한다. 자연 성장세를 가정하여 로지스틱 함수 꼴을 이용하였다.

 

기본적인 로지스틱 함수 1 / (1+e^(-x)) 의 개형이다

 

기본적인 개형에서 1), 2)의 변화를 줘서 조금 더 유연한 모형을 만들어보자.

 

1) C의 경우 시간이 변함에 따라 곱해주는 상수 값이다. 전체적인 타겟 볼륨이 시간이 흐름에 따라 커질 수 있으므로(예시) 타겟이 금액인 경우, 시간이 흐름에 따라 사람들의 전체적인 소비액이 증가한다고 가정) 상수 대신 변수 $C(t)$로 대체한다.

 

2) k의 경우 값이 커질수록 기울기가 가파르게 변한다. 변화율을 조절하는 파라미터인 k도 시간에 따라 바뀔 수 있다. 다시 말해서 시점에 따라서 성장세(그래프 기울기)를 조절해주는 방식이다. 따라서 k대신 $k+a(t)^T\mathbf{\delta}$가 들어가야한다. 

 

시점에 따라 성장세를 조절해주는 방식이 중요한 컨셉인 것 같은데, 그러한 시점들은 어떻게 정할 수 있을까?

 

Trend changepoints

$$a_j(t) =  \begin{cases} 1, \quad if \; t\geq s_j \\ 0, \quad o.w. \end{cases}$$

 

$\mathbf{\delta}$는 벡터인데 각 성분은 $s_j$시간에 벌어진 사건의 변화율(change in rate)을 의미한다. 즉 $a(t)^T \mathbf{\delta}$는 현재 시점(t)기준으로 이전 시점들($s_j$)에 벌어진 사건의 영향력(비율에 미치는 영향 $\delta$)의 합을 의미한다.

 

$\delta$는 sparse prior로 $\delta_j \sim Laplace(0,\tau)$ 를 준다. (라플라스 분포는 정규 분포가 더 뾰족해진 형태이다.)

이 때 $s_j$는 성장률을 변화시키는 사건이나 여러 이벤트들을 반영해주는 식이다. 디폴트 세팅으로는 라플라스 분포를 사전 분포로 부여하고 hierarchical bayesian 방식으로 학습한다.

 

위의 사항을 반영한 최종 모형은 다음과 같다.

$$g(t) \dfrac{C(t)}{1+exp(-(k+a(t)^T\delta)(t-(m+a(t)^T\gamma))}$$

 

그리고 모형에서는 오버피팅을 방지하기 위해 디폴트로 전체 시계열 자료의 앞 80%에서 변화점을 찾는다. 그리고 디폴트로는 최대 25개의 점을 찾는다. 이 점들은 결국 함수의 기울기가 급격하게 변화하는 점을 의미하는 것이다. 

 

Seasonality

주기성 반영을 위해 아래 푸리에 급수를 이용한다. 

$$s(t) = \sum^N_{n=1} \Big( a_n cos \Big( \dfrac{2\pi nt}{P} \Big) + b_n sin\Big( \dfrac{2\pi nt}{P} \Big) \Big)=X(t)\beta$$

P는 우리가 반영을 원하는 주기이다(7이면 주별 효과, 365.25이면 연도별효과)

$\beta = [a_1,b_1,\cdots,a_N,b_N]^T$ 

 

예를 들어, N=10이고 연도 효과를 반영하는 경우는 다음과 같다.

$$s(t) = \Big[ cos\Big( \dfrac{2\pi (1)t }{365.25} , \cdots, sin \Big( \dfrac{2\pi (10)t}{365.25}\Big) \Big]$$

이 때 계수는 generative model로서 $\beta \sim N(0, \sigma^2)$를 취한다.

 

Holidays and Events

$$Z(t) = [ 1(t\in D_1) , ..., 1(t\in D_L) ] \\ h(t)=Z(t)\kappa, \quad \kappa \sim N(0,\nu ^2)$$

$D_i$는 과거 이벤트들 목록이다. 각 이벤트는 $\kappa_i \sim N(0,\nu_i ^2)$를 가진다. 해당 파트가 prophet의 큰 장점 중 하나다. prophet은 결국 GAM 모형이기에 이런 공휴일 효과를 회귀로 반영해 줄 수 있다. 그 반영 효과는 $\nu$에 의해 좌우되는데 이는 각 이벤트가 가지고 있을거라고 예상되는 효과(정규분포의 분산)이다. 

 

깃헙 코드구조는 다음과 같다.

construct_holiday_dataframe(self,dates): make_holidays_df

make_holiday_features(self, dates, holidays): holidays는 construct_holiday_dataframe return all_holidays

 

make_holiday_features를 돌린 결과 서로 다른 연도의 공휴일이더라도 이름이 같으면 같은 공휴일로 인식되고 window에 따라 다르게 반영됨을 확인할 수 있다.(추후 확인 한 번 더 필요)

 

Implementation

https://facebook.github.io/prophet/docs/quick_start.html#python-api

실제 파이썬에서의 동작은 다음과 같이 이루어진다.

Pythono 3.7.9 

 

# Setting
import pandas as pd
from prophet import Prophet

# 외부 url로 데이터 불러올 때 에러 방지
import ssl
ssl._create_default_https_context = ssl._create_unverified_context

# 데이터 수입
df = pd.read_csv('https://raw.githubusercontent.com/facebook/prophet/main/examples/example_wp_log_peyton_manning.csv')

# 모델링
m = Prophet()
m.fit(df)

# periods 기준으로 기존 데이터 프레임에 날짜를 더하여 날짜 테이블을 만듬
future = m.make_future_dataframe(periods=365)
future.tail()

forecast = m.predict(future)
'''
forecast의 열은 다음으로 이루어짐
Index(['ds', 'trend', 'yhat_lower', 'yhat_upper', 'trend_lower', 'trend_upper',
       'additive_terms', 'additive_terms_lower', 'additive_terms_upper',
       'weekly', 'weekly_lower', 'weekly_upper', 'yearly', 'yearly_lower',
       'yearly_upper', 'multiplicative_terms', 'multiplicative_terms_lower',
       'multiplicative_terms_upper', 'yhat']
'''

# 시각화 (plotly로 동적시각화도 가능)
m.plot(forecast)
m.plot_components(forecast)

# 성장 총량

m = Prophet(growth='logistic')
# growth 파라미터가 결국 위 식의 g(t) 개형을 결정
m.fit(df)

# 최소 최대를 설정: 로지스틱 회귀의 경우에만 적용이 가능
df['floor'] = 1.5
future['cap'] = 6
future['floor'] = 1.5
m = Prophet(growth='logistic')
m.fit(df)
fcst = m.predict(future)
fig = m.plot(fcst)

# Trend changepoints
from prophet.plot import add_changepoints_to_plot
fig = m.plot(forecast)
a = add_changepoints_to_plot(fig.gca(), m, forecast)

# 변화점 직접 명시하는 방식
m = Prophet(changepoints=['2014-01-01'])
forecast = m.fit(df).predict(future)
fig = m.plot(forecast)

 

- cap, floor가 어떻게 범위를 결정짓는지?

로지스틱으로 함수 개형이 설정되어있고, cap이 지정된 경우 위와 같은 방식으로 범위 밖의 변수를 조절한다.