2021-03-29(시계열데이터 심화9_ARMR&ARIMA)

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

 

 

MA모형(Moving Average)

위 처럼 과거 시점의 오차를 사용하여 Yt를 구한다. 오차는 실측치 와 예측치를 뺀 값이다. 어떤 변화가 보일때는 이 오차 값이 커진다. 반대로 변동이 적고 충분히 예측이 가능하다면 오차는 0에 가깝게 됩니다. MA는 과거 데이터의 변동을 가지고 미래가 어떤 식으로 변할건가를 알아보는 방법입니다. 

 

MA의 ACF는 q시점까지 상관성이 있고, PACF는 완만한 감소 패턴을 보인다. 이를 아래 실습으로 확인해보자

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

MA의 q가 5일때
plt.figure(figsize=(10,4))
np.random.seed(123)
ar_params= np.array([])

#t-1의 과거 시점 오차의 계수 -1, t-2의 과거 시점 오차의 계수 -1.6, ...,t-5의 과거 시점 오차의 계수 0.7,
ma_params= np.array([-1,1.6,09.,-1.5,0.7])
ar, ma = np.r_[1, -ar_params], np.r_[1,ma_params]
y = sm.tsa.ArmaProcess(ar,ma).gerate_sample(500, burnin=50)

#샘플링한 데이터 시각화
plt.plot(y, '0-')
plt.tight_layout()
plt.show()

#ACF 시각화
plt.figure(figsize=(10,6))
plt.subplot(311)
plt.stem(sm.tsa.ArmaProcess(ar,ma).acf(lags=10)) #lags는 10개 까지 실행(5개이후로는 0으로 나온다)
plt.xlim(-1,11) # X범위 설정
plt.ylim(-1.1, 1.1) # Y범위 설정 
plt.title("Theoretical autocorrelation function of an ARMA process")


#PACF 시각화
plt.figure(figsize=(10.6))
plt.subplot(312)
plt.stem(sm.tsa.ArmaProcess(ar,ma).pcaf(lags=10))
plt.xlim(-1,11)
plt.ylim(-1.1, 1.1)
plt.title("Sampling autocorrelation function of an ARMA process")
plt.tight_layout()
plt.show()

#plot_acf로 ACF 시각화
sm.graphics.tsa.plot_acf(y, lags=10, ax=plt.subplot(313))
plt.xlim(1,11)
plt.ylim(-1.1, 1.1)
plt.title("Sampling autocorrelation function of an ARMA process")
plt.tight_layout()
plt.show()

 

AR(Auto-Regressive)

변수의 과거 값의 선형 조합을 이용하여 관심 있는 변수를 예측합니다.

1) a <0 일때(음수): 부호를 바꿔가면서 진동하는 패턴을 보이며 감소한다. 제곱은 +로, 세제곱은 -로 위아래 반복

2) a >0 일때(양수): 차수(q)가 증가하면서 자기상관계수는 지수적으로 감소(과거 관측값으로 갈수록 상관성이 감소)

3) a=1 일때(단위근): a=1일때를 단위근이라 말하며, 이는 비정상성으로써 발산하는 형태를 보이는 random walk

q가 늘어 날수록 a값들의 제한조건들이 추가가 됩니다. 여기서 et는 현재 시점의 잔차이자 백색잡음입니다.

여기서 a가 𝜙1라고 가정할때, 이 값이 0.9라고 한다면 첫번째 corr값은 0.9 그대로 나온다. 두번째 corr값은 0.9^2=0.81로 출력된다. 이렇게 세제곱, 네제곱, ..., 무한대로 간다

위 무한대로 뻗어가면 q가 몇개여야 하는지 추론을 하지 못하게 된다. 이때 PACF를 통해 확인하게 되면 영향을 주는 계수만 값이 나오고, 나머지 과거 시점은 다 0으로 떨어진다. 아래 AR(1)의 실습으로 확인해본다.

###AR(1)
plt.figure(figsize=(10,4))
np.random.seed(123)

ar_params = np.array([0.9])
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,6))
plt.subplot(311)
plt.stem(sm.tsa.ArmaProcess(ar, ma).acf(lags=10))
plt.xlim(-1,11)
plt.ylim(-1.1, 1.1)
plt.title('Theoretical autocorrelation function of an ARMA process')

#PACF에서는 계수추론이 가능하다
plt.subplot(312)
plt.stem(sm.tsa.ArmaPrcess(ar, ma).pacf(lags=10))
plt.xlim(1,11)
plt.ylim(-1.1,1.1)
plt.title('Theoretical partial autocorrelation function of an ARMA process')

sm.graphics.tsa.plot_pacf(y, lags=10, ax=plt.subplot(313))
plt.xlim(-1,11)
plt.ylim(-1.1,1.1)
plt.title('Sampling partial autocorrelation function of an ARMA process')
plt.tight_layout()
plt.show()

정리

1) MA 무한대 = AR(1)

2) MA의 계수추론은 ACF, AR의 계수추론은 PACF

ARMA(Auto-Regressive Moving Average)

AR+MA를 결합된 모델로써 두개의 패턴이 같이 나올 가능성이 있다. 

• ϕϕ: 자기회귀 모형의 모수
• θθ: 이동평균 모형의 모수
• c: 상수
• ϵϵ: 오차항 (백색잡음; white noise)

### ARMA(1,1)
#데이터의 트렌드를 한번 간단하게 뽑아낼때 ARMA(1,1)와 같은 걸로 많이 적합한다.
# Setting

np.random.seed(123)
ar_params = np.array([0.75])
ma_params = np.array([0.65])
ahead=100
ar, ma = np.r_[1,-ar_params], np.r_[1, maparams]
ar_order, ma_order = len(ar)-1, len(ma)-1 #1,1

#generator
y = statsmodels.tsa.arima_process.arma_generate_sample(ar,ma, nsample=1000, burnin=500)
fit = statsmodels.tsa.arima_model.ARMA(y, (ar_order, ma_order)).fit(trend='c', disp=0)

pred_ts_point = fit.forecast(steps=ahead)[0] #점추정값
pred_ts_interval = fit.forecast(steps=ahead)[2] #구간추정
ax = pd.DataFrame(y).plot(figsize=(12,5))
forecast_index =[i for i in range(pd.DataFrame(y).index.max()+1, pd.DataFrame(y).index.max()+ahead+1)]
ax.fill_between(pd.DataFrame(pred_ts_interval, index = forecast_index).index,
				pd.DataFrame(pred_ts_interval, index=forecast_index).iloc[:,0],
                pd.DataFrame(pred_ts_interval, index=forecast_index).iloc[:,1], color='k',alpha=0.15)
plt.legend(['observed', 'forecast'])
display(fit.summary2())
plt.figure(figsize=(12,3))
statsmodels.graphics.tsaplots.plot_acf(y, lags=50, zero=True, use_vlines=True, alpha=0.05, ax=plt.subplot(121))
statsmodels.graphics.tsaplots.plot_pacf(y,lags=50, zero=True, use_vlines=True, alpha=0.05, ax=plt.subplot(122))
plt.show()

'''
MA의 PACF 특성(무한대로 가는 꼴) & AR의 ACF 특성(무한대 꼴)이 
따라서 두개의 arma를 사용하기 시작하면 어떠한 계수가 적합한지 쉽게 의사결정하기 어렵다
-> 통계량, 가설검정으로 종합적으로 판정한다
'''

 

파라미터를 위 실습내용보다 많이 늘리게 되면 오히려 파라미터 끼리 관계가 짙어져 과적합을 일으켜서 자기 상관성이 커지는 현상이 발생 할수 있다. 그래서 모델링할때는 파라미터 사용에 유의해야한다.

 

 

ARMAX(ARMA with eXogenous)

기존의 arma는 자기자신의 과거 관측치 또는 오차만을 가지고 했기때문에 정확성 측면에서는 좋은 성능을
발휘하지 못했다. 이를 개선시키고자 적합할때 X(설명변수)까지 고려한 모델이 armax이다.

 

강의 자료

위 자료를 보면 ARMAX  Yt=f(Y(t-p), 𝜖(t-q),x)+𝜖t이다.

즉, ARMA와 X(설명변수)까지 고려한 추정값에 현재시점의 오차를 더해준 값이다. 이를 통해 정확성을 개선시킬수 있다.

 

import requests
from io import BytesIO
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import statsmodels.api as sm

# 데이터 로딩 및 확인
# m2:통화량, 화페량
# 실제 국민이 소비한 화폐량
source_url = requests.get('http://www.stata-press.com/data/r12/friedman2.dta').content
raw = pd.read_stata(BytesIO(source_url))
raw.index = raw.time
raw_using = raw.loc['1960':'1981',["consump", "m2"]]
raw_using.plot()
plt.show()

# 모델링
## 회귀분석
fit = sm.OLS(raw_using.consump, sm.add_constant(raw_using.m2)).fit()
display(fit.summary())
'''
Prob(JB):0.716으로 정규분포가 아니라는 결과가 나옴
Durbin-Watson:0.094 positive correlation 경향이 있다.
DW: 오차항에 자기상관이 있느니 없느니를 판단하기 위한 검정법
value=1 : 양의 자기상관이 있음(통계량은 0)
value=0 : 자기상관이 없다(통계량은 2)
value=-1 : 음의 자기상관이 있음(통계량은4)
'''
## 잔차 확인
fit.resid.plot()
plt.show()

## 잔차 ACF/PACF
plt.figure(figsize=(10,8))
sm.graphics.tsa.plot_acf(fit.resid, lags =50, ax=plt.subplot(211))
plt.xlim(-1,51)
plt.ylim(-1.1,1.1)
plt.title('Residual ACF')
plt.tight_layout()
plt.show()

sm.graphics.tsa.plot_pacf(fit.resid, lags=50, ax=plt.subplot(212))
plt.xlim(-1, 51)
plt.ylim(-1.1, 1.1)
plt.title("Residual PACF")
plt.tight_layout()
plt.show()
'''
ACF :MA기준에서는 네개의 파라미터 
PACF : AR기준에서는 한 개의 파라미터
즉, arma(1,4)를 써야겠구나 생각하기
'''
==================================================
#ARIMAX 모델링
#데이티 X를 raw_using.m2로 직접 지정하기 -> exog
fit = sm.tsa.ARMA(raw_using.consump, (1,1), exog= raw_using.m2).fit()
display(fit.summary())

## 잔차 확인
'''
추세는 없지만 최근일수록 분산이 커진다
'''
fit.resid.plot()
plt.show()

## 잔차 ACF/PACF
plt.figure(figsize=(10, 8))
sm.graphics.tsa.plot_acf(fit.resid, lags=50, ax=plt.subplot(211))
plt.xlim(-1, 51)
plt.ylim(-1.1, 1.1)
plt.title("Residual ACF")

sm.graphics.tsa.plot_pacf(fit.resid, lags=50, ax=plt.subplot(212))
plt.xlim(-1, 51)
plt.ylim(-1.1, 1.1)
plt.title("Residual PACF")
plt.tight_layout()
plt.show()
'''
앞서 arma(1,4)가 적합하다고 했지만 (1,1)만으로도 상당히 안정적인 형태를 취하고 있다.

핵심은..
회귀분석과 비교해서 AIC,BIC가 200개 정도 줄어들었고
화이트 노이즈 경향도 훨씬 강해졌다

왜이러는가?
회귀분석 y= C+ax+e
arma yt=C+a*y(t-1)+b*e(t-1)+ax+e 
회귀분석의 C자리에 추가적으로 arma 모형을 넣어서 좀더 많은 패턴을 잡을 수 있도록 해준다
회이트 노이즈가 될 확률을 높인다.
'''

회귀분석_summay

회귀분석_잔차plot

ARIMA_summary

ARIMA_잔차 plot

ARIMA_ACF,PACF

 

위 결과에서 상당히 개선된 모습을 봤다. 그런데 ARIMA의 plot들을 보면 대체로 뒤로 갈수록 값이 커짐을 볼수있다. 특히 잔차 plot과 PACF에서 뚜렷히 볼수가 있는데 이 현상의 원인 여러 원인이 있겠지만 주 원인은 제거 되지 않은 계절성 때문이다. ARMA는 추세패턴을 뽑아내는 모형이다. 그래서 적합 시킬때 추세패턴을 없앨 수 있었지만 계절성 패턴은 남아 있어서 위 처럼 영향을 줄 수 있다. 이를 해결 하기 위한 모형이 SARIMAX이다.

S가 seasonal를 뜻한다. 실습을 통해 자세히 알아본다

 

'''
lag값이 뒤로 갈수록 corr이 커지는 현상은 왜 나타는가?
ARMA는 추세를 제거하는게 주 목적이다. 자기자신의 값과 화이트 노이즈를 쓰기때문에
추세가 제거가 된거다. 하지만 계절성은 제거하지 못했다.

따라서, 궁극적으로 이 둘다 줄일수 있는 것으로써 sarimax를 사용한다.
sarimax의 s는 seasonal
'''

# SARIMAX 모델링
#추세 arma => (1, 0.0)
#계절성 s => (1,0,1,4)
fit = sm.tsa.SARIMAX(raw_using.consump, exog=raw_using.m2, order=(1,0,0), seasonal_order=(1,0,1,4)).fit()
display(fit.summary())
'''
ar.S.L4,  ma.S.L4
S가 붙은건 계절성을 뽑아서 적합을 시킨것
'''

'''summary 해석
등분산성(Heteroskedasticity )은 아니지만 자기상관(Ljung-Box)은 없고 정규분포(자크베라)에 가까우며 왜도,첨도로 정규분포의 모양을 확인할수 
있으며 마지막으로 ACF,PACF로 잔차가 화이트노이즈에 가깝다는 걸 확인
'''

#잔차확인
fit.resid.plot()
plt.show()
#추세/계절성 패턴이 없으며 분산이 조금은 있어보인다 -> 통계량을 통해 등분산성 인지 확인

#잔차 ACF/PACF
plt.figure(figsize=(10,8))
sm.graphics.tsa.plot_acf(fit.resid, lags=50, ax=plt.subplot(211))
plt.xlim(-1,51)
plt.ylim(-1.1,1.1)
plt.title("Residual PACF")
plt.tight_layout()
plt.show()
'''
잔차그림에서는 위 시각화에 비해 안정적으로 바꿨다.
ACF,PACF가 완벽하게 화이트 노이즈임을 보인다
'''

 

 

ARIMA(Auto-Regressive Integrated Moving Average)

지금까지는 비정상 데이터(추세+계절성)를 차분을 통해 제거 후 정상성이 된 데이터를 ARMA로 적합시키는 과정을 배웠다. ARIMA는 추세+계절성 제거 &  ARMA를 동시에 하는 모형이다. 주의할건 seasonal은 없기 때문에 오직 추세차분만제거할수 있다.

위 정의에서 차분이라는 말이 나오는데 이는 추세 차분만을 말한다. 식에서 d가 바깥에서 ^d 형태를 가진다. 이는 추세 차분을 뜻한다. 따라서 ARIMA는 추세차분을 1회 이상 한 것을 말하며 차분이 없다면 ARMA 모형이 된다는 걸 추론할수 있다. 그리고 차분을 해야 한다는 것은 데이터가 비 정상성이라는 뜻이다. d=0이라면 애초에 데이터가 정상성을 가진 데이터라는 것이다. 

따라서 ARIMA(p,d,q) 모형에 따라서 의미가 다르다. 정리한 내용은 아래와 같다.

*오류 부분이 있는데 맨 아래 0,0,1는 음수를 뜻하므로 -1<세타<0이다.