2021-03-30(시계열데이터 심화10_ARIMA&SARIMA)

*본 내용은 시계열데이터 강의 내용 중 일부분을 요약한 내용입니다

 

 

 

ARIMA의 한계를 실습을 통해 확인해 보기

# ARIMA 모형의 한계
# 데이터 로딩 및 시간변수 반영
raw_set = datasets.get_rdataset("accdeaths", package='MASS')
raw = raw_set.data
raw.time = pd.data_range("1973-01-01", periods= len(raw), freq='M')
raw['month'] = raw.time.dt.month

#데이터 확인 
display(raw.tail())
plt.plot(raw.time, raw.value)
plt.show()

# 정상성 확인
display(stationarity_adf_test(raw.value, [])) #비정상(추세가 아직 있다고 판단)

#계절성이 있는 것처럼 보이는데 계절성보단 사이클이라서 이런 결과를 받을수있다
display(stationarity_kpss_test(raw.value, [])) #정상(계절성이 없다고 판단)

#계절성 있음
sm.graphics.tsa.plot_acf(raw.value, lags= 50, use_vlines=True)
plt.tight_layout()
plt.show()

'''
kpss는 계절성이 없다는 결과가 나왔지만 ACF그래프에선 계절성이 뚜렷하게 보인다.
ARIMA는 추세차분는 가능하지만 계절성 차분은 포함하지 않는다. 따라서 계절성까지 커버할 모델이 바로
SARIMA이다. 
'''

 

 

 

SARIMA(Seasonal ARIMA)

SARIMA는 두개의 ARMA를 가지고 있습니다. 추세영역의 ARMA 과 계절성 영역의 ARMA입니다. 

Non-seasonal of the model의 p,d,q는 추세 차분(d) 및 p(AR), q(MA) 입니다. 이렇게 ARMA와 추세차분d로 구성되있고 계절성도 같은 구성이지만 하나 더 추가 된 점은 빈도를 뜻하는 m입니다. 계절성이 만일 12개월마다 나타나면 m=12가 됩니다. 

두 영역으로 구성된 SARIMA는 데이터를 다룰때 자유자재로 추세와 계절성에 대해서 유연하게 대처할수 있습니다.

 

 

위 이미지에서는 오직 계절성의 MA(1)과  m=12를 다룹니다. MA는 화이트노이즈 영역이므로 작년 동월 백색잡음이 더해집니다. 

import pandas as pd
import numpy as np
from statsmodels import datasets
import matplotlib.pyplot as plt
import statsmodels.api as sm

# SARIMA(0,0,0)(0,0,1,12) : MA(1)->ACF보기
plt.figure(figsize=(10,4))
np.random.seed(123)
ar_params = np.array([])
ma_params = np.array([0,0,0,0,0,0,0,0,0,0,0,0.95]) #m=12이므로 12번째에서 값을 주고 나머지는 0
ar,ma = np.r_[1, -ar.params], np.r_[1, ma_params]
y = sm.tsa.ArmaProcess(ar,ma).generate_sample(500, burnin=50)

#12배수마다 peak가 친다
plt.plot(y, 'o-')
plt.tight_layout()
plt.show()

#이론적인 그래프
plt.figure(figsize=(10,9))
plt.subplot(411)
plt.stem(sm.tsa.ArmaProcess(ar,ma).acf(lags=50))
plt.xlim(-1,51)
plt.ylim(-1.1,1.1)
plt.title("Theoretical ACF of an SARIMA process")

plt.subplot(412)
plt.stem(sm.tsa.ArmaProcess(ar, ma).pacf(lags=50))
plt.xlim(-1, 51)
plt.ylim(-1.1, 1.1)
plt.title("Theoretical PACF of an SARIMA process")

#샘플링한 그래프
sm.graphics.tsa.plot_acf(y, lags=50, ax=plt.subplot(413))
plt.xlim(-1, 51)
plt.ylim(-1.1, 1.1)
plt.title("Sampling ACF of an SARIMA process")
'''
m=12이므로 12번째에서 peak가 발생
'''

sm.graphics.tsa.plot_pacf(y, lags=50, ax=plt.subplot(414))
plt.xlim(-1, 51)
plt.ylim(-1.1, 1.1)
plt.title("Sampling PACF of an SARIMA process")
plt.tight_layout()
plt.show()
'''
이론과 다르게 PACF에서 튀는 값들이 보인다
이를 통해 알고 넘어가야 하는 부분은 완벽하게 이론처럼 나오진 않지만 이론과 비슷한
패턴은 실제 데이터에서 나올수 있겠구나 정도로 이해하기
'''

 

 

 

이번엔 MA(1) 대신 AR(1)을 입력합니다. AR는 과거 시점의 관측값으로 현재 값을 표현하기때문에

작년 동월 자료값이 더해집니다.

Yt를 분배해주면 Yt-Y(t-12) = 𝜖𝑡 가 되고 이를 이항하면 위 그림처럼 정리가 됩니다. 

# SARIMA(0,0,0)(1,0,0,12) AR(1)->PACF
plt.figure(figsize=(10, 4))
np.random.seed(123)
ar_params = np.array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.95])
ma_params = np.array([])
ar, ma = np.r_[1, -ar_params], np.r_[1, ma_params]
y = sm.tsa.ArmaProcess(ar, ma).generate_sample(500, burnin=50)

#12의 배수가 될때마다 0.99*Y(t-12)+현재백색잡음+(t-12)백색잡음이 더해진다
plt.plot(y,'0-')
plt.tight_layout()
plt.show()

plt.figure(figsize=(10,9))
plt.subplot(411)
plt.stem(sm.tsa.ArmaProcess(ar,ma).acf(lags=50))
plt.xlim(-1,51)
plt.ylim(-1.1,1.1)
plt.title("Theoretical ACF of an SARIMA process")

plt.subplot(412)
plt.stem(sm.tsa.ArmaProcess(ar, ma).pacf(lags=50))
plt.xlim(-1, 51)
plt.ylim(-1.1, 1.1)
plt.title("Theoretical PACF of an SARIMA process")

sm.graphics.tsa.plot_acf(y, lags=50, ax=plt.subplot(413))
plt.xlim(-1, 51)
plt.ylim(-1.1, 1.1)
plt.title("Sampling ACF of an SARIMA process")

sm.graphics.tsa.plot_pacf(y, lags=50, ax=plt.subplot(414))
plt.xlim(-1, 51)
plt.ylim(-1.1, 1.1)
plt.title("Sampling PACF of an SARIMA process")
plt.tight_layout()
plt.show()

'''
Theoretical PACF of an SARIMA process의 그래프를 보면 12번째의 preak가 나오고 나머지는 0이 나왔습니다.
이는 이론적인 내용입니다. 실제 샘플링한 PACF를 보게되면 12번째 이외에도 튀는 값들이 여러 존재합니다.
현실은 이와 같은 패턴이 나올수 있음을 숙지 합니다. 
'''

 

 

계절성 차분이 1회 & Ma(1)가 들어가 있습니다. 차분이 있다는 얘기는 원본 데이터가 비정상성 이라는 뜻을 내포합니다. 이 내용을 식으로 표현하면 위와 같습니다. 왼쪽항은 계절성 차분에 대한 내용이고 오른쪽 항은 MA(1)입니다. 차분을 하게되면 AR에서 나왔던 관측값이 만들어집니다. Y(t-12) 

# SARIMA(0,0,0)(0,1,1,12)
plt.figure(figsize=(10, 4))
np.random.seed(123)
ar_params = np.array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.99])
ma_params = np.array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.95])
ar, ma= np.r_[1, -ar.params], np.r_[1, ma_params]
y = sm.tsa.ArmaProcess(ar, ma).generate_sample(500, burnin=50)
plt.plot(y, '0-')
plt.tight_layout()
plt.show()

plt.figure(figsize=(10,9))
plt.subplot(411)
plt.stem(sm.tsa.ArmaProcess(ar,ma).acf(lags=50))
plt.xlim(-1, 51)
plt.ylim(-1.1, 1.1)
plt.title("Theoretical ACF of an SARIMA process")
'''
계절성 차분이 있기 때문에 값이 누적되었다
12 시차 마다 누적이 된다.. 왜 누적이 될까?.. (좀더 공부가 필요하다)
'''
plt.subplot(412)
plt.stem(sm.tsa.ArmaProcess(ar, ma).pacf(lags=50))
plt.xlim(-1, 51)
plt.ylim(-1.1, 1.1)
plt.title("Theoretical PACF of an SARIMA process")

sm.graphics.tsa.plot_acf(y, lags=50, ax=plt.subplot(413))
plt.xlim(-1, 51)
plt.ylim(-1.1, 1.1)
plt.title("Sampling ACF of an SARIMA process")

sm.graphics.tsa.plot_pacf(y, lags=50, ax=plt.subplot(414))
plt.xlim(-1, 51)
plt.ylim(-1.1, 1.1)
plt.title("Sampling PACF of an SARIMA process")
plt.tight_layout()
plt.show()
'''
앞쪽은 stable한 모습을 보이는데 뒤로 갈수록 튀는 값들이 많아졌다.
이는 추세 적합만되고 계절성 적합되지 않을때 발생->사리마 적용해야한다
'''

 

 

식을보면 Yt-24가 있기 때문에 PACF를 보게되면 Yt-12, Yt-24 자리에 corr값이 존재 할 것임을 알수 있습니다. 

실제 샘플링으로 plot한 그래프는 이론적인 모습처럼 깔끔하게 나와있지 않고 다영한 peak들이 존재합니다. 

SARIMA가 다양한 ACF,PACF의 특성을 반영하고 있기 때문입니다.

# SARIMA(0,0,0)(1,1,0,12)
plt.figure(figsize=(10,4))
np.random.seed(123)
ar_params= np.array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.99, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.95])
ma_params = np.array([])
ar, ma = np.r_[1, -ar_params], np.r_[1, ma_params]
y = sm.tsa.ArmaProcess(ar,ma).generate_sample(500, burnin=50)
plt.plot(y,'-0')
plt.tight_layout()
plt.show()

plt.figure(figsize=(10,9))
plt.subplot(411)
plt.stem(sm.tsa.ArmaProess(ar, ma).acf(lags=50)
plt.xlim(-1, 51)
plt.ylim(-1.1, 1.1)
plt.title("Theoretical ACF of an SARIMA process")

plt.subplot(412)
plt.stem(sm.tsa.ArmaProcess(ar, ma).pacf(lags=50))
plt.xlim(-1, 51)
plt.ylim(-1.1, 1.1)
plt.title("Theoretical PACF of an SARIMA process")

sm.tsa.plot_acf(y, lags=50, ax=plt.subplot(413))
plt.xlim(-1, 51)
plt.ylim(-1.1, 1.1)
plt.title("Sampling ACF of an SARIMA process")

sm.graphics.tsa.plot_pacf(y, lags=50, ax=plt.subplot(414))
plt.xlim(-1, 51)
plt.ylim(-1.1, 1.1)
plt.title("Sampling PACF of an SARIMA process")
plt.tight_layout()
plt.show()

'''
12배수 이외에도 자잘한 peak가 존재한다.
SARIMA부터는 다양한 ACF,PACF의 특성을 반영하고 있다. 
'''

 

위 경우에는 추세 MA(1)의 영향으로 ACF plot은 12번째 앞뒤인 11번째,13번째가 같은 크기로 peak를 칩니다. 

(11번째의 상관성 증가에 대한 계산 방식이나 원인은 자세히 모르겠습니다ㅠ 자세한 설명을 못드린 점에서 양해 부탁드리겠습니다)

실습을 통해 내용을 확인해봅니다.

# SARIMA(0,0,1)(0,0,1,12) 
plt.figure(figsize=(10,4))
np.random.seed(123)
phi, Phi = 0, 0 
theta, Theta = 0.5, 0.8
ar_params = np.array([])
ma_params = np.array([0,0,0,0,0,0,0,0,0,0,0,Theta, theta*Theta])
ar, ma = np.r_[1, -ar_params], np.r_[1, ma_params]
y = sm.tsa.ArmaProcess(ar, ma).generate_sample(500, burnin=50)
plt.plot(y, '-0')
plt.tight_layout()
plt.show()

plt.figure(figsize=(10,9))
plt.subplot(411)
plt.stem(sm.tsa.ArmaProcess(ar, ma).acf(lags=50))
plt.xlim(-1,51)
plt.ylim(-1.1,1.1)
plt.title("Theoretical ACF of an SARIMA process")
'''
ACF에서 11,12,13번째 값이 튀는 현상을 볼수 있다
'''

plt.subplot(412)
plt.stem(sm.tsa.ArmaProcess(ar, ma).pacf(lags=50))
plt.xlim(-1, 51)
plt.ylim(-1.1, 1.1)
plt.title("Theoretical PACF of an SARIMA process")
'''
12배수마다 앞뒤 값이 함께 튀고 있다
'''
sm.graphics.tsa.plot_acf(y, lags=50, ax=plt.subplot(413))
plt.xlim(-1, 51)
plt.ylim(-1.1, 1.1)
plt.title("Sampling ACF of an SARIMA process")

sm.graphics.tsa.plot_pacf(y, lags=50, ax=plt.subplot(414))
plt.xlim(-1, 51)
plt.ylim(-1.1, 1.1)
plt.title("Sampling PACF of an SARIMA process")
plt.tight_layout()
plt.show()

 

 

 

계절주기에는 AR(1), 비계절에는 MA(1)인 경우입니다.

계절주기는 ACF를 보는데, 위 같이 앞뒤가 같이 상관성을 보이는 경우가 ACF에서 계절성 패턴처럼 주기적으로(12배수) 나타나게 됩닌다. 추세의 MA(1)는 PACF를 봐야 하므로 11~13번 튀는 현상을 보인 후 나머지 계수는 0으로 나타나게 됩니다. 

실습을 통해 확인해 봅니다

# SARIMA(0,0,1)(1,0,0,12)
plt.figure(figsize(12,4))
np.random.seed(123)
phi,Phi = 0,0.75
theta, theta = 0.5,0 
ar_params = np.array([0,0,0,0,0,0,0,0,0,0,0,Phi])
ma_param = np.array([theta])
ar, ma = np.r_[1, -ar_params], np.r_[1, ma_params]
y = sm.tsa.ArmaProcess(ar,ma).generate_sample(500, burnin=50)
plt.plot(y, '-0')
plt.tight_layout()
plt.show()

plt.figure(figsize=(10,9))
plt.subplot(411)
plt.stem(sm.tsa.ArmaProcess(ar,ma).acf(lags = 50))
plt.xlim(-1,51)
plt.ylim(-1.1,1.1)
plt.title("Theoretical ACF of an SARIMA process")

plt.figure(figsize=(10,9))
plt.subplot(412)
plt.stem(sm.tsa.ArmaProcess(ar,ma).pacf(lags=50))
plt.xlim(-1, 51)
plt.ylim(-1.1, 1.1)
plt.title("Theoretical PACF of an SARIMA process")

sm.graphics.tsa.plot_acf(y, lags=50, ax = plt.subplot(413))
plt.xlim(-1,51)
plt.ylim(-1.1,1.1)
plt.title("Sampling ACF of an SARIMA process")

sm.graphics.tsa.plot_pacf(y, lags = 50, ax=plt.subplot(414))
plt.xlim(-1, 51)
plt.ylim(-1.1, 1.1)
plt.title("Sampling PACF of an SARIMA process")
plt.tight_layout()
plt.show()

 

실습: 호흡기질환 사망자수 SARIMA 모델링

SARIMA 모델을 적용하는 것을 순서대로 실습을 하기에 앞서서 적용 순서를 정리해 봤습니다

 

step1. 데이터가 비정상성이라면 transformation, differencing(차분)을 통해 정상성으로 바꿔준다
step2. ACF,PACF을 통해 어떤 모델을 사용할건지를 판단한다.(차분이 된 상태에서 어떤 패턴이 있는지를 추론 하는것)
step3. ACF,PACF을 통해 파라미터를 판단하기 -> p or q를 추론하면서 AIC,BIC 성능이 가장 좋은 것을 찾는다. 하지만 현실에서는 for문을 통해 p,q의 경우의 수를 모두 돌려본다(보통 1~20까지 돌려본다 )
step4. residual(실제값-예측값)을 가지고 ACF,PACF 그려본다. 이때 범위가 생기는 이 범위는 residual 분산의 3시그마이다. 대략 40개중 2~3개 까지는 OK 입니다.
step5. 진단 결과 범위안에 잘 들어왔으면 최종모델로 선택, 범위를 벗어난다면 2번 단계로 돌아가서 모델을 변경

 

import warnings
warnings.filterwarnings('always')
warnings.filterwarnings('ignore')

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import statsmodels.api as sm

#데이터 로딩 및 확인
data = sm.datasets.get_rdataset("deaths","MASS")
raw= data.data
raw.value= np.log(raw.value)
raw.plot(x='time', y='value')
plt.show()

#ACF/PACF 확인
plt.figure(figsize=(10,8))
sm.graphics.tsa.plot_acf(raw.value, lags=50, ax=plt.subplot(211))
plt.xlim(-1,51)
plt.ylim(-1.1,1.1)
plt.title("ACF")
plt.show()

plt.figure(figsize=(10,8))
sm.graphics.tsa.plot_pacf(raw.value, lags=50, ax=plt.subplot(211))
plt.xlim(-1,51)
plt.ylim(-1.1,1.1)
plt.title("PACF")
plt.show()
'''
위 시각화를 통해 모델및 파라미터를 추론할수 있다
ACF에서 계절성이 뚜렷하게 보인다. PACF를 보면 뒤로 갈수록 불안정해지는 특성을 봐서는 
계절성 패턴이 있음을 짐작케 한다. -> 계절성 차분이 필요해 보임

ACF에서 lag1이 가장 큰 값을 가지므로 MA(1) + PACF도 lag1에서 가장 큰 값을 가지므로 AR(1)
모델을 시험삼아 돌려본 후 진단을 해보도록한다
'''

#ARMA(1.1)모델링
fit = sm.tsa.SARIMAX(raw.value, trend='c', order=(1,0,1), seasonal_order=(0,0,0,0)).fit()
display(fit.summary())
'''
시그마: 잔차&상수
Jarque-Bera:정규분포가 아니다
왜도첨도: 정규분포에 가깝다
Heteroskedasticity(H): 등분산이다. 분산이 일정하다
Ljung-Box:자기상관성이 있다.
요약: 정규분포와 비슷하나 정규분포는 아니고 자기상관성은 있고 분산은 일정하다
'''

 

 

 

계절성 order를 추가하여 실습을 진행한 결과 아래와 같다

# SARIMA 모델링
fit = sm.tsa.SARIMAX(raw.value, trend='c', order =(1,0,1), seasonal_order=(1,1,1,12)).fit()
display(fit.summary())

# 잔차진단
fit.plot_diagnostics(figsize=(10,8))
plt.tight_layout()
plt.show()
"""
정규분포가 아니라고 또 다시 나왔는데 이는 아마도 
아웃라이어 때문일 가능성이있다. 
그리고 변수가 유의하지 않는 것들이 많은 것도 이 원인일 가능성이 있다.
"""