시계열(Time series) > Forecasting
2021. 1. 2. 10:54ㆍ계량경제학
반응형
지금까지 다루어 왔던 ARIMA, ARIMAX, SARIMA, SARIMAX를 활용하여, 예측 문제를 해결해보도록 하겠습니다.
먼저 ARIMAX와 SARIMAX를 추정하기 위해선, Univariate 시계열 데이터뿐만 아니라 추가적인 Exogenous 변수가 필요하게 됩니다. 따라서 이번 포스팅에서 활용한 데이터는 기존의 S&P 500 index와 더불어 Nikkei index를 활용하고자 합니다.(이전 포스팅에서와 동일하게 yahoo finance- historical data에서 다운로드 받았습니다.)
최종적으로는, Nikkei index의 return을 예측해보도록 해보겠습니다.
먼저 필요한 Module을 import하고, 함수를 정의하도록 하겠습니다.
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import scipy.stats
import itertools
import warnings
import os
import statsmodels.api as sm
import statsmodels.graphics.tsaplots as sgt
import statsmodels.tsa.stattools as sts
from statsmodels.tsa.seasonal import seasonal_decompose
from statsmodels.tsa.arima_model import ARIMA
from scipy.stats.distributions import chi2
org_path = "../Time series"
os.chdir(org_path)
sns.set()
warnings.filterwarnings("ignore")
## LLR test
def LLR_test(mod_1, mod_2, DF=1):
L1 = mod_1.llf
L2 = mod_2.llf
LR = (2*(L2-L1))
p = chi2.sf(LR, DF).round(3)
return p
데이터를 로딩하고, S&P 500과 Nikkei 데이터의 Length가 다르기 때문에, Length를 같게 맞춰준 후 하나의 데이터프레임으로 합칩니다. 그 다음에, " return_snp", "return_nikkei" 칼럼을 만들게 되는데 reutrn은 아래와 같은 식으로 계산됩니다.
nikkei.Date = pd.to_datetime(nikkei.Date, dayfirst = False)
## Setting the start date and end date
snp = snp.loc[snp["Date"] != "1994-01-03"]
nikkei=nikkei.loc[nikkei["Date"] < "2020-01-01"]
## Setting index to make two dataframe with same length
snp.set_index("Date", inplace = True)
snp=snp.asfreq(freq="B")
nikkei.set_index("Date", inplace = True)
nikkei=nikkei.asfreq(freq="B")
print("snp:", len(snp),"nikkei:", len(nikkei))
## Merge
df = pd.merge(snp, nikkei, on="Date")
# df.set_index("Date", inplace = True)
df=df.fillna(method='ffill')
df = df[["Close_x", "Close_y"]]
df.columns = ["SNP", "Nikkei"]
## diff
df["diff_SNP"]=df["SNP"].diff()
df["diff_Nikkei"]=df["Nikkei"].diff()
## return
df["return_snp"]=df["SNP"].pct_change(1)*100
df["return_nikkei"]=df["Nikkei"].pct_change(1)*100
## Train, test
size = int(len(df)*0.8)
df_train = df.iloc[:size]
df_test = df.iloc[size:]Correlogram을 활용하여 우리가 가지고 있는 시계열 데이터의 특성을 살펴 보도록 해보겠습니다. Retrun 데이터의 경우 Stationary하다고 볼 수 있으며, 살펴보았을 때 MA, AR모드 Lag 6까지 살펴볼 필요가 있을 것으로 보입니다.
##Decompose
decomposition =seasonal_decompose(df_train["return_nikkei"][1:], period =1)
fig = decomposition.plot()
fig.set_size_inches(10,10)
plt.show()
## ACF, PACF
fig = plt.figure()
ax1 = fig.add_subplot(2, 1, 1)
ax2 = fig.add_subplot(2, 1, 2)
sgt.plot_acf(df_train["return_nikkei"][1:], lags = 20, zero = False, ax=ax1)
ax1.set_title("ACF Nikkei") ##MA6
sgt.plot_pacf(df_train["return_nikkei"][1:], lags = 20, zero = False, method = ('ols'), ax=ax2)
ax2.set_title("PACF Nikkei") ##AR6
fig.set_size_inches(10,10)
plt.show()
## ADF
sts.adfuller(df_train["return_nikkei"][1:]) ##p-value : 0.0
이제 AR과 MA 시차를 6까지 하여, Gridsearch를 하도록 해보겠습니다. Likelihood와 AIC를 기준으로 검색해보겠습니다.
## Selecting lag
p = range(0,7)
d = range(0,1)
q = range(0,7)
pdq = list(itertools.product(p,d,q))
dict_model ={}
for i in pdq:
try:
model = ARIMA(df["return_nikkei"][1:], order=(i))
print(i)
model_fit = model.fit( maxiter = 100)
dict_model[i]=[model_fit.llf, model_fit.aic]
except:
print("error lag:",i)
information=pd.DataFrame.from_dict(dict_model, orient ="index", columns =["llf", "AIC"])
information.loc[information["llf"] == information["llf"].max()]
'''
llf AIC
(6, 0, 6) -11985.306923 23998.613847
'''
information.loc[information["AIC"] == information["AIC"].min()]
'''
llf AIC
(6, 0, 4) -11986.521433 23997.042866
'''
information.sort_values(by=["AIC"], ascending = True)
'''
llf AIC
(6, 0, 4) -11986.521433 23997.042866
(6, 0, 6) -11985.306923 23998.613847
(1, 0, 1) -11996.113118 24000.226236
(4, 0, 3) -11991.984544 24001.969088
(2, 0, 1) -11996.110825 24002.221650
(0, 0, 2) -11997.178110 24002.356220
(2, 0, 2) -11995.267446 24002.534892
(6, 0, 3) -11990.375607 24002.751214
(5, 0, 5) -11989.407329 24002.814659
(2, 0, 0) -11997.432949 24002.865898
...
'''
##(6,0,4)
model_ar_6_ma_4 = ARIMA(df_train["return_nikkei"][1:], order=(6,0,4))
results_ar_6_ma_4 = model_ar_6_ma_4.fit()
results_ar_6_ma_4.summary()
'''
=======================================================================================
coef std err z P>|z| [0.025 0.975]
---------------------------------------------------------------------------------------
const 0.0076 0.018 0.414 0.679 -0.028 0.043
ar.L1.return_nikkei -0.0827 0.027 -3.063 0.002 -0.136 -0.030
ar.L2.return_nikkei 1.3293 0.028 46.777 0.000 1.274 1.385
ar.L3.return_nikkei -0.0702 0.028 -2.527 0.012 -0.125 -0.016
ar.L4.return_nikkei -0.8979 0.028 -31.801 0.000 -0.953 -0.843
ar.L5.return_nikkei -0.0307 0.014 -2.176 0.030 -0.058 -0.003
ar.L6.return_nikkei -0.0488 0.014 -3.453 0.001 -0.077 -0.021
ma.L1.return_nikkei 0.0482 0.024 2.047 0.041 0.002 0.094
ma.L2.return_nikkei -1.3711 0.026 -53.278 0.000 -1.422 -1.321
ma.L3.return_nikkei 0.1093 0.021 5.203 0.000 0.068 0.150
ma.L4.return_nikkei 0.9488 0.021 45.092 0.000 0.908 0.990
'''
##(6,0,6)
model_ar_6_ma_6 = ARIMA(df_train["return_nikkei"][1:], order=(6,0,6))
results_ar_6_ma_6 = model_ar_6_ma_6.fit()
results_ar_6_ma_6.summary()
'''
=======================================================================================
coef std err z P>|z| [0.025 0.975]
---------------------------------------------------------------------------------------
const 0.0076 0.018 0.418 0.676 -0.028 0.043
ar.L1.return_nikkei 0.2479 0.312 0.795 0.427 -0.363 0.859
ar.L2.return_nikkei 1.2646 0.275 4.593 0.000 0.725 1.804
ar.L3.return_nikkei -0.5410 0.447 -1.210 0.226 -1.418 0.336
ar.L4.return_nikkei -0.7390 0.409 -1.806 0.071 -1.541 0.063
ar.L5.return_nikkei 0.2834 0.294 0.964 0.335 -0.293 0.860
ar.L6.return_nikkei -0.1365 0.274 -0.498 0.618 -0.673 0.400
ma.L1.return_nikkei -0.2828 0.313 -0.904 0.366 -0.896 0.330
ma.L2.return_nikkei -1.2958 0.279 -4.640 0.000 -1.843 -0.748
ma.L3.return_nikkei 0.5804 0.449 1.293 0.196 -0.299 1.460
ma.L4.return_nikkei 0.7755 0.416 1.864 0.062 -0.040 1.591
ma.L5.return_nikkei -0.3099 0.291 -1.066 0.286 -0.880 0.260
ma.L6.return_nikkei 0.0962 0.275 0.350 0.726 -0.442 0.635
'''
LLR_test(results_ar_6_ma_4,results_ar_6_ma_6, DF=2)
Likelihood관점에선 ARIMA(6,6)이 AIC관점에서는 ARIMA(6,4)이 가장 좋은 지표를 나타냈습니다. ARIMA(6,6)과 ARIMA(6,4) LLR test를 통해 유의하게 Likelihood가 다른지를 검증해본 결과, ARIMA(6,6)이 ARIMA(6,4) 유의하게 Likelihood가 높다고 보기 어렵기 때문에 최종적으로 ARIMA(6,4)로 선택하도록 하겠습니다.
또한 추가적으로 추정된 모델의 Parameter의 P-value를 봤을 때 ARIMA(6,4)일 때 계수가 대부분 유의하였습니다. 그렇기 때문에 ARIMA(6,4)가 더 적절한 모델이라고 판단할 수 있을 것입니다.
이제 차수 (6,0,4)를 기본으로 하여 ARIMA, ARIMAX, SARIMA, SARIMAX 모형을 수립하여 예측해보도록 하겠습니다.
## ARIMA
model_ar_6_ma_4 = ARIMA(df_train["return_nikkei"][1:], order=(6,0,4))
results_ar_6_ma_4 = model_ar_6_ma_4.fit()
results_ar_6_ma_4.summary()
## ARIMAX(6,4)
model_ar_6_ma_4_X = ARIMA(df_train["return_nikkei"][1:], exog = df_train["return_snp"][1:],order=(6,0,4))
results_ar_6_ma_4_X = model_ar_6_ma_4_X.fit()
results_ar_6_ma_4_X.summary()
'''
=======================================================================================
coef std err z P>|z| [0.025 0.975]
---------------------------------------------------------------------------------------
const 0.0015 0.016 0.089 0.929 -0.031 0.034
return_snp 0.1855 0.019 9.889 0.000 0.149 0.222
ar.L1.return_nikkei -1.1859 nan nan nan nan nan
ar.L2.return_nikkei -0.0211 nan nan nan nan nan
ar.L3.return_nikkei 0.8509 nan nan nan nan nan
ar.L4.return_nikkei 0.2901 0.280 1.035 0.301 -0.259 0.839
ar.L5.return_nikkei 0.0655 nan nan nan nan nan
ar.L6.return_nikkei 0.0278 nan nan nan nan nan
ma.L1.return_nikkei 1.0876 nan nan nan nan nan
ma.L2.return_nikkei -0.1330 nan nan nan nan nan
ma.L3.return_nikkei -0.9187 nan nan nan nan nan
ma.L4.return_nikkei -0.2262 nan nan nan nan nan
'''
## ARIMAX(1,1)
model_ar_1_ma_1_X = ARIMA(df_train["return_nikkei"][1:], exog = df_train["return_snp"][1:],order=(1,0,1))
results_ar_1_ma_1_X = model_ar_1_ma_1_X.fit()
results_ar_1_ma_1_X.summary()
## SARIMA & SARIMAX
from statsmodels.tsa.statespace.sarimax import SARIMAX
model_sarima = SARIMAX(df_train["return_nikkei"][1:], order=(6,0,4), seasonal_order = (1,0,1,10),initialization='approximate_diffuse')
results_sarima = model_sarima.fit(maxiter=50)
results_sarima.summary()
model_sarimax = SARIMAX(df_train["return_nikkei"][1:], exog = df_train["return_snp"][1:], order=(6,0,4), seasonal_order = (1,0,1,10),initialization='approximate_diffuse')
results_sarimax = model_sarimax.fit(maxiter=50)
results_sarimax.summary()여기서 ARIMAX(6,4)의 경우에는 Summary를 보면 std err가 nan으로 표기되는 것을 볼 수 있습니다. 이는 모형의 과적합을 뜻할 수 있습니다. 따라서 차수를 임의로 (1,1)로 나누어 추정하였습니다.
Seaonality는 임의로 10으로 정의하를 하여 추정하였습니다. 사실 우리가 가지고 있는 데이터를 살펴보면 주기성을 가진 Pattern이 나타난다고 보기는 어려울 것 같지만 실습하는 의미로 SARIMA와 SARIMAX를 추정해보도록 하겠습니다.
이제 Train date set을 활용하여 Fitting을 해보고 RMSE기준으로 어떤 모형이 Train data를 잘 설명하는지를 살펴보도록 하겠습니다.
## RMSE
def rmse_(real, pred):
diff = real-pred
diff2=np.sqrt(diff**2)
return diff2.mean()
## Fitting
start_date = '1994-01-05'
end_date = '2014-10-16'
df_train["prediction_arima"] = results_ar_6_ma_4.predict(start = start_date, end = end_date)
df_train["prediction_arimax"] = results_ar_1_ma_1_X.predict(start = start_date, end = end_date, exog = df_train["return_snp"])
df_train["prediction_sarima"] = results_sarima.predict(start = start_date, end = end_date)
df_train["prediction_sarimax"] = results_sarimax.predict(start = start_date, end = end_date , exog = df_train["return_snp"])
fig= plt.figure(figsize = (20,20))
plt.plot(df_train["return_nikkei"][-100:], "b-", label = "Actual", alpha = 0.5, lw=5)
plt.plot(df_train["prediction_arima"][-100:] , "r-", label = "ARIMA" , alpha = 0.5, lw=5)
plt.plot(df_train["prediction_arimax"][-100:] , "g-" , label = "ARIMAX" , alpha = 0.5, lw=5)
plt.plot(df_train["prediction_sarima"][-100:], "y-" , label = "SARIMA" , alpha = 0.5, lw=5)
plt.plot(df_train["prediction_sarimax"][-100:] , "k-" , label = "SARIMAX" , alpha = 0.5, lw=5)
plt.legend()
plt.show()
rmse_(df_train["return_nikkei"],df_train["prediction_arima"]) ##1.0304148654167324
rmse_(df_train["return_nikkei"],df_train["prediction_arimax"]) ##1.0283191038651487
rmse_(df_train["return_nikkei"],df_train["prediction_sarima"]) ##1.0319786474733221
rmse_(df_train["return_nikkei"],df_train["prediction_sarimax"]) ##1.02870513544557
RMSE 관점에서는 ARIMAX(1,1)이 가장 좋은 Fitting을 보여주었습니다.
이제 Test dataset를 바탕으로 예측력을 평가해 보겠습니다.
## Forecasting
steps = len(df_test)
df_test["prediction_arima"] = results_ar_6_ma_4.forecast(steps)[0]
df_test["prediction_arimax"] = results_ar_6_ma_4_X.forecast(steps, exog = df_test["return_snp"])[0]
df_test["prediction_sarima"] = results_sarima.forecast(steps).values
df_test["prediction_sarimax"] = results_sarimax.forecast(steps, exog = df_test[["return_snp"]]).values
fig= plt.figure(figsize = (20,20))
plt.plot(df_test["return_nikkei"][-100:], "b-", label = "Actual", alpha = 0.9, lw=3)
plt.plot(df_test["prediction_arima"][-100:] , "r-", label = "ARIMA" , alpha = 0.9, lw=3)
plt.plot(df_test["prediction_arimax"][-100:] , "g-" , label = "ARIMAX" , alpha = 0.9, lw=3)
plt.plot(df_test["prediction_sarima"][-100:] , "y-" , label = "SARIMA" , alpha = 0.9, lw=3)
plt.plot(df_test["prediction_sarimax"][-100:] , "k-" , label = "SARIMAX" , alpha = 0.5, lw=3)
plt.legend()
plt.show()
rmse_(df_test["return_nikkei"],df_test["prediction_arima"]) ##0.7961663716793427
rmse_(df_test["return_nikkei"],df_test["prediction_arimax"]) ##0.7882616502694494
rmse_(df_test["return_nikkei"],df_test["prediction_sarima"]) ##0.7957024431703473
rmse_(df_test["return_nikkei"],df_test["prediction_sarimax"]) ##0.787914126856654
Test date set기준으론 SARIMAX가 가장 RMSE가 낮았습니다. Univariate 시계열 모형만 사용하기 보다는 Exogenous 변수를 활용할 때 더 좋은 예측력을 가지고 있는 것 같습니다.
마지막으로 pmdarima 패키지를 활용하여 자동으로 시차를 결정해보도록 하겠습니다. 해당 패키지에 대한 자세한 내용은 alkaline-ml.com/pmdarima/modules/generated/pmdarima.arima.auto_arima.html를 참조하시면 좋을 것 같습니다.
## auto-arima
from pmdarima.arima import auto_arima
model_auto = auto_arima(df_train["return_nikkei"][1:], exogenous = df_train[['return_snp']][1:], m = 10,
max_order = None, max_p = 6, max_q = 6, max_P = 4, max_Q = 4,
maxiter = 50, alpha = 0.05, n_jobs = -1, with_intercept=True,
information_criterion = 'oob', out_of_sample = int(len(df_train)*0.2))
model_auto.summary()
'''
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
return_snp 0.1857 0.012 16.126 0.000 0.163 0.208
ar.L1 0.3776 0.059 6.358 0.000 0.261 0.494
ma.L1 -0.4775 0.056 -8.598 0.000 -0.586 -0.369
sigma2 2.1245 0.021 102.071 0.000 2.084 2.165
'''
df_test["auto_arima"] =model_auto.predict(n_periods = steps, exogenous=df_test[["return_snp"]])
rmse_(df_test["return_nikkei"],df_test["auto_arima"])Intercept가 없는 ARIMAX(1,1)을 최적의 모델로 찾아주었습니다. 여기서 Criterion으로 활용한것은 Validation Set을 활용한 MSE를 기준으로 모델을 탐색하였습니다.
이번 포스팅은 사실 인과관계를 분석하고 통계적 엄밀성을 추구하는 학문적 관점에서의 계량경제학과는 조금 거리가 있지만 실무적으로 시계열분석 방법론을 적용하는 방안에 대해서 정리해보았습니다.
반응형
'계량경제학' 카테고리의 다른 글
| 시계열(Time series) > Break (0) | 2021.01.05 |
|---|---|
| 시계열(Time series) > ARIMAX, SARIMA, SARIMAX (2) | 2021.01.01 |
| 시계열(Time series) > Diagnosing Models(ARMA, ARIMA)(2/2) (0) | 2020.12.27 |
| 시계열(Time series) > Diagnosing Models(ARMA, ARIMA)(1/2) (0) | 2020.12.27 |
| 시계열(Time series) > Moving average model(이동평균모형) (0) | 2020.12.23 |