카테고리 없음

[시계열] 시계열 모델 기초 (2) 불안정 시계열 - ARIMA (세상에서 제일 쉬운 설명)

허니비 honeybee 2023. 3. 16. 16:32

우리가 지난 포스팅에서 살펴본 AR, MA, ARMA 모델은 안정 시계열이었으나, 대부분의 시계열 모델은 정상성 등을 파괴해버리는 특성을 가지는 불안정 시계열인 경우가 크다. 따라서 차분이나 로그 변환을 통해 해결할 수 있는 불안정 시계열 모델을 이용하는게 일반적이다. 

 

게 중 ARIMA 모델은 가장 기본적이고도 유명한 시계열 모델인데,, 오늘은 이 ARIMA 모델의 내부를 들여다 보고 간단한 실습까지 해보도록 하겠다. 앞선 포스팅을 보지 않았다면 꼭 보고 와주시길! 

 

 

 

 

 

1.  ARIMA 모델 설명 

 

 ARIMA(p, d, q) : d차 차분한 데이터에 AR(p) 모형과 MA(q) 모형을 합친 모형

 

ARIMA 모델은 이름에서도 알 수 있듯 AR 모델의 파라미터 p 와  MA 모델의 파라미터 q 를 Inetegrate 한 모델이다. d 만큼 차분을 해서 나타낸 것. 사실 시계열에서 AR 과 MA 를 봤다는 뜻은 자기회귀와 오차 두가지를 본 것이다. 그런데 저 두 가지만 이용하면 차분이라는 중요한 과정이 빠져 불안정 시계열 데이터에서는 쓸 수 없다. 따라서 두 가지를 보면서 차분을 추가해서 합친 것이 바로 ARIMA 라고 이해하면 되겠다. 

 

요약: ARIMA 는 비정상 시계열을 정상 시계열로 만들기 위해 합치는 과정에서 차분을 추가한 모형 

 

 

 

 

2.  ARIMA 의 파라미터 

 

ARIMA(p, 0, 0) = AR(p)
ARIMA(0, 0, q) = MA(q)

 

아리마의 파라미터는 p, d, q 세 가지가 있다. AR 의 파라미터인 p, 차분의 d, MA 의 파라미터인 q 가 있는데 따라서 ARIMA (p, 0, 0) 은 MA 는 고려하지 않고 AR만 p 만큼 고려한 모델이며, ARIMA(0, 0, q)  역시 AR 은 고려하지 않고 MA 만 (q) 만큼 고려한 모델이라고 볼 수 있다. 

 

 

이미 최근에는 이 파라미터도 다 자동화되고 파이썬 패키지가 알아서 계산을 해주기 때문에 굳이 하나하나 테스트하며 돌릴 필요는 없지만... 파라미터들이 의미하는 것을 알고 있는 것은 중요하니까! 이 부분도 언급해 보았다 

 

 

 

 

 

3. 파라미터 튜닝 

 

p, d, q 값은 AIC, BIC, 그리고 경험적인 자기상관 등 다양한 방법을 통해서 선택할 수 있다
자기상관은 ACF 와 PACF 그래프를 그려서 확인한다

시계열 데이터가 AR의 특성을 띄는 경우, ACF는 천천히 감소하고 PACF 는 처음 시차를 제외하고 급격히 감소

 

-> 나(현재 데이터) 와의 상관관계를 띄기 때문에 

 

시계열 데이터가 MA의 특성을 띄는 경우, ACF는 급격히 감소하고 PACF는 천천히 감소

 

-> 백색잡음 (noise) 와의 관계성을 띄기 때문에 

 

 

 

ACF 와 PACF 가 무엇을 의미하는지 역시 따로 포스팅을 파서 설명해 두었으니 참고 바란다

 

 

 

 

 

 

 

 

 

4. ARIMA 파이썬 실습 

 

4.1.  전처리 + EDA + 목표 확인 

 

이번 실습에서 사용하게 될 데이터는 케글의 날씨 데이터, 위치 데이터, 공습 정보 데이터를 짜집기한 데이터이다. 세 데이터를 자르고 시간 순대로 이어붙인 후, EDA 를 통해 알 수 있었던것들은 다음과 같은데, 그냥... 기본적인 전처리이고 다른건 다 value counts 찍어보고 멧플로립 그래프 만들어본거니 ㅋㅋㅋ 코드나 자세한 설명은 생략하겠다. 

 

 

- (목표 수립) 날씨 데이터를 통해 이상 기후를 파악하고 싶다 

- (데이터 편중 문제) 공습 지역은 미국이 제일 , 영국이 2등으로 많았음 

- (데이터 편중 문제) Target country 역시 지역별로 편중이 심했다. 이탈리아, 리비아 등등이 많았음 

- (불안정 시계열 문제) 날씨 데이터는 너무나 당연히 주기성을 띄는 데이터였음 (주기성이 크고, 진폭이 큼, 상한선 하안선은 일정한 편) 

 

 

 

4.2 시계열 데이터의 정상성 검정: ADF 

 

 

 

그래 저 그래프는 한눈에 봐도 정상성을 만족하지 않는 모습을 보이나... ㅋㅋㅋ 예의상 ADF (Dicky Fuller test) 를 이용해서 정상성 검정을 해보도록 하자. ADF test 는 계속 언급했지만 정상성이 있는지 보는 test 다. 파이썬 패키지로 구현되어 있다. 

 

귀무가설: 시계열 자료가 정상 시계열이 아니다. (비정상 시계열이다)
대립가설: 시계열 자료가 정상 시계열이다  (reject 시 채택할 가설) 
 
 
 
- ADF 테스트 돌려보기 
# adfuller library 
from statsmodels.tsa.stattools import adfuller

# Check adfuller
def check_adfuller(ts):
    # Dickey-Fuller test
    result = adfuller(ts, autolag='AIC')
    print('Test statistic: ' , result[0])
    print('p-value: '  ,result[1])
    print('Critical Values:' ,result[4])
    
#out 
Test statistic:  -1.4095966745887756
p-value:  0.5776668028526357
Critical Values: {'1%': -3.439229783394421, '5%': -2.86545894814762, '10%': -2.5688568756191392}

 

 

ADF 테스트를 볼 때는 Test statistic 이 Critical Value 내에 들어가는지 보면 된다. p-value 는 Critical value 안에 들어갈 확률인데, 지금은 0.57 수준으로 절반 이상의 확률을 보이기 때문에 유의수준을 기각할 확률 (일반적으로 5프로나 1프로) 내에 들지 않아 귀무가설을 채택하게 된다. 따라서 이 시계열은 불안정 시계열이다 

 

 

* p-value 는 귀무가설이 맞을 때, 검정통계량과 같거나 더 극단적인 값을 얻을 확률을 나타내는 값인데, p-value 가 작을수록 귀무가설을 기각할 증거가 강해진다. 주로 유의수준 0.05나 0.01 하에서 검정을 진행하는데, p-value 가 유의수준보다 작다면 앞서 세운 귀무가설을 기각하고 대립가설을 채택하게 된다. 

 

 

* test statistic 은 특정 가설에 대해 검정할 때, 측정된 표본과 가설에서 예상되는 값의 차이? 라고 한다. 검정통계량이 커질수록 귀무가설을 기각할 확률이 커지는데, 여기서는 값이 음수이기 때문에 측정된 표본통계량이 가설에서 예상되는 값보다 작다는걸 의미한다. 그러나 판정은 p-value 로 내리는 것이 일반적인듯? 

 

 

 

 

 

- rolling mean 과 rolling std 그래프 그려보기 

# Check mean & std
def check_mean_std(ts):
    # Rolling statistics
    rolmean = ts.rolling(6).mean()
    rolstd = ts.rolling(6).std()

    plt.figure(figsize=(22,10))   
    orig = plt.plot(ts, color='red',label='Original')
    mean = plt.plot(rolmean, color='black', label='Rolling Mean')
    std = plt.plot(rolstd, color='green', label = 'Rolling Std')
    plt.xlabel("Date")
    plt.ylabel("Mean Temperature")
    plt.title('Rolling Mean & Standard Deviation')
    plt.legend()
    plt.show()

 

 

이동평균과 이동 표준편차를 그려보았을 때 이런 모양새가 나왔다. 표준편차가 확 오르는 구간 (1945-04 부근) 에서 변동성이 심하고 이상치가 될 법 하다? 정도의 정보를 얻고 갈 수 있겠다. 원본 데이터를 보니 (빨간선) 확 기온이 꺾이는 이상함이 보이는 것.. 맞다 

 

 

 

* Rolling mean 과 Rolling std 는 시계열 데이터에서 이동평균과 이동표준편차를 의미한다. 이 둘이 우리가 흔히 알고 있는 그냥 평균, 그냥 표준편차와의 차이는 주어진 시계열데이터에서 일정한 기간 (주로 window 라고 부름) 만큼의 평균과 표준편차를 계산한다는거다. data.rolling(6).mean -> 6일간의 이동평균을 구해줘 라고 이해하면 된다. 

 

굳이 왜 이걸 쓰냐면.. 시계열 데이터에서 나타나는 추세와 변동성을 더 쉽게 알고 파악할 수 있기 때문이다. 기간을 정해서 평균과 표편을 구하는게 더 낮은 변동성을 가지게 되니까! 일시적인 변동에 대한 영향력을 줄여주니깐!! 

 

 

4.3. 차분을 통해 정상시계열 만들기 

 

 

window size 는 6일간의 데이터를 기준치로 보겠다는 것이다. 그렇게 이동평균인 moving avg 를 구했을 때 원래의 시계열 데이터보다 훨씬 smoothing 된 형태로 oulier 가 약간 제거된 형태의 시계열을 볼 수 있다 

# Moving Average method
window_size = 6
moving_avg = ts.rolling(window_size).mean()

plt.figure(figsize=(22,10))
plt.plot(ts, color = "red",label = "Original")
plt.plot(moving_avg, color='black', label = "moving_avg_mean")
plt.title("Mean Temperature of Bindukuri Area")
plt.xlabel("Date")
plt.ylabel("Mean Temperature")
plt.legend()
plt.show()

 

 

이번에는 데이터를 차분하고 그래프를 찍어보겠다 

 

데이터를 차분하는 과정은 현재값에서 직전값을 빼는 것이다. 여기서는 window size 를 6으로 둔 상태의 이동평균을 원본값에서 빼는 방법으로 차분을 진행했다. 

#원본데이터 - 이동평균 
ts_moving_avg_diff = ts - moving_avg

# 처음 6개는 window size 때문에 nan 임으로 drop 처리 
ts_moving_avg_diff.dropna(inplace=True) 

# 정상성 검증: mean, variance(std)and adfuller test
check_mean_std(ts_moving_avg_diff)
check_adfuller(ts_moving_avg_diff.MeanTemp)


# differencing method2 직전값 빼기
ts_diff = ts - ts.shift()
plt.figure(figsize=(22,10))
plt.plot(ts_diff)
plt.title("Differencing method") 
plt.xlabel("Date")
plt.ylabel("Differencing Mean Temperature")
plt.show()

Test statistic: -11.1385143351385

p-value: 3.150868563164087e-20

Critical Values: {'1%': -3.4392539652094154, '5%': -2.86546960465041, '10%': -2.5688625527782327}

 

측정된 표본통계량이(Test Stat) 가설에서 예상되는 값보다(Critical Values) 보다 크고 (그리고 굉장한 음의 값을 띄고) , 이는 결국 귀무가설을 반박할 강한 근거가 된다 (대립가설을 채택하게 된다 -> 정상 시계열이다) 또한 p-value 역시 아주 작은 값을 띄고 있기 때문에 (0.05 이하) 귀무가설을 기각해도 된다. 

 

뭐.. 이런 것 없이 딱 눈으로 봐도 정상성을 띄는 시계열인게 보이지만 말이다 

 

 

* 차분을 하는 방법도 두가지인데, 하나는 moving avg 를 빼는 방법. 주로 이 방법은 트렌드나 천천히 움직이는 요소들을 제거할 때 쓰인다. 두번째는 직전값을 빼는 법. 직전값을 빼는 방법은 계절성을 제거할 때 쓰이고, 단기적인 변화를 제거할 때 유용하다. Additionally, subtracting a moving average can result in smoother data than taking the first difference, which can be useful in some applications. However, subtracting the value before the value can provide more information about the rate of change in the data and can be more useful in identifying specific points of inflection or turning points in the time series.

 

 

 

 

 

이제 이때의 데이터는 완전한 random noise 의 값을 보이는 ! 완전한 백색잡음의 데이터라고 볼 수 있겠다. 

 

 

 

 

 

 

 

 

 

4.4. ARIMA 로 평균기온 예측하기 

 

먼저 ARIMA 를 이용하기 전에 차분을 해서 변형한 데이터의 정상성을 알아보기 위해!  ACF 와 PACF 그래프를 statsmodel 패키지를 통해 그려보자. lag 를 20을 주고 그렸다

 

*기본적으로 ACF 는 자기상관, t 시점 기준으로 t-1 시점, t-2 시점... t-100 시점 등 여러 다른 지점의 관측값을 보겠다는 검정 ACF 는 정상 시계열인 경우 0으로 빠르게 떨어지고, 비정상 시계열은 천천히 떨어진다 (0차일 경우는 당연히 1이고 1차부터 확 떨어지는 모양)

 

*PACF 는 부분자기상관, 산정하려는 딱 두 지점의 상관관계를 보는 방법. 두 지점 사이의 다른 시점의 영향력은 제거. 정상 시계열인 경우 신뢰한계의 범위 내에 존재함. 

 

# ACF and PACF 
from statsmodels.tsa.stattools import acf, pacf
lag_acf = acf(ts_diff, nlags=20)
lag_pacf = pacf(ts_diff, nlags=20, method='ols')
# ACF
plt.figure(figsize=(16,8))

plt.subplot(121) 
plt.plot(lag_acf)
plt.axhline(y=0,linestyle='--',color='gray')
plt.axhline(y=-1.96/np.sqrt(len(ts_diff)),linestyle='--',color='gray')
plt.axhline(y=1.96/np.sqrt(len(ts_diff)),linestyle='--',color='gray')
plt.title('Autocorrelation Function')

# PACF
plt.subplot(122)
plt.plot(lag_pacf)
plt.axhline(y=0,linestyle='--',color='gray')
plt.axhline(y=-1.96/np.sqrt(len(ts_diff)),linestyle='--',color='gray')
plt.axhline(y=1.96/np.sqrt(len(ts_diff)),linestyle='--',color='gray')
plt.title('Partial Autocorrelation Function')
plt.tight_layout()

 

 

여기서는 ACF, PACF 그래프의 신뢰한계가 점선으로 표현되는데, 차분을 한 뒤의 데이터를  돌려봤더니 1차 이후에 값이 뚝 떨어지는 걸 볼 수 있다. 여기서 ACF 는 1차에서 값이 제일 낮고, PACF 는 2차에서 값이 제일 낮은 점을 고려해, MA 의 차수는 2, AR 의 차수는 1로 정하는게 best 임을 고려할 수 있다 ! (물론 지금은 다 autoarima 를 돌리겠지만 말이다..)

 

 

 

또한 값들이 잘 떨어졌으니 1차 차분을 한 데이터가 정상성을 띄게 되어 이제 ARIMA 모델에 넣어볼 수 있는 것이다! 

ARIMA 모델 패키지가 바뀌어서 저 아이를 불러와줬다. 우선은 짧은 기간만 예측을 시겨보도록 하자. 

from statsmodels.tsa.arima_model import ARIMA
import statsmodels.api as sm
from pandas import datetime

# fit model
model = sm.tsa.arima.ARIMA(ts, order=(1,0,1))
model_fit = model.fit()

# predict
#특정 기간을 지정해주고 싶은 경우
start_index = datetime(1944, 6, 25)
end_index = datetime(1945, 5, 31)
forecast = model_fit.predict(start=start_index, end=end_index) # start=start_index, end=end_index

# visualization
plt.figure(figsize=(22,10))
plt.plot(weather_bin.Date,weather_bin.MeanTemp,label = "original", c= 'red')
plt.plot(forecast,label = "predicted")
plt.title("Time Series Forecast")
plt.xlabel("Date")
plt.ylabel("Mean Temperature")
plt.legend()
plt.show()

 

 

 

 

이번에는 전체값을 예측해 보겠다 

 

# Predict all path

# Fit model
model2 = sm.tsa.arima.ARIMA(ts, order=(1,0,1))
model_fit2 = model2.fit()
forecast2 = model_fit2.predict()

# Calculate error
from sklearn.metrics import mean_squared_error
error = mean_squared_error(ts, forecast2)
print("Mean Squared Error" ,error)

# Visualization
plt.figure(figsize=(22,10))
plt.plot(weather_bin.Date,weather_bin.MeanTemp,label = "original", c= 'red')
plt.plot(forecast2,label = "predicted")
plt.title("Time Series Forecast")
plt.xlabel("Date")
plt.ylabel("Mean Temperature")
plt.legend()
plt.show()

 

전반적으로 잘 smoothing 된 값을 예측하는 양상을 보인다. 사실 mse 찍어봤을 때 저정도면 높은 정확성을 보이는 모형이기 때문에 여기서 멈출 수도 있지만, autoarima (파라미터 튜닝) 을 아래 코드를 통해 할 수 있다 

 

from itertools import product

p = range(0,3)
d = range(1,2)
q = range(0,3)
pdq = list(product(p,d,q))  #조합 만들기

#최적의 파라미터 찾기
# AIC score 는 낮을수록 좋음. (2,1,2)
aic =[]
for i in pdq:
  model = ARIMA(y_train.TotalExport, order=(i))
  model_fit = model.fit()
  print(f'ARIMA:{i} >> AIC : {round(model_fit.aic,2)}')
  aic.append(round(model_fit.aic,2))

 

 

 

4.5. 평가지표 

 

 

시계열 데이터는 MAE, MSE, RMSE  등의 지표를 평가지표로 삼는다. 왜냐면! 이건 구하는 공식 자체가 예측값과 실제값의 오차를 통해 구하기 때문이다. 반면에, R2 score 는 linear 회귀 모델에 적합한 평가지표다. 시계열 데이터는 기본적으로 자기회귀를 띄고, 이는 r2 score 의 변수 독립성 가정을 파괴하기 때문이다. 때문에 r2 는 참고 자료로만 사용하고, 모델을 튜닝하는 과정에서는 MAE, MSE 를 쓰는게 좋다. 

 

from sklearn import metrics

def mae(y_true, y_pred):
    return metrics.mean_absolute_error(y_true,y_pred) #MAE
def mse(y_true, y_pred):
    return metrics.mean_squared_error(y_true,y_pred) # MSE
def rmse(y_true, y_pred):    
    return np.sqrt(metrics.mean_squared_error(y_true,y_pred))  # RMSE
def r2(y_true, y_pred):    
    return metrics.r2_score(y_true,y_pred) # R2
def mape(y_true, y_pred):
    return np.mean(np.abs((y_pred - y_true) / y_true)) * 100 # MAPE

def get_score(model, y_true, y_pred):
    model = model
    mae_val = mae(y_true, y_pred)
    mse_val = mse(y_true, y_pred)
    rmse_val = rmse(y_true, y_pred)
    r2_val = r2(y_true, y_pred)
    mape_val = mape(y_true, y_pred)
    
    score_dict = {"model": model,
                  "mae" :  mae_val,
                  "mse" :  mse_val,
                  "rmse" : rmse_val,
                  "r2":    r2_val, 
                  "mape" : mape_val
                 }
    return score_dict

 

 

 

 

 

 

5. 요약 

 

 

- ARIMA 는 불안정 시계열의 대표적인 모델 . (p, d, q) 세가지 차수로 구성되어있으며, AR 모델, MA 모델, 차분 세 가지가 모델을 결정하는 요소 

 

- 차수에 대한 결정은 ACF, PACF 를 이용해 상관성을 검증한 뒤 할 것

 

- 평가지표로는 시계열 데이터에 적합한 MAE, MSE 를 이용

728x90
반응형
저작자표시 (새창열림)