[5] ARIMA - auto_arima를 활용한 예제 실습

 

 

 

Steps¶

1. 시계열 시각화 - 추세, 계절성, 주기 파악¶

1-1. 추세

: 추세가 있다면 차분이 필요, 차분 횟수를 알 수 있는 단위근 검정 실시 (pmdarima.arima의 ndiffs 함수로 구할 수 있음)

1-2. 계절성

: 계절성이 있다면, seasonal=True로 지정해주어 SARIMA로

1-3. 주기

: 주기가 있다면, auto_arima 함수의 m의 인자에 넣어주기

m=7 - daily m=52 - weekly (1년=52주) m=12 - monthly m=1 - 비계졀성 (default)

2. 모형 적합 - 적절한 p,d,q를 auto_arima로 추정¶

3. 잔차 검정 - 잔차가 정상성을 만족하느지 ,정규성 및 등분산성을 만족하는지 파악¶

4. 모형 refresh & 예측¶

• 한번에 테스트 데이터를 예측하는 것이 아니라, 하나씩 예측하고 관측치를 업데이트 해주어야함
• 가장 최신의 데이터들로 모형을 적합해야 미래의 값을 예측하기 좋을 것 #### 방법

1. 새 관측치가 왔을 때, 주기적으로 ARIMA 모형의 모수들을 재추정하기 (P,D,Q값을 바꿔주고, 계수도 다시 추정)
2. 아직 모수를 재추정할 때가 아니라면, 가장 최신의 데이터를 모형에 추가해서 새로운 예측치가 최신의 데이터를 반영할 수 있도록 해주기 → 1번은 모형을 새로 바꾸는 작업이기에 자주 하지 않음

→ 2번은 데이터가 쌓일 때마다 해주는 것이 좋음

5. 모형 평가 - MAPE¶

 

In [5]:

#pip install finance_datareader

 

In [8]:

#pip install pmdarima

 

In [10]:

import FinanceDataReader as fdr
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
from pmdarima.arima import ndiffs
import pmdarima as pm

 

In [11]:

df_krx = fdr.StockListing('KRX') # 한국거래소 상장종목 전체

 

In [13]:

df_krx.head()

 

Out[13]:

	Code	ISU_CD	Name	Market	Dept	Close	ChangeCode	Changes	ChagesRatio	Open	High	Low	Volume	Amount	Marcap	Stocks	MarketId
0	005930	KR7005930003	삼성전자	KOSPI		73000	1	400	0.55	72800	73000	72200	9861960	716958458106	435794126150000	5969782550	STK
1	373220	KR7373220003	LG에너지솔루션	KOSPI		428500	2	-2000	-0.46	434000	434000	427000	139380	59879431000	100269000000000	234000000	STK
2	000660	KR7000660001	SK하이닉스	KOSPI		128900	1	1400	1.10	128700	129100	127300	3696790	475006760500	93839504848500	728002365	STK
3	207940	KR7207940008	삼성바이오로직스	KOSPI		703000	2	-6000	-0.85	709000	710000	702000	48422	34105546000	50035322000000	71174000	STK
4	005935	KR7005931001	삼성전자우	KOSPI		58800	1	400	0.68	58700	58800	58400	2241939	131507020526	48385737960000	822886700	STK

 

In [19]:

def stockCode(val):
    df_krx = fdr.StockListing('KRX')
    df_name = df_krx[df_krx['Name'] == val]
    symbol = df_name.iloc[0,0]
    name = df_name.iloc[0,2]
    return symbol, name

 

In [21]:

stockCode('삼성전자')

 

Out[21]:

('005930', '삼성전자')

 

In [22]:

def stockMain(val, strFromDay, strToDay):
    print("stockMain")  
    symbol, name = stockCode(val)
    
    title = name + "("+symbol+")"
    
    df = fdr.DataReader(symbol, strFromDay, strToDay)
    
    return df

 

In [29]:

df = stockMain('삼성전자','2020-01-01','2021-08-30')

 

 

stockMain

 

In [30]:

df

 

Out[30]:

	Open	High	Low	Close	Volume	Change
Date
2020-01-02	55500	56000	55000	55200	12993228	-0.010753
2020-01-03	56000	56600	54900	55500	15422255	0.005435
2020-01-06	54900	55600	54600	55500	10278951	0.000000
2020-01-07	55700	56400	55600	55800	10009778	0.005405
2020-01-08	56200	57400	55900	56800	23501171	0.017921
...	...	...	...	...	...	...
2021-08-24	73900	75700	73900	75600	21016913	0.031378
2021-08-25	76200	76600	74900	75700	22319664	0.001323
2021-08-26	76100	76200	74600	74600	16671494	-0.014531
2021-08-27	74300	75000	73800	74300	15172748	-0.004021
2021-08-30	75400	75500	74200	74600	12686999	0.004038

413 rows × 6 columns

 

 

1. 시계열 자료 시각화¶

 

In [31]:

y_train = df['Close'][:int(0.7*len(df))]
y_test = df['Close'][int(0.7*len(df)):]
y_train.plot()
y_test.plot()

 

Out[31]:

<Axes: xlabel='Date'>

 

 

 

• 정상성을 만족하지 않아 차분이 필요함
• 상승 추세가 있으나 최근에는 주춤하는 모습
• 계절성, 주기성이 보이지 않음

 

In [32]:

kpss_diffs = ndiffs(y_train, alpha=0.05, test='kpss', max_d=6)
adf_diffs = ndiffs(y_train, alpha=0.05, test='adf', max_d=6)
n_diffs = max(adf_diffs, kpss_diffs)

print(f"추정된 차수 d = {n_diffs}")

 

 

추정된 차수 d = 1

 

 

2. 모형 적합¶

 

In [33]:

model = pm.auto_arima(y = y_train        
                      , d = 1          
                      , start_p = 0 
                      , max_p = 3   
                      , start_q = 0 
                      , max_q = 3   
                      , m = 1       
                      , seasonal = False 
                      , stepwise = True
                      , trace=True
                      )

 

 

Performing stepwise search to minimize aic
 ARIMA(0,1,0)(0,0,0)[0] intercept   : AIC=4934.821, Time=0.06 sec
 ARIMA(1,1,0)(0,0,0)[0] intercept   : AIC=4936.780, Time=0.02 sec
 ARIMA(0,1,1)(0,0,0)[0] intercept   : AIC=4936.799, Time=0.02 sec
 ARIMA(0,1,0)(0,0,0)[0]             : AIC=4934.428, Time=0.01 sec
 ARIMA(1,1,1)(0,0,0)[0] intercept   : AIC=4938.798, Time=0.05 sec

Best model:  ARIMA(0,1,0)(0,0,0)[0]          
Total fit time: 0.168 seconds

 

In [34]:

model = pm.auto_arima (y_train, d = 1, seasonal = False, trace = True)
model.fit(y_train)

 

 

Performing stepwise search to minimize aic
 ARIMA(2,1,2)(0,0,0)[0] intercept   : AIC=4942.230, Time=0.31 sec
 ARIMA(0,1,0)(0,0,0)[0] intercept   : AIC=4934.821, Time=0.01 sec
 ARIMA(1,1,0)(0,0,0)[0] intercept   : AIC=4936.780, Time=0.02 sec
 ARIMA(0,1,1)(0,0,0)[0] intercept   : AIC=4936.799, Time=0.02 sec
 ARIMA(0,1,0)(0,0,0)[0]             : AIC=4934.428, Time=0.01 sec
 ARIMA(1,1,1)(0,0,0)[0] intercept   : AIC=4938.798, Time=0.06 sec

Best model:  ARIMA(0,1,0)(0,0,0)[0]          
Total fit time: 0.435 seconds

Out[34]:

 ARIMA(0,1,0)(0,0,0)[0]          

In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

 

 

ARIMA(0,1,0)은 1차 차분을 했을 때 백색 잡음이 되는, 임의 보행 모형을 따르는 모형

 

 

3. 잔차 검정¶

 

In [35]:

print(model.summary())

 

 

                               SARIMAX Results                                
==============================================================================
Dep. Variable:                      y   No. Observations:                  289
Model:               SARIMAX(0, 1, 0)   Log Likelihood               -2466.214
Date:                Mon, 11 Dec 2023   AIC                           4934.428
Time:                        21:47:09   BIC                           4938.091
Sample:                             0   HQIC                          4935.896
                                - 289                                         
Covariance Type:                  opg                                         
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
sigma2      1.599e+06   9.63e+04     16.608      0.000    1.41e+06    1.79e+06
===================================================================================
Ljung-Box (L1) (Q):                   0.02   Jarque-Bera (JB):                46.72
Prob(Q):                              0.90   Prob(JB):                         0.00
Heteroskedasticity (H):               1.47   Skew:                             0.53
Prob(H) (two-sided):                  0.06   Kurtosis:                         4.67
===================================================================================

Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).

 

 

결과 해석 (잔차 검정)

• Ljung-Box(융-박스): 0.90으로 0.05를 넘어 귀무가설인 '잔차가 백색잡음을 따른다'를 기각할 수 없다. 시계열 모형이 잘 적합되었고, 낭믄 잔차는 더이상 자기상관을 갖지 않는 백색잡음
• Jarque-Bera(자크-베라): 0.00으로 0.05를 넘지 않아 귀무가설인 '잔차가 정규성을 만족한다'를 기각한다. 정규성을 따르지 않음
• Heteroskedasticity: 0.06으로 0.05를 넘어 귀무가설인 '잔차가 이분산을 띠지 않는다'를 기각한다. 이분산을 띠지 않음
• skew(비대칭도)는 0에 가까워야 하고, Kurtosis(첨도)는 3에 가까워야한다. 비대칭도는 0에 가깝지만, 첨도는 3보다 높은 값을 가짐

 

In [36]:

model.plot_diagnostics(figsize=(16, 8))
plt.show()

 

 

 

 

• Standardized Residual: 잔차를 시계열로 그린 것, 백색잡음이기에 평균 0을 중심으로 무작위하게 움직임 (=정상시계열)
• Correlogram: 잔차에 대한 ACF Plot, 자기상관을 띠지 않음
• Histogram, Normal Q-Q: 약간 벗어나는 모습
• 잔차가 백색잡음이지만, 정규성을 따르지 않음
  -> 모형 설정을 바꾸거나, P,D,Q의 MAX를 늘려보거나, train데이터를 더 최근인 데이터만 가져오거나, 분산 안정화를 시도해볼 수 있음

 

 

4. 모형 Refresh 및 예측¶

 

In [48]:

# 테스트 데이터 개수만큼 예측
y_predict = model.predict(n_periods=len(y_test)) 
y_predict = pd.DataFrame(y_predict,index = y_test.index,columns=['Prediction'])

 

In [49]:

# 그래프
fig, axes = plt.subplots(1, 1, figsize=(12, 4))
plt.plot(y_train, label='Train')        # 훈련 데이터
plt.plot(y_test, label='Test')          # 테스트 데이터
plt.plot(y_predict, label='Prediction')  # 예측 데이터
plt.legend()
plt.show()

 

 

 

 

모형이 상수항이 없는 임의보행 모형이기 때문에, 예측치들은 가장 마지막 관측치가 된다. 그래서 일직선을 얻게 되는 것

• auto_arima가 임의보행 모형을 적합한 모델로 선택한 이유는, 데이터에 특정한 주기나 추세가 없기 때문
• aic로 모형을 최적화하는 과정에서 의미있는 자기상관(AR)이나 이동평균(MA)를 찾기 어려웠기 대문에, 최선책으로 임의보행 모형(=어제의 값이 오늘의 값을 가장 잘 설명한다는 모형)을 선택한 것
• 즉, 데이터에서 특정한 구조를 보기 어려워 가장 마지막 관측치가 가장 좋은 예측치다라고 말하고 있다.
• 위처럼 예측되지 않게 하기 위해서는, 한 스텝씩 예측하고, 테스트 데이터를 관측할 때마다 모형을 업데이트해주는 REFRESH 과정이 필요

 

In [42]:

def forecast_one_step():
    fc, conf_int = model.predict(n_periods=1 # 한 스텝씩
        , return_conf_int=True)              # 신뢰구간 출력
    return (
        fc.tolist()[0],
        np.asarray(conf_int).tolist()[0]
    )
# 가장 마지막 관측치가 업데이트 되어 관측치가 추가됨으로써 모형에서의 예측치도 업데이트 됨

 

In [43]:

forecasts = []
y_pred = []
pred_upper = []
pred_lower = []

for new_ob in y_test:
    fc, conf = forecast_one_step()
    y_pred.append(fc)
    pred_upper.append(conf[1]) 
    pred_lower.append(conf[0])

    ## 모형 업데이트
    model.update(new_ob)

 

In [45]:

pd.DataFrame({"test": y_test, "pred": y_pred})

 

Out[45]:

	test	pred
Date
2021-03-05	82100	74600.0
2021-03-08	82000	82100.0
2021-03-09	81400	82000.0
2021-03-10	80900	81400.0
2021-03-11	82000	80900.0
...	...	...
2021-08-24	75600	73300.0
2021-08-25	75700	75600.0
2021-08-26	74600	75700.0
2021-08-27	74300	74600.0
2021-08-30	74600	74300.0

124 rows × 2 columns

 

 

오늘의 실제 값이 내일의 예측 값임을 확인할 수 있다.

 

In [47]:

print(model.summary())

 

 

                               SARIMAX Results                                
==============================================================================
Dep. Variable:                      y   No. Observations:                  537
Model:               SARIMAX(0, 1, 0)   Log Likelihood               -4529.559
Date:                Mon, 11 Dec 2023   AIC                           9061.118
Time:                        22:00:53   BIC                           9065.402
Sample:                             0   HQIC                          9062.794
                                - 537                                         
Covariance Type:                  opg                                         
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
sigma2      1.286e+06   4.18e+04     30.775      0.000     1.2e+06    1.37e+06
===================================================================================
Ljung-Box (L1) (Q):                   1.34   Jarque-Bera (JB):               651.12
Prob(Q):                              0.25   Prob(JB):                         0.00
Heteroskedasticity (H):               0.77   Skew:                             1.02
Prob(H) (two-sided):                  0.08   Kurtosis:                         8.00
===================================================================================

Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).

 

In [50]:

from plotly.subplots import make_subplots
import plotly.graph_objects as go

fig = go.Figure([
    # 훈련 데이터
    go.Scatter(x = y_train.index, y = y_train, name = "Train", mode = 'lines'
              ,line=dict(color = 'royalblue'))
    # 테스트 데이터
    , go.Scatter(x = y_test.index, y = y_test, name = "Test", mode = 'lines'
                ,line = dict(color = 'rgba(0,0,30,0.5)'))
    # 예측값
    , go.Scatter(x = y_test.index, y = y_pred, name = "Prediction", mode = 'lines'
                     ,line = dict(color = 'red', dash = 'dot', width=3))
    
    # 신뢰 구간
    , go.Scatter(x = y_test.index.tolist() + y_test.index[::-1].tolist() 
                ,y = pred_upper + pred_lower[::-1] ## 상위 신뢰 구간 -> 하위 신뢰 구간 역순으로
                ,fill='toself'
                ,fillcolor='rgba(0,0,30,0.1)'
                ,line=dict(color='rgba(0,0,0,0)')
                ,hoverinfo="skip"
                ,showlegend=False)
])

fig.update_layout(height=400, width=1000, title_text="ARIMA(0,1,0)모형")
fig.show()

 

 

 

 

 

 

5. 모형 평가¶

 

In [51]:

def MAPE(y_test, y_pred):
	return np.mean(np.abs((y_test - y_pred) / y_test)) * 100 
    
print(f"MAPE: {MAPE(y_test, y_pred):.3f}")

 

 

MAPE: 0.862