다변량 시계열 데이터 1 (Multivariate Time Series Data)

220928 작성

<본 블로그는 today-1님의 블로그를 참고해서 공부하며 작성하였습니다 :-) >

https://today-1.tistory.com/38?category=886697 

다변량 선형 확률과정(VAR/Granger Causality/Cointegration)

 

 

 

 

1️⃣ 다변량 시계열 데이터 (Multivariate Time Series Data)

: 각 시간 단위마다 여러 개의 값을 가지는 데이터

: 다중 시간 종속 변수로 구성

: 다변량 분석에서 예측할 변수의 과거의 데이터를 고려해야할 뿐만 아니라 여러 변수들 사이의 의존성을 고려

 

 

2️⃣ 다변량 시계열 모델

다변량 선형&nbsp;확률과정

 

💕 벡터 자동 회귀 분석 VAR(Vector Auto Regression)

: 예측할 변수의 과거 값뿐만 아니라 예측할 변수와 의존성이 있는 변수들까지 고려하여 선형 함수로 나타내는 확률적 과정

: 종속 변수와 독립 변수는 상호 영향을 받는 존재

: 두 변수들 중 어떤 변수가 종속변수로 적합한지에 대한 문제를 해결하기 위해 활용

 

💕 그래인저 인과관계 (Granger Causality)

: 정상성 데이터 입력 (차분 필요)

: '닭이 먼저냐 달걀이 먼저냐' 문제를 해결할 때 사용

• 추론 불가한 문제: "닭이 먼저인가 달걀이 먼저인가?" (인과관계)
• 추론 가능한 문제: "닭과 달걀의 생성순서 별 서로의 영향력은 어떤가?" (Granger 인과관계)

: 원인과 인과 관계를 규명하는 어렵기 때문에 상대적으로 두 요인 중 먼저 영향을 미치는 변수를 알아보고자 할 때 사용

• 귀무가설(Null Hypothesis, 𝐻0H0): 한 변수가 다른 변수를 예측하는데 도움이 되지 않는다
• 대립가설(Alternative Hypothesis, 𝐻1H1): 한 변수가 다른 변수를 예측하는데 도움이 된다

 

💕 공적분 (Cointegration)

: 비정상성 데이터 입력

: 공적분 상태 = 두 비정상성 시계열을 선형조합하여 생성한 시계열의 적분 차수가 낮아지거나 정상상태가 되는 경우

: 공적분 시계열은 서로 상관관계를 가지고 있지 않더라도 장기적으로 같은 방향으로 움직이는 특성을 지님

: 페어 트레이딩 전략에 활용

 

 

 

3️⃣ 코드 구현

🔴 Library & data load

• statsmodels : 사용자가 데이터를 탐색하고 통계적 모델을 추정하며 통계적 테스트를 수행할 수 있게 도와주는 API

import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
import pandas as pd

import statsmodels.api as sm
from statsmodels.tsa.api import VAR
from statsmodels.tsa.stattools import adfuller

• 실질GDP 샘플

data = sm.datasets.macrodata.load_pandas().data
data.head()

• year, realgdp, realdpi 세가지만 사용

mydata = data[["realgdp", 'realdpi']]
mydata.index = data["year"]
mydata.head()

mydata.plot(figsize = (8,5))

상승 추세

• Stationary time series
  • 데이터가 정상성을 가진다는 의미는 데이터의 평균과 분산이 안정되어 있어 분석하기 쉬움
  • VAR을 적용시켜보기 전에 두 시계열 변수가 모두 stationay 상태이어야 함

 

🔴 AIC 기준을 사용

: 데이터의 stationarity를 찾기 위해 ADF(Advanced Dickey-Fuller test)와 같은 통계적 테스트를 수행

---

from statsmodels.tsa.stattools import adfuller

• p-value
• test-statistics
• critical value : 통계 검정에서 검정 값의 판단 기준으로 작용하는 값 (ADF Test 검정 값이 보다 낮게 되면 p-value가 낮아 귀무 가설 기각 가능)
• lag, observation

dftest = adfuller(x, maxlag, regression, autolag)

• x : 시계열 데이터
• maxlag : ADF Test 에서 p 지정 (지정하려면 autolag를 None으로 설정)
• regression
  • c : 추세는 없고 상수항 존재
  • nc : 상수항과 추세가 없음
  • ct : 추세와 상수항 둘다 존재한다고 가정
  • ctt : 상수항과 일차, 이차 추세가 모두 존재한다고 가정
• autolag
  • ADF Test 에서 p를 자동으로 지정
    • AIC, BIC : 둘중 가장 낮게 나오는 p를 자동 설정 (여기서 지정하면 maxlag 무시)
    • None : maxlag 지정값 설정
    • t-stat : maxlag에서 지정한 값부터 regression을 수행하면서 통계 검정 p-value가 5% 미만 될 때의 래그 p값으로 설정

---

adfuller_test = adfuller(mydata['realgdp'], autolag= "AIC")

print("ADF test statistic: {}".format(adfuller_test[0]))
print("p-value: {}".format(adfuller_test[1]))

adfuller_test = adfuller(mydata['realdpi'], autolag= "AIC")

print("ADF test statistic: {}".format(adfuller_test[0]))
print("p-value: {}".format(adfuller_test[1]))

• 두 경우 모두 p-value가 충분히 유의미한 값을 가지지 않아 시계열 데이터가 non-stationary

 

🔴 differencing (차분)

adfuller_test = adfuller(mydata['realgdp'], autolag= "AIC")

print("ADF test statistic: {}".format(adfuller_test[0]))
print("p-value: {}".format(adfuller_test[1]))

adfuller_test = adfuller(mydata_diff['realdpi'], autolag= "AIC")
print("ADF test statistic: {}".format(adfuller_test[0]))
print("p-value: {}".format(adfuller_test[1]))

• realgdp, realdpi 모두 p-value 값이 작아짐 => sationary 

 

🔴 모델링

• 마지막 10일은 test 나머지는 train

train = mydata_diff.iloc[:-10,:]
test = mydata_diff.iloc[-10:,:]

 

• VAR모델의 최적 순서
  • 최적의 모델을 찾기 위한 기준 AIC(Akaike's Information Criterion)를 모델 선택 기준
  • 최상의 AIC점수를 바탕으로 VAR의 순서(p)를 선택
  • AIC는 일반적으로 모델이 너무 복잡하다는 이유로 불이익을 주곤 하는데 복잡한 모델은 일부 다른 모델 선택 기준에서 약간 더 나은 성능을 보여 줄 수 있음
  • 순서(p) 검색 시 변곡점이 예상되는데, 이는 일정 순서가 될 때까지 순서 p가 커지면 AIC점수가 감소하고, 이후 점수가 높아지기 시작한다

 

• grid-search를 수행해서 최적의 p
  • fit으로 VAR 모델을 학습
  • 1부터 10까지 적합한 순서에 대한 AIC 점수를 찾기 위해 반복문을 통해 grid-search

forecasting_model = VAR(train)
results_aic = []

for p in range(1,10):
  results = forecasting_model.fit(p)
  results_aic.append(results.aic)

• 결과 그래프에서 가장 낮은 AIC점수는 2이고, 그 이후 p가 커짐에 따라 증가 추세
  • VAR모델의 최적 순서는 2

sns.set()
plt.plot(list(np.arange(1,10,1)), results_aic)
plt.xlabel("Order")
plt.ylabel("AIC")
plt.show()

• 모형에 순서 2로 fit 시키고 요약 결과

results = forecasting_model.fit(2)
results.summary()

 

 

🔴 예측하기

• 학습된 모델에 2일 동안의 훈련을 넣어 향후 10일 동안의 테스트 데이터를 예측

laaged_values = train.values[-2:]
forecast = pd.DataFrame(results.forecast(y= laaged_values, steps=10), index = test.index, columns= ['realgdp_1d', 'realdpi_1d'])
forecast

 

• 언급한 예측이 차분(diffencing)에 대한 모델에 대한 것
• 차분을 더하여 우리가 예측해야 할 값으로 만들기
  • 왼쪽(_1d)은 차분에 대한 예측값
  • 오른쪽 (_forcasted)은 원래 시리즈에 대한 예측값

forecast["realgdp_forecasted"] = mydata["realgdp"].iloc[-10-1] + forecast['realgdp_1d'].cumsum()
forecast["realdpi_forecasted"] = mydata["realdpi"].iloc[-10-1] + forecast['realdpi_1d'].cumsum() 
forecast

 

• 실제 test 셋과 합쳐서 시각화
  • realdpi와 realdpi_forecasted는 비슷한 패턴
  • realgdp와 realgdp_forecasted는 절반 정도는 비슷하다가 다른 패턴

test = mydata.iloc[-10:,:]
test["realgdp_forecasted"] = forecast["realgdp_forecasted"]
test["realdpi_forecasted"] = forecast["realdpi_forecasted"]
test.plot()