시계열 딥러닝은 구간추정값 없이 오직 점추정값을 출력하고 연산과정에서의 설명력이 떨어지므로 신뢰성이 낮다는 단점이 있습니다. 그런데도 불구하고 시계열 딥러닝을 이용하는 이유는 뭘까요? 대략 세가지의 장점이 있습니다.
1. Feature Engineering을 자동으로 구해줍니다.
히든레이어를 통해서 데이터 안에 숨겨진 패턴을 찾아낼 수 있습니다. 히든레이어의 층이 깊어 질수록 1차,2차,3차 쭉죽 늘어나게 됩니다. 그래서 3차식에서 알수 있는 특징이 있고, 4차식에서 알수 있는 특징을 찾아 낼수 있습니다.
2. 입력& 출력의 자유도 높습니다.
입력을 어떻게 줄건지, 출력은 어떤 형태로 할건지 사용자가 정할수 있는 자유도가 있습니다.
3. 시퀀스 길이를 제한없이 다양한 패턴으로 적용할수있다.
시퀀스길이를 사용자 마음대로 조절이 가능합니다. 4일치 데이터로 미래를 예측할건지 60일치 데이터로 예측할건지 정할수가 있다는 거죠 게다가 높은 정확도도 딥러닝을 사용 할수 밖에 없는 이유 중 하나입니다.
본격적으로 시계열 딥러닝의 알고리즘에 대한 소개를 하겠습니다.
Recurrent Neural Network: RNN
입력층 > 은닉층 > 출력층으로 연결된 단방향 신경망 외에 이전 출력값이 다시 입력으로 연결되는 순환신경망
RNN은 AR을 적용한 딥러닝 버전입니다. AR은 자기 자신의 이전 값을 입력으로 사용합니다. 이와 마찬가지로 RNN도 자기 자신의 이전의 값을 가져오는데, 여기서 차이점은 AR은 데이터를 가져오고 RNN은 히든레이어 정보를 가져옵니다.
왜 차이가 나는걸까요? RNN의 입력값은 batch사이즈 만큼 들어옵니다. 따라서 매번 층을 돌릴때마다 데이터를 누적시키는건 상식적으로 생각해도 비효율적입니다. 이 batch사이즈 만큼 들어온 데이터를 대표할 만한 값이 있어야 하는데 이를 히든레이어가 합니다. 히든레이어에 데이터의 특징 정보가 있기 때문에 다음 히든레이어를 계산 시 함께 연산하면 됩니다. 위 이미지를 보면 녹색 동그라미가 이에 해당됩니다. 주황색의 f() 는 비선형으로 만들어 주는 활성함수 입니다.
학습자료
역전파(Backpropagation Through Time: BPTT)
앞으로 전진하는 Feed Forward Network(FFN)과 달리 각 시점마다 Loss Function에서 추정된 Error를 이전 시점으로 전파함
공통 파라미터: 각 시점마다 출력 희망값(Y^)과 실제 출력값(Y)의 Error로 가중치(w)들을 업데이트 함
학습자료
역전파를 사용해서 가중치를 업데이트 할수 있습니다. 하지만 이는 또 다른 문제를 야기 시킵니다. 저번 포스팅에서 언급했었듯이 미분은 곱셈의 연산이 있기 때문에 0.5 *0.5 = 0.25 가 되는 것처럼 점점 값이 작아져서 잔차가 0이 됩니다.
그렇기 때문에 전파가 되는 양이 점점 작아지는 문제가 생깁니다.
LeRu를 사용을 지양하는 이유
-> RNN계열은 이전 값들이 사용되므로 한번 0으로 빠지면 이전 정보가 손실이 되기때문에 tanh를 여전히 사용한다
Long Short-Term Memory: LSTM
학습자료
곱하기로 인한 정보 손실을 막고자 덧셈을 이용한 알고리즘입니다.
C라는 기억셀을 추가하여 이를 바탕으로 외부 계층에 은닉상태H를 출력합니다. C는 LSTM 계층 내에서만 전달되며 과거부터t까지 모든 정보 저장됩니다.
Ft = sigmoid(새로운 X데이터* 가중치 W + 이전 히든레이어의 정보값*가중치 U)
이 Ft에다가 이전의 모~든 정보를 기억하고있는 Ct-1를 곱해 줌으로써 정보를 잊게 합니다.
정리: RNN은 시간이 지나면 이전 값을 잊어버리는 반면, LSTM은 각 게이트를 통해 이전 입력값이 계속 저장되어 필요한 시점에 출력 및 반영됨
장기기억셀(C)는 덧셈으로 되어있기 때문에 곱셈에 의해 작아지는 전파력이 감소하지 않게 됩니다.
1-U는 위 설명처럼 얼마나 잊을지 결정하는데, U가 1이면 (새로운 입력정보를 모두 사용하겠다) ht 게이트는 작동을 할 것이고 forget을 담당하는 1-U는 0 이 되므로 기능을 발휘 하지 못합니다. 즉 잊을 부분이 없다는 겁니다. U을 구성요소를 보면 이전 히든레이어값이 들어 있습니다. 이 히든레이어에 t-1정보가 있을것입니다. 따라서 과거 정보와 현재 입력으로 들어가는 정보를 모두 사용한다면 1이 될것이고, 버릴 부분이 없다는 의미로 1-U =0이 되어서 forget할 부분이 없다는 의미 입니다.
만일 U의 비율이 0.5라면 어떤 해석이 가능할까요
U에서는 sigmoid(새로운 정보+이전 정보)의 결과가 0.5 이므로 새 정보+ 이전 정보에서 절반은 반영하고 절반은 버린다는 의미로 해석될것입니다. 이때 ht는 절반만 들어간 U * reset gate + 새로운 정보
절반이 남은 1-U는 ht-1과 곱해집니다. 만일 1-U상태로 출력이 된다면 vanishing Gradients 문제는 해결되지 않습니다.
왜냐면 다음 t+1 시점에서 또 forget이 작동하기 때입니다. 따라서 새로운 정보를 '더하기' 해줘야 하는데 이때 ht로 더해줍니다.
오늘은 딥러닝에 대해서 포스팅 하겠습니다. 발전 개요는 생략했고 중요하다고 생각하는 것들만 정리해 보았습니다.
단층퍼셉트론(Single-Layered Perceptron: SLP)
위 이미지는 단층 구조입니다. H0와 이를 가리키는 X들만 따로 구별해서 본다면 우리가 알고 있는 회귀분석 y=f(x) 구조와 같습니다. 즉 ho=f(x) 구조 형식을 띕니다. 즉 각 노드가 하나의 회귀분석이며 딥러닝은 이러한 회귀분석을 여러개 계산하는 것입니다. H0 -> Y으로 넘어 갈때 활성함수가 들어가게됩니다. 활성함수는 비선형으로 만들어 주는 장치입니다. 따라서 정리하면 뉴런(노드)의 수만큼 회귀분석 한 다음에 이를 결합 하고 또다시 회귀 분석을 합니다
Y= g(H0, H1, ... , Hk-2) 이때 g는 활성함수입니다.
히든레이어가 1개밖에 없다는 의미는 오직 선형으로만 구분 할수 있다는 것입니다. 아래 이미지를 보겠습니다
위 예시로는 한 직선으로 구분이 가능합니다. 하지만 ...
XOR 같은 문제는 한 직선으로 구분을 못하게 됩니다. 이 문제를 구분 할려면 비선형이 필요하게 되고 이를 만들기 위해서는 히든레이어를 2개 이상 쌓아야 합니다. 아래 1번 문항이 이에 해당됩니다.
2번의 의미는 가중치 업데이트가 안된다는 것입니다. 어찌저찌 해서 회귀분석 수십개를 돌린 후 Y라는 결과를 한번 딱 냈다고 했을때, 실제 Y와 예측치 Y 간의 잔차를 토대로 잔차가 줄어드는 방향으로 가중치를 조절해야 하는데, 단층에서는 이 조절할 장치가 없습니다. 초기에 정해진 가중치가 업데이트가 되지않고 유지가 되어서 최적의 값을 찾지 못합니다.
다층퍼셉트론(Multi-Layered Perceptron: MLP)
XOR문제를 해결하기 위해서 2개 이상의 히든 레이어를 쌓은 형태입니다. 하지만 여전히 가중치 조절은 해결하지 못한 상태입니다.
히든 레이어가 1개면 1차식, 2개면 2차식, 3개면 3차식 형태로 증가합니다. 따라서 다양한 비선형 특징을 추출할수 있게 됩니다. 이로써 XOR 문제는 해결이 가능해 졌습니다. 이제 가중치 업데이트를 어떻게 하면 될건가 문제인데
이는 역전파를 통해서 해결이 가능합니다.
역전파(Back Propagation: BP)
실제 Y와 예측치 Y의 차이인 error를 학습 시 사용했던 파라미터의 수만큼 잘게 쪼갭니다. 그러면 위 이미지의 노란색 화살표에 그려진 e값이 분산이 됩니다. 정리를 하자면, 단층 퍼셉트론 설명할때 노드의 수만큼 회귀분석이 실행이 된다고 했는데, 이 각각의 회귀분석에게 잔차를 나눠주고 각자 할당 받은 잔차를 가장 낮게 가중치를 조절하는 것입니다. 마치 회사 대표가 이번년도 매출은 100억원이다 라고 정해 놓으면 각 부서의 차장이나 책임자들이 목표를 소분해서 가장효율적인 방법을 찾아 일하는 것과 비슷해 보입니다. 역전파를 통해 학습 업데이트를 해결 할수 있었습니다.
하지만 이는 또다른 문제를 야기 했습니다. 이때 사용되는 활성함수는 sigmoid인데, 바로 vanishing Gradients가 발생하는 것입니다.
히든레이어가 많아지면 잔차 값이 0이 됩니다. 그러면 가장 처음부분에서 계산된 파라미터, 여러개의 회귀분석은 잔차값을 받지 못하여 업데이트를 못하게 됩니다. 역전파 과정에서 미분끼리 연산과정이 발생합니다. 예를들어 0.5, 0.5 이 둘을 곱하면 0.25로 줄어 들듯 뒤로 갈수록 잔차가 줄어들고 결국 0밖에 남지 않는다는 것이죠
이 문제는 LeRu라는 활성함수를 통해 해결할수 있습니다.
딥러닝(Deep Learning: DL)
인공신경망의 다층구조가 심화(Deep)된 네트워크 구조의 알고리즘
(시계열)회귀분석과 (딥러닝)인공신경망의 차이
y=f(x)의 형태를 취하고 있는 반면 인공신경망은 중간에 히든레이어가 추가된 형태입니다.
Y를 구하기 위한 연관 과정을 위 식에서 확인 할 수 있습니다. Vi * h(x) 에서 h(x)는 단층에서 언급했었던 내용처럼 히든레이어를 회귀식으로 취했을때 형태입니다. h(x) 안에는 X들과 가중치 값들을 포함하고 있습니다.
위 내용의 이유로 현재 딥러닝을 이용한 비지니스 수익창출하기가 아직은 어렵다고 합니다. 딥러닝은 사람의 손을 많이 탑니다. 자유도가 높아서 값이 매번 바뀌고 가중치의 설명력이 부족하기때문에 정확성이 높아도 신뢰하기가 어렵습니다.
두 비정상성 시계열을 선형조합(Linear Combination)하여 생성한 시계열의 적분(Integration) 차수가 낮아지거나 정상상태가 되는 경우입니다.
위 정의가 어떤 의미인가를 생각해 보자
일단 두 비정상성 시계열 데이터가 필요해 보이고 이 둘 사이에서 뭔가 ''공통''인 부분이 있습니다. 결론부터 말씀드리자면 이 공통이어야 하는 것은 방향입니다. 보통 술에 취한 사람과 개로 예를 듭니다.
술 취한 사람은 비틀비틀 걸어 다녀서 어디로 튈지 모릅니다. 가다가 갑자기 넘어질수도 있고 중심을 잃어서 갑자기 왼쪽으로 휙 방향을 틀 수도 있습니다. 만일 술 취한 사람이 두 명이라면 이 두 명은 각자 가는 방향이 다를 것입니다. 여기서 술 취한 사람이 목줄 하고 있는 개가 있다면 얘기가 달라집니다. 개도 자기 페이스대로 이리저리 왔다 갔다 하겠지만 어찌 됐든 술 취한 주인을 따라가게 되어 있습니다. 이때 애완견의 위치를 나타내는 시계열과 취한 사람의 위치를 나타내는
시계열은 공적분이라고 합니다.
공적분이 되기위한 조건은 아래와 같습니다
각각의 시계열들이 모두 같은 차분의 수(order of integraion)를 가진다. (조건1)
시계열들의 선형결합으로 만들어진 새로운 시계열은 기존의 시계열들보다 더 낮은 차분의 수를 가져야 한다.(조건 2)
추가적인 설명을 하자면
Y1(술 취한 사람 위치 시계열)과 Y2(애완견 위치 시계열)을 각각 봤을 때 비정상 시계열 데이터입니다. 정상성으로 만들기 위해서 차분이 각각 2씩 필요해야 합니다. (조건 1: 차분의 수가 같아야 한다)
또한 Y1과 Y2를 선형 결합을 했을 때 만들어지는 비정상 시계열을 차분하게 될 때 차분의 수가 2보다 미만이어야 합니다(조건 2: 기존 시계열 차분의 수 보다 작아야 한다)
만일 애초에 차분의 수가 하나씩 필요했다면 선형 결합으로 인해 차분의 수가 0이어야 한다. 즉 선형 결합을 하면 정상성 데이터가 된다는 뜻입니다.
차분의 수가 각자 1일 경우에 대해서 좀 더 자세히 말씀드리겠습니다
강의 자료
둘 다 비 정상성인 랜덤워크가 있을 때 각자 방정식은 위와 같습니다
이 두 개를 선형 결합한 대문자 Y의 수식 풀이를 보면 결국 잔차만 남게 될 때 이 잔차가 WN(화이트 노이즈)라면 차분이 필요가 없게 되므로 Y1, Y2는 공적분 관계입니다.
이를 실습으로 확인해 보겠습니다
실습: 공적분의 이해
# 데이터 생성
np.random.seed(1)
Y1= np.random.randn(500).cumsum() #차분을 하면 정상성 -> diff, 누적을 하면 비정상성 -> cumsum
Y2= 0.6*Y1 + np.random.randn(500) #랜덤워크의 60% + 추가적인 랜덤을 추가(np.random.randn(500)) => Y1과 스케일 차이가 발생
Y = 0.6 * Y1 - Y2 #선형결합 (왼쪽 식 계산결과 e만 남는다. 즉 잔차는 WN)
#시각화
plt.figure(figsize=(10,4))
plt.plot(Y1, 'g', label ='Y1')
plt.plot(Y2, 'b', label ='Y2')
plt.legend()
plt.show()
plt.figure(figsize=(10,4))
plt.plot(Y, label='Y * 0.6*Y1 - Y2')
plt.legend()
plt.show()
#VECM 공적분 확인
coint_result = sm.tsa.coint(Y1,Y2)
pd.DataFrame([coint_result[0], coint_result[1]], index=['statistics', 'p-value'], columns=['output'])
오늘의 포스팅 내용은 다변량 선형확률입니다. 지금껏 진행 한 것은 Y가 1개인 경우였습니다. 다변량은 Y가 2개 이상인 것을 말합니다. 매출과 광고비가 있을때 단변량에서는 매출이 Y라고 정해놓고 시작했습니다. 하지만 반대로 광고비가 종속변수로써 영향을 더 많이 받는다면 말이 달라지게 됩니다. 따라서 어떤 변수가 가장 영향을 많이 받는 종속변수 인지를 알고자 할때 필요한게 단변량입니다.
단변량 알고리즘 중에서 대표적으로 벡터자기회귀 모형을 먼저 보겠습니다.
벡터자기회귀 모형(Vector Autoregressive Model)
교재에서는 위 내용을 VAR의 정의로 적혀 있는데 직관적으로 이해하기 힘들어서 구글링을 해봤습니다.
나의 과거에만 영향을 받는 모델: 자기회귀
나의 과거 뿐만 아니라 다른 변수의 과거에서도 영향을 받는 모델: 벡터자기회귀
다른 변수의 과거까지 수식으로 함께 넣다보니까 벡터 형식으로 표현이 가능합니다.
VAR(1)인 경우: 각 변수 한 타임 과거
Yt-1안에 A라는 변수 At-1과 B라는 변수 Bt-1이 함계 들어 있어서 이를 간단하게 표현한게 위 식입니다.
At-1,Bt-1 계수를 벡터형식으로 모아 놓은게 A1입니다. 참고로 꺽쇠가 있다는 것은 그 안에 여러개가 들어 있다는 뜻입니다.
import numpy as np
import matplotlib.pyplot as plt
import statsmodels
import statsmodels.api as sm
#차수 입력
intercept = np.array([5,3])
matrix_A = np.array([[0.2,0.3], [-0.6,1.1]])
residual_covarinace = np.array([[1,0.8], [0.8, 2]])
# VAR 데이터 생성
fit = statsmodels.tsa.vector_ar.var_model.VARProcess(matrix_A, intercept, residual_covariance)
#시뮬레이션 시각화1
simul_num = 100
fit.plotsim(steps=sumul_sum, sseed=123) #100개 만큼 임의로 생성 후 시각화
plt.tight_layout()
plt.show() #첫번째, 두번째 그래프
# 시뮬레이션 시각화2_fit의 계수를 기반으로 시뮬레이션
# 임의로 데이터 두개를 시각화한 사실이 중요하다
# 예를들어 하나는 한국주식시장, 다른 하나는 미국 주식시장
simul_num =100
#sigma_u = corr과 같음
simul_values = statsmodels.tsa.vector_ar.util.varsim(fit.coefs, fit.intercept, fit.sigma_u, steps=simul_num)
plt.figure(figsize=(10,5))
plt.plot(simul_values)
plt.tight_layout()
plt.show() #두개가 겹친 그래프
#ACF 시각화
fit.plot_acorr()
plt.tight_layout()
plt.show()
'''
대각: Y1(1,1)에 대한 자기상관, Y2(2,2)에 대한 자기상관 (둘다 자기상관이 있음)
반 대각: Y1,Y2 서로에 대한 corr
'''
#VAR 모형적합
fit = sm.tsa.VAR(simul_values).fit()
display(fit.summary())
# 예측 및 시각화
# endog: 각 시점 마다의 계수가 출력된다
forecast_num = 20
pred_var =fit.forecast(fit.model.endog[-1;], steps = forecast_num)
pred_var_ci = fit.forecast_interval(fit.model.endog[-1:], steps=forecast_num)
fit.plot_forecast(forecast_num)
plt.tight_layout()
plt.show()
#잔차진단
fit.plot_acorr()
plt.tight_layout()
plt.show()
'''
각각의 autocorr이 작기때문에 화이트노이즈에 가깝다는 추론 할수있다.
'''
Results for equation이 y1과 y2 두개가 있습니다. y1을 한국시장, y2를 미국 시장이라고 했을때 방정식이 두개로 나뉘어서 어떤 변수가 종속변수로써 더 적합한지를 분석하게 됩니다. (equation=방정식)
한국시장은 0.085의 유의하지 않는게 하나 발견이 되었고 미국시장은 모두 유의하다고 나왔습니다. 따라서 미국시장이 다른 변수에 의해 많은 영향을 받고 있기 때문에 한국시장보다 미국시장이 종속변수로써 더 적합하다고 할수 있습니다.
실선+도트는 신뢰구간이고 실선만 있는 것은 점추정 값입니다. 마지막 값이후로 steps 만큼 예측을 하게 됩니다.
임펄스 응답 함수(Impulse Response Function)
VAR 모형은 여러개의 시계열 상호상관관계를 기반으로 각각의 변수가 다른 변수에 어떤 영향을 주는지 임펄스 반응 함수로 알 수 있습니다. 즉 위 실습에서 y1은 한국주식시장, y2는 미국주식시장 이라는 가정에서 y1이 변화를 했을때 y2에 미치는 영향을 표시한 것입니다. 실습코드는 위 실습 내용과 이어짐을 알려 드립니다
# 임펄스반응함수 추정
#irf: impulse responsible function
#20라는 건 20lag을 뜻함
fit.irf(forecast_num).plot()
plt.tight_layout()
plt.show()
y1이 변화할때 y2에 끼치는 영향은 마이너스를 보이고 반대로 y2가 변했을때 y1에 끼치는 영향은 플러스입니다.
만일 주식시장아니라 소비(y1), 매출(y2)라고 했을때 소비가 줄면 매출도 감소하는 효과가 있다 라고 해석할수 있습니다.
실습2: 거시경제 VAR 모형화 및 해석
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import statsmodels
import statsmodels.api as sm
#데이터 로딩
raw = sm.datasets.macrodata.load_pandas().data
dates_info =raw[['year','quarter']].astype(int).astype(str)
raw.index = pd.DateTimeIndex(sm.tsa.datatools.dates_from_str(dates_info['year']+'Q'+date_info['quarter']))
raw_use = raw.iloc[:,2:5]
#데이터 시각화
raw_use.plot(subplots=True, figsize(12,5))
plt.tught_layout()
plt.show()
raw_use.diff(1).dropna().plot(subplots=True, figsize(12,5))
plt.tight_layout()
plt.show()
# VAR 모형적합
#추세를 제거한게 정상화가 될 확률이 높기때문에 이걸로 VAR모형적합을 한다
raw_use_return = raw_use.diff(1).dropna()
fit = sm.tsa.VAR(raw_use_return).fit(maxlags=2) #p=2
display(fit.summary())
'''
realcons(소비)이 가장 유의한 변수가 많다. 6가지 변수가 realcons에 영향을 준다는 뜻이며
가장 종속성이 짙은걸 알수있다.
'''
#예측 및 시각화
forecast = 20
fit.plot_forecast(forecast_num)
plt.tight_layout()
plt.show()
#임펄스반응함수 추정
fit.irf(forecast_num).plot()
plt.tight_layout()
plt.show()
# 잔차진단
fit.plot_acorr()
plt.tight_layout()
plt.show()
'''
튀는 값이 있긴하지만 전반적으로 WN에 가깝다는 것만 알고 넘어간다다
'''
차분 전 그래프차분 후 그래프
1.realcons, realinv는 특정 시차동안 꾸준하게 양수의 효과를 보인다. 2.realgdp->realcons: gdp가 1만큼증가시 realcons가 오히려 감소한다. 즉 gdp가 증가한다고 해서 바로 소비로 이어지지 않는다 오히려 줄인다 3. realinv->realcons: 투자가 늘면 소비도 증가한다. 4.realcons -> realgdp: 소비를 늘리면 gdp가 1이상 늘어나는 긍정적인 효과를 보인다. 그래서 gdp를 늘리려면 realinv를 늘리는 것 보다는 소비를 늘리는게 효과가 좋다
X,Y 둘간의 관련성이 없는데 단순히 상관성 있다는 것이 마치 인과성이 있다고 해석하는게 Spurious(거짓된, 겉으로만 그럴싸한) Regression 용어이다.
그래프상에서 아이스크림과 상어공격이 언뜻 상관성이 있어보인다. 아이스크림 매출 증가시 상어 공격 건수도 증가하므로 이 둘은 인과관계가 있다고 착각하는게 Spurious Regression이다.
단지 여름이라서 아이스크림 판매와 상어 공격이 많아진것이 올바른 인과관계이다.
이를 더 확장해서 생각해 보면 실생화의 징크스도 해당이 된다.
징크스: 머리를 염색했더니 시험점수가 100점이 나오더라
머리염색을 한 후에 시험을 치뤘더니 100점이 나왔다는 징크스얘기는 머리염색이 원인이라는 인과관계 형태를 보인다.
머리염색과 시험성적의 인과관계를 우리는 규명할수가 없다. 하지만 정말로 머리염색을 한 후 (선행) 의 시험성적은 영향을 받는지에 대한 인과관계는 확인해 볼수 있다. 다시 말하면 인과관계는 파악 못하지만 시간 흐름에 따른 인과관계를 파악하고자 할때를 Granger Causality라는 개념을 사용해서 확인할수 있다. 아래 예를들어본다.
추론불가한 문제:"닭이 먼저인가 달걀이 먼저인가?" (인과관계)
추론가능한 문제:"닭과 달걀의 생성순서 별 서로의 영향력은 어떤가?" (Granger 인과관계)
Granger Causality를 사용하기 위한 조건들이 있다.
강의 자료
just use Y식을 보면 X(다른 변수)가 없고 오직 자기자신의 과거만 해당이 된다. 현재는 과거의 나에게 얼만큼 영향을 주는지를 알수있고
use X and Y에서 베타X 가 있는데 위 이미지처럼 Y가 달걀 X가 닭이라고 했을때 달걀의 생산량은 닭의 생상량에 얼마큼 영향을 받고 있는지를 알수가 있다.
이해를 돕기위해서 아래 실습을 실행한다.
실습: 닭과 달걀의 생산량을 통한 Granger 인과관계 이해
# 닭과 달걀을 시간적 순서로 봤을때 어떤게 선행성이 있는가
import pandas as pd
import matplotlib.pyplot as plt
import statsmodels.api as sm
#데이터 로딩
location = './Data/ChickenEggProduction/Chicken_Egg_Production.txt'
raw_all = pd.read_csv(location, sep='\t')
raw_all.head()
#데이터 시각화
plt.figure(figsize=(10,4))
plt.plot(raw_all['Chicken'], 'g', label='Chicken')
plt.legend()
plt.show()
plt.figure(figsize=(10,4))
plt.plot(raw_all['Egg'], 'b', label='Egg')
plt.legend()
plt.show()
# 정상성 변환시각화
# 정상성을 가져야 실행이 되므로 diff를 적용
plt.figure(figsize=(10,4))
plt.plot(raw_all['Chicken'].diff(1), 'g', label='Diff(Chicken)')
plt.legend()
plt.show()
plt.figure(figsize=(10,4))
plt.plot(raw_all['Egg'].diff(1),'b',label='Diff(Egg)')
plt.legend()
plt.show()
# Granger Causality 테스트
print('\n[Egg -> Chicken]')
#maxlag 몇개의 파라미터를 적용할것이가(lag값)
#addcons 상수항 적용 여부
#맨 처음 x의 입력은 array 형식이여야 하므로 values
'''
알고리즘 설계상 두번째 컬럼이 첫번째 컬럼에 영향을 주는지를 분석하므로 달걀이 두번째 컬럼으로 와야한다.
iloc[:,1:] -> 컬럼이 닭,달걀순이므로 두번째 컬럼인 달걀이 닭에 Granger Causality가 있는지를 본다
'''
granger_result1 = sm.tsa.stattools.grangercausalitytests(raw_all.diff(1).dropna().iloc[:,1:].values, maxlag=3, verbose=True)
'''
lag1~3 통계량을 봤을때 p-value가 5%보다 작으므로 X가 Y에 영향을 준다
'''
print('\n[Chicken -> Egg]')
#순서 바꾸기 닭->달걀
granger_result2= sm.tsa.stattools.grangercausalitytests(raw_all.diff(1).iloc[:,[2,1]].values, maxlag=3, verbose=True)
'''
전혀 도움이 되고 있지 않다.
즉, 달걀->닭이다.
'''
# 의사결정
# 닭이 달걀을 낳으면 그 수는 약 3년후까지 닭의 수에 영향을 준다
# 왜 3년? lag3까지 유의하게 나왔기 때문
# 닭의 수가 많아진다고해서 달걀을 많이 낳지는 않는다
# 달걀 -> 닭 (Granger Causality) 달걀이 닭에 Granger Causality 경향을 보인다.
차분 후차분 후
달걀->닭 인 경우 lag3까지 모두 유의하며, 닭-> 달걀 인 경우 lag1까지만 유의하다
import timeit
start = timeit.default_timer()
from fbprophet import Prophet
from fbprophet.plot import plot_plotly, plot_components_plotly
from fbprophet.plot import add_changepoints_to_plot
# Rearrange of Data
Y_train_feR_prophet = Y_train_feR.reset_index() #시간축을 인덱스로 가지고 있으면 안된다.
#Y값과 시간정보가 들어간 컬럼이 있어야 한다.
#컬럼명도 아래처럼 입력해야 작동된다.
Y_train_feR_prophet.columns =['ds','y'] #Prophet실행 시 무조건 해야 하는것
#하나의 데이터프레임으로 만들어서 fit 함수에 적용
YX_train_prophet = pd.concat([Y_train_feR_prophet, X_train_feRSM.reset_index().iloc[:,1:]],
axis=1)
# Applying Prophet Model
#change point: 그래프에 변화를 보일때 시점(CP)
#선형적인 효과로 추세가 있다.
fit_reg1_feRSM_prophet = Prophet(growth='linear',#추세가 선형인지, 2차함수인지 로그함수인지 임의로 지정
'''
여기서 말하는 trend를 이해하기 쉽게 설명.
한 회사의 사장님이 구속되는 일이 발생하면 주식(trend)은 상향에서 하향으로 변화 할것이다.
변화하는 시점을 CP로 본다
'''
'''
trend라는 건 잘개 쪼개진 것들의 합으로 구성된 것
잘개 쪼깨는 기준이 CP가 된다
이 알고리즘이 원하는 것: CP 정보만 주면 예측하겠다는 것이다.
이때, CP의 정보가 아래이다.
'''
#미래든 과거든 이벤트가 발생한 날을 적용
changepoints=None, # CP가 발생하는 시점들의 list ['2012-01-01']
#예상되는 cp를 넣으면 해당 범위에서 분석한다(max 값)
n_changepoints=25, # CP의 수
#CP값을 비율로 나타낸것
# 1000개 중 0.8, 즉 최대로 나올수 있는 CP는 800개 라는 뜻
changepoint_range=0.8, # CP의 기존 데이터 수 대비 최대 비율
# 퍼센트 개념x, 많을수록 CP를 많이 추정하려고 한다.
changepoint_prior_scale=0.05, # CP 추정 민감도로 높을수록 민감
'''
Seasonality
'''
#계절성이 y값에 더하기로 표현되는지 곱하기로 표현되는지 지정해준다.
seasonality_mode='additive', # 계절성 모델: 'additive' or 'multiplicative'
#계절성 가중치를 결정
seasonality_prior_scale=10.0, # 계절성 추정 민감도로 높을수록 민감
#데이터를 보고 자동으로 일~연계절성을 자동으로 추정해준다
#없으면 추정 안한다.
#auto: 있으면 추정 없으면 추정 안함
yearly_seasonality='auto', # 연계절성
weekly_seasonality='auto', # 주계절성
daily_seasonality='auto', #일계절성
# Holiday
# 2021-12-25 크리스마스 <-이런형식으로 데이터프레임으로 입력
holidays=None, # 휴일 또는 이벤트 시점 dataframe
holidays_prior_scale=10.0, # 휴일 추정 민감도로 높을수록 민감
'''
Others
'''
#미래에 나올 예측구간을 설정한다. -> 기존의 범위(train)에서 80% 반영한 범위를 사용하겠다
interval_width=0.8, # 추세 예측 정확도 구간범위
mcmc_samples=0) # 계절성 예측 정확도 제어(몬테카를로) changepoints =None
fit_reg1_feRSM_prophet = fit_reg1_feRSM_prophet.fit(YX_train_prophet)
#예측할 데이터의 freq가 뭔지, 몇개 시점을 예측해야 하는지를 입력
#실제 예측한다는 용도보다는 freq가 적당한지를 검증하기 위한 정도->도메인 전문가에게 검증을 받는 용도
forecast = fit_reg1_feRSM_prophet.make_future_dataframe(freq='H',
periods=Y_test_feR.shape[0])
pred_reg1_feRSM_prophet = fit_reg1_feRSM_prophet.predict(forecast) #train,test 예측값을 함께 출력한다.
#붙어있는 train,test를 분리한다
pred_tr_reg1_feRSM_prophet = np.ravel(pred_reg1_feRSM_prophet.loc[:Y_train_feR.shape[0]-1,['yhat']])
pred_te_reg1_feRSM_prophet = np.ravel(pred_reg1_feRSM_prophet.loc[Y_train_feR.shape[0]:, ['yhat']])
# Result
display(pred_reg1_feRSM_prophet[['ds'.'yhat','yhat_lower','yhat_upper']])
#trend 시각화
fig =fit_reg1_feRSM_prophet.plot(pred_reg1_feRSM_prophet) #fit한 것과 예측한 것을 같이 시각화fit한 것과 예측한 것을 같이 시각화
#gca: fig안에 라벨값, 폰트 값 등 셋팅 값을 가지고있음
add_chnagepoints_to_plot(fig.gca(), fit_reg1_feRSM_prophet, pred_reg1_feRSM_prophet)
plt.xticks(rotation=90)
plt.show()
#seasonal 시각화
# plot_components: 일,연,달 계절성, 예측값과 관련된 값이 나온다
fig = fig_reg1_feRSM_prophet.plot_components(pred_reg1_feRSM_prophet)
ax = fig.get_axes()
#아래처럼 하나씩 작업해야 시각화가 됨
for i in range(len(ax)):
ax[i].tick_params(axis='x', labelsize=15, rotation=90)
plt.show()
## using plotly
# display(plot_plotly(fit_reg1_feRSM_prophet, pred_reg1_feRSM_prophet))
#plot_components_plotly: 핏팅 값과 예측값을 각각 시각화
# display(plot_components_plotly(fit_reg1_feRSM_prophet, pred_reg1_feRSM_prophet))
# Evaluation
Score_reg1_feRSM_prophet, \
Resid_tr_reg1_feRSM_prophet, \
Resid_te_reg1_feRSM_prophet = evaluation_trte(Y_train_feR,
pred_tr_reg1_feRSM_prophet,
Y_test_feR,
pred_te_reg1_feRSM_prophet,
graph_on=True)
display(Score_reg1_feRSM_prophet)
# Error Analysis
error_analysis(Resid_tr_reg1_feRSM_prophet, ['Error'], X_train_feRSM, graph_on=True)
stop = timeit.default_timer()
print('Time: ', stop - start)
검정색이 실제 train값, 파란색이 train,test 예측값 주황색: 점추정값 빨간색 점선: CP
각 해당되는 빈도로 계절성을 뽑은 패턴입니다. 이를 통해 각 빈도의 흐름을 파악할수 있습니다.
본 강의에서는 trend와 seasonal 파라미터 조정 시 어떻게 변화하는지에 대한 실습이 많이 있지만 그중에서 최종적으로
적용된 최종 실습 하나만 올리도록 하겠습니다
Final
import timeit
start = timeit.default_timer()
from fbprophet import Prophet
from fbprophet.plot import plot_plotly, plot_components_plotly
from fbprophet.plot import add_changepoints_to_plot
# Rearrange of Data
Y_train_feR_prophet = Y_train_feR.reset_index()
Y_train_feR_prophet.columns = ['ds', 'y']
YX_train_prophet = pd.concat([Y_train_feR_prophet,
X_train_feRSM.reset_index().iloc[:,1:]],
axis=1)
# Applying Prophet Model_비선형적 추세를 넣기 위한 작업
### saturating ###
saturating = 'logistic' #logistic를 쓸때는 Y값의 상한과 하한을 컬럼으로 생성해야한다.
if saturating == 'logistic':
#예상되는 train 상한과 하한을 정해준다
YX_train_prophet['cap'] = 900
YX_train_prophet['floor'] = 0
### holiday ### -> 수동으로 넣을때 방법
newyear = pd.DataFrame({'holiday': 'newyear',
'ds': pd.to_datetime(['2011-01-01', '2012-01-01']),
'lower_window': 0, # 이전날짜 포함범위
'upper_window': 1}) # 이후날짜 포함범위
christmas = pd.DataFrame({'holiday': 'christmas',
'ds': pd.to_datetime(['2011-12-25', '2012-12-25']), # 미래 시점도 포함
'lower_window': 0, # 이전날짜 포함범위
'upper_window': 1}) # 이후날짜 포함범위
holidays = pd.concat([newyear, christmas])
fit_reg1_feRSM_prophet = Prophet(growth=saturating,
# Trend
changepoints=None, # CP가 발생하는 시점들의 list ['2012-01-01']
n_changepoints=25, # CP의 수
changepoint_range=0.9, # CP의 기존 데이터 수 대비 최대 비율
changepoint_prior_scale=0.5, # CP 추정 민감도로 높을수록 민감
# Seasonality
seasonality_mode='additive', # 계절성 모델: 'additive' or 'multiplicative'
seasonality_prior_scale=20.0, # 계절성 추정 민감도로 높을수록 민감
yearly_seasonality='auto', # 연계절성
weekly_seasonality='auto', # 주계절성
daily_seasonality='auto', #일계절성
# Holiday
holidays=holidays, # 휴일 또는 이벤트 시점 dataframe
holidays_prior_scale=5.0, # 휴일 추정 민감도로 높을수록 민감 (이전결과참고반영)
# Others
interval_width=0.8, # 추세 예측 정확도 구간범위
mcmc_samples=0) # 계절성 예측 정확도 제어
### holiday ###
#US를 적으면 미국의 휴일,공유일이 자동으로 인식된다.(한국은 없음)
fit_reg1_feRSM_prophet.add_country_holidays(country_name='US')
### monthly seasonality(계절성 추가하기)###
#성능이 안좋아져서 주석처리 후 위에 holiday 수동으로 직접 넣음
# fit_reg1_feRSM_prophet.add_seasonality(name='monthly', period=24*7*4, fourier_order=1)
### extra feature ###
# X(설명변수)들을 고려하여 예측하기
target_name = X_train_feRSM.columns[:5]
for col in target_name:
fit_reg1_feRSM_prophet.add_regressor(col)
fit_reg1_feRSM_prophet = fit_reg1_feRSM_prophet.fit(YX_train_prophet)
forecast = fit_reg1_feRSM_prophet.make_future_dataframe(freq='H',
periods=Y_test_feR.shape[0])
### saturating ###
# test 상하한 범위 지정
if saturating == 'logistic':
forecast['cap'] = 700
forecast['floor'] = 0
### extra feature ###
forecast = pd.concat([forecast,
pd.concat([X_train_feRSM[target_name],
X_test_feRSM[target_name]], axis=0).reset_index().iloc[:,1:]],
axis=1)
pred_reg1_feRSM_prophet = fit_reg1_feRSM_prophet.predict(forecast)
pred_tr_reg1_feRSM_prophet = np.ravel(pred_reg1_feRSM_prophet.loc[:Y_train_feR.shape[0]-1, ['yhat']])
pred_te_reg1_feRSM_prophet = np.ravel(pred_reg1_feRSM_prophet.loc[Y_train_feR.shape[0]:, ['yhat']])
# Result
display(pred_reg1_feRSM_prophet[['ds', 'yhat', 'yhat_lower', 'yhat_upper']])
fig = fit_reg1_feRSM_prophet.plot(pred_reg1_feRSM_prophet)
add_changepoints_to_plot(fig.gca(), fit_reg1_feRSM_prophet, pred_reg1_feRSM_prophet)
plt.xticks(rotation=90)
plt.show()
fig = fit_reg1_feRSM_prophet.plot_components(pred_reg1_feRSM_prophet)
ax = fig.get_axes()
for i in range(len(ax)):
ax[i].tick_params(axis="x", labelsize=15, rotation=90)
plt.show()
## using plotly
# display(plot_plotly(fit_reg1_feRSM_prophet, pred_reg1_feRSM_prophet))
# display(plot_components_plotly(fit_reg1_feRSM_prophet, pred_reg1_feRSM_prophet))
# Evaluation
Score_reg1_feRSM_prophet, \
Resid_tr_reg1_feRSM_prophet, \
Resid_te_reg1_feRSM_prophet = evaluation_trte(Y_train_feR,
pred_tr_reg1_feRSM_prophet,
Y_test_feR,
pred_te_reg1_feRSM_prophet,
graph_on=True)
display(Score_reg1_feRSM_prophet)
# Error Analysis
error_analysis(Resid_tr_reg1_feRSM_prophet, ['Error'], X_train_feRSM, graph_on=True)
stop = timeit.default_timer()
print('Time: ', stop - start)
위 실습은 어떻게 적용하는지에 대한 교육 초점을 뒀기 때문에 성능 개선을 위한 정밀한 작업은 없었음을 알려드립니다.
