이번 포스팅에서는 자기 회귀 모형(Autoregressive Model)과 이동 평균 모형(Moving Average Model)을 결합한 ARMA 모형에 대해서 알아보려고 한다.
1. ARMA란?
- 정의 -
시계열 데이터 $X_t$가 정상성을 만족하고 $ARMA(p, q)$ 모형을 따른다고 한다면 아래의 관계식을 만족한다.
$$X_t = c + \sum_{i=1}^p\phi_iX_{t-i} + Z_t - \sum_{i=1}^q\theta_iZ_{t-i} \tag{1.1}$$
여기서 $Z_t$는 정규분포를 따르는 백색 잡음이다. 즉, $Z_t \text{ i.i.d. } \sim N(0, \sigma^2)$ 이다.
이때 $B$를 Back Shift 연산자라고 하자. 즉, $BX_t = X_{t-1}$이고 $B^pX_t = X_{t-p}$이다. 이 연산자를 이용하면 식 (1.1)을 아래와 같이 표현할 수 있다.
$$(1-\phi_1B-\cdots \phi_pB^p)X_t = c + (1-\theta_1B - \cdots - \theta_qB^q)Z_t\tag{1-1}$$
- 필요성 -
자기 회귀 모형 또는 이동 평균 모형 하나만 가지고 시계열 데이터를 설명하기 위하여 많은 파라미터를 사용해야 할 수 있지만 ARMA를 이용하면 모수 절약 효과를 볼 수 있다. 왜 그런지 살펴보자.
먼저 $\phi (B) = 1-\sum_{i=1}^p\phi_iB^i$, $\theta (B) = 1-\sum_{i=1}^q\theta_iB^i$라 하자. 이를 이용하면 식 (1-1)은 다음과 같이 표현할 수 있다.
$$\phi (B)X_t = c + \theta(B)Z_t \tag{1-2}$$
그리고 복소수 $z$에 대하여 $\phi(z) = 1-\phi_1z - \cdots - \phi_pz^p$라 하고 $\theta(z)$도 비슷하게 정의하자.
정리 1-1
이때 $\theta$와 $\phi$의 공통근은 없다고 하자. 또한 $|z| \leq 1$인 모든 복소수 $z$에 대하여 $\phi(z)\neq 0$이라 하자. 그러면 다음을 만족하는 $\psi_1, \psi_2, \ldots$가 존재한다.
$$\frac{\theta(B)}{\phi(B)} = 1+\psi_1B+\psi_2B^2 +\cdots \equiv \psi(B)$$
이때 정리 1-1을 만족하는 $X_t$, 즉 $ARMA(p, q)$ 모형을 Causal 모형이라고 한다.
정리 1-2
이번엔 $|z| \leq 1$인 모든 복소수 $z$에 대하여 $\theta(z)\neq 0$이라 하자. 그러면 다음을 만족하는 $\pi_1, \pi_2, \ldots $가 존재한다.
$$\frac{\phi(B)}{\theta(B)} = 1 - \pi_1B - \pi_2B^2 - \cdots \equiv \pi(B)$$
이때 정리 1-2를 만족하는 $X_t$, 즉 $ARMA(p, q)$ 모형을 Invertible 하다고 한다.
또한 임의의 상수 $c$는 시간에 따라 변하지 않으므로 $Bc = c$인 것을 알 수 있다. 이제 식(1-1)을 자기 회귀(AutoRegressive : AR) 표현식과 이동 평균 (Moving Average : MA) 표현식으로 나타내 보자.
1) ARMA의 AR 표현식
식 (1-2)에서 양변에 $\theta (B)^{-1}$를 적용하자.
$$\begin{align}\frac{\phi(B)}{\theta(B)}X_t &= \frac{c}{\theta(B)\cdot 1}+Z_t \\ &= \frac{c}{1-\theta_1 -\cdots - \theta_q}+Z_t \end{align}$$
정리 1-2를 이용하면 위 식을 다음과 같이 정리할 수 있다.
$$X_t = \frac{c}{1-\theta_1 -\cdots - \theta_q} +\pi_1X_{t-1}+\pi_2X_{t-2}+\cdots+Z_t$$
Invertible 성질에 의하여 $\pi_j$ 계수들의 값은 $j$가 커질수록 0에 가까워지며 이는 과거 데이터가 현재 데이터의 미치는 영향력이 감소함을 뜻한다.
여기서 Invertibility는 $ARMA(p, q)$ 모형의 유일한 AR 표현식을 얻기 위한 필요충분조건이라고 생각하면 된다.
2) ARMA의 MA 표현식
식 (1-2)에서 양변에 $\phi (B)^{-1}$를 적용하자.
$$\begin{align}X_t &= \frac{c}{\phi(B)\cdot 1}+Z_t \\ &= \frac{c}{1-\phi_1 -\cdots - \phi_p}+\frac{\theta(B)}{\phi(B)} Z_t \end{align}$$
정리 1-1과 $(1-\phi_1-\cdots-\phi_p)\mu = c$를 이용하면 위 식을 다음과 같이 정리할 수 있다.
$$X_t = \mu +Z_t + \psi_1Z_{t-1}+\psi_2Z_{t-2}+\cdots = \mu+\psi(B) \tag{1-3}$$
여기서 $\mu = E(X_t)$이다.
위 식에서 $\psi_j$를 충격 반응 함수(Impulse Response Function)이라고 한다. Causality 성질에 의하여 $\psi_j$ 계수들의 값은 $j$가 커질수록 0에 가까워지며 이는 과거 Random Shock $Z_{t-j}$가 현재 데이터 $X_t$의 미치는 영향력이 감소함을 뜻한다.
여기서 Causality는 $ARMA(p, q)$ 모형의 유일한 MA 표현식을 얻기 위한 필요충분조건이라고 생각하면 된다.
앞에서 살펴본 ARMA의 표현식을 통해서 알 수 있는 것은 다음과 같다. 이론적으로 ARMA와 같은 설명력을 지닌 AR 모형 또는 MA 모형으로는 파라미터의 개수가 무한 개다. 즉, ARMA 모형은 같은 설명력을 지니더라도 파라미터의 개수를 획기적으로 줄여줄 수 있다.
2. ARMA 모형 추정
이제 ARMA 모형을 어떻게 추정하는지 알아보자. 추정 과정은 1) 차수 $p, q$를 정하고 2) 파라미터들을 추정해야 한다. 각 과정에 대해서 알아보자.
2.1 차수 정하기
ARMA의 차수를 정하는 두 가지 방법을 소개한다. 첫 번째는 Extended Sample Autocorrelation Function(ESACF)을 이용하는 방법과 Akaike Information Criterion(AIC)를 이용하는 방법이다.
2.1.1 ESACF
ESACF는 Tsay와 Tiao의 1984년 논문 "Consistent Estimates of Autoregressive Parameters and Extended Sample Autocorrelation Function for Stationary and Nonstationary ARMA Models"에서 소개하였다. 간단히 말해서 ESACF는 현재 데이터를 과거의 데이터 영향력을 걸러낸 부분 $W_t$의 자기 회귀 함수이다. 즉,
$$W_t = X_t - \phi_1X_{t-1} - \cdots - \phi_pX_{t-p}$$
의 Sample Autocorrelation Function이라고 생각하면 된다. ESACF에 대한 설명과 파이썬으로 구현하는 방법은 아래에 포스팅해두었으니 참고하기 바란다.
[논문 리뷰] 5. Consistent Estimates of Autoregressive Parameters and Extended Sample Autocorrelation Function for Stationary
본 포스팅에서는 수식을 포함하고 있습니다. 티스토리 피드에서는 수식이 제대로 표시되지 않을 수 있으니 PC 웹 브라우저 또는 모바일 웹 브라우저에서 보시기 바랍니다. 이번 포스팅에서는 ARM
zephyrus1111.tistory.com
ESACF를 이용한 방법은 먼저 ARMA 차수 후보를 선정한다. 이때 AR 차수 후보는 행으로, MA 차수 후보를 칼럼으로 하는 ESACF 테이블을 그린다. 테이블 하나하나의 값은 ESACF 값이다. 이때 아래와 같이 ESACF 값이 통계적으로 0에 가깝다면 'O', 아니라면 'X'로 표시한다.
이때 'O'로 표시된 영역의 패턴을 살펴본 뒤 삼각형 패턴이 나오는지 확인하고 이때 좌측 상단 'O'에 해당하는 AR, MA 차수가 ARMA의 차수로 결정된다.
2.1.2 AIC
AIC를 이용한 방법도 ESACF와 마찬가지로 여러 ARMA의 차수 후보군을 정한 뒤 가장 작은 AIC 값을 갖는 차수를 최종 ARMA의 차수로 결정한다.
AIC는 파라미터를 다소 많이 갖는 모형을 선택하는 경향성이 있어서 과적합(Overfitting)의 위험이 있다고 한다.
2.2 파라미터 추정하기
식 (1.1)에서 우리가 추정해야 할 파라미터는 $\boldsymbol{\theta} = (c, \phi_1,\ldots,\phi_p,\theta_1, \ldots, \theta_q, \sigma^2)$이다. 추정 방법은 조건부 우도 추정과 정확 우도 추정법이 있다. 정확 우도 추정은 칼만 필터(Kalman Filter)를 사용하면 된다는데 이 부분은 칼만 필터를 다룰 때 다루도록 하겠다. 여기에서는 조건부 우도 추정법에 대해서 소개한다.
조건부 우도 추정은 아래의 로그 조건부 우도 함수를 최소화하는 파라미터를 찾는다.
$$\log L(x_T, \ldots, x_{p+1}|x_p, \ldots, x_1, z_p=\cdots=z_{p-q+1}=0) \\ = -\frac{T-p}{2}\log (2\pi) - \frac{T-p}{2}\log (\sigma^2) - \sum_{t=p+1}^T\frac{z_t^2}{2\sigma^2}\tag{2.1}$$
이때 $z_t$는 아래와 같이 계산된다.
$$z_t = x_t-c-\phi_1x_{t-1}-\cdots-\phi_px_{t-p}+\theta_1z_{t-1} +\cdots+\theta_qz_{t-q}$$
조건부 우도 추정은 아래 방정식의 해 집합이 단위원 바깥에 있을 경우에 유효하다.
$$1+\theta_1z+\cdots+\theta_qz^q=0$$
(2.1)을 최소화하는 것은 Newton-Rhapson 방법과 같은 수치적인 방법을 이용한다.
3. 예측(Forecasting)
3.1 예측값
$ARMA(p, q)$ 모형을 추정했다면 이를 이용하여 예측을 해보자. $X_t$가 $ARMA(p,q)$를 따르므로 식(1.1)과 같이 아래 관계식이 성립한다.
$$X_t = c + \sum_{i=1}^p\phi_iX_{t-i} + Z_t - \sum_{i=1}^q\theta_iZ_{t-i}$$
$I_n = (X_1, X_2, \ldots, X_n)$라 하자. 먼저 $n$시점까지의 시계열 데이터가 주어졌을 때 $n+k$ 시점의 예측값 $X_{n+k, n}$은 다음과 같이 예측 오차 제곱을 최소화하는 것으로 한다.
$$E(X_{n+k}-f(I_n))^2$$
계산해보면 $X_{n+k,n}=E(X_{n+k}|I_n)$인 것을 알 수 있다. 따라서 우리는 $E(X_{n+k}|I_n)$을 계산해야한다.
$$\begin{align} E(X_{n+k}|I_n) &= c + \sum_{i=1}^p\phi_iE(X_{n+k-i}|I_n) + E(Z_{n+k}|I_n) - \sum_{i=1}^q\theta_iE(Z_{n-i}|I_n) \\ &= c+\sum_{i=1}^p\phi_iE(X_{n+k-i}|I_n) - \sum_{i=1}^q\theta_iE(Z_{n-i}|I_n) \end{align}$$
이때 조건부 기대값 $E(X_j|I_n)$은 $j\leq n$인 경우 관측값 $X_j$가 되고 $J>n$인 경우 $X_j$의 예측값으로 놓고 계산한다. 또한 $E(Z_j|I_n)$은 $j\leq n$인 경우 관측값 $Z_j$가 되고 $J>n$인 경우 0으로 계산한다.
실제로 예측할 때에는 $(c, \phi_1, \ldots, \phi_p, \theta_1, \ldots, \theta_q)$의 추정 값을 사용하며 $Z_j$의 경우 $Z_0 = Z_{-1} = \cdots = Z_{-q} = 0$으로 초기값을 설정하고 식 (1.1)의 관계식을 이용하여 계산한다.
3.2 예측 구간
먼저 예측 구간을 구하기 위해서는 $n+k$시점의 예측 오차 $e_{n+k,n} = X_{n+k} - X_{n+k, n}$과 예측 오차 분산 $\text{Var} (e_{n+k,n})$을 구해야 한다.
예측 오차의 분산은 식 (1.1)보다 ARMA의 MA 표현식을 이용하면 쉽게 계산할 수 있다. 먼저 예측 오차를 백색 잡음 $Z_t$로 표현해보자.
$$\begin{align}e_{n+k,n} &= X_{n+k} - X_{n+k, n} \\ &= \mu+Z_{n+k} + \psi_1 Z_{n+k-1} + \psi_2 Z_{n+k-2} + \cdots \\ &- (\mu + E(Z_{n+k}|I_n) + \psi_1 E(Z_{n+k-1})|I_n) + \psi_2 E(Z_{n+k-2}|I_n) +\cdots) \\ &= Z_{n+k}+\psi_1Z_{n+k-1}+\psi_2Z_{n+k-2}+\cdots+\psi_{k-1}Z_{n+1}\end{align}$$
백색 잡음은 서로 독립이기 때문에 분산 계산이 매우 쉬워진다.
$$\text{Var} (e_{n+k,n}) = (1+\psi_1^2+\cdots+\psi_{k-1}^2)\sigma_a^2$$
$n+k$ 시점의 값 $X_{n+k}$의 예측 구간은 다음과 같다.
$$X_{n+k,n}\pm z_{\alpha/2}\sqrt{\text{Var}(e_{n+k, n})}$$
$z_{\alpha/2}$는 표준 정규 분포의 확률 밀도 함수에서 오른쪽 면적이 $\alpha/2$가 되게 하는 값이다.
4. 파이썬 예제
이제 실제 데이터를 이용하여 ARMA 모형을 적합해보자. 여기서는 statsmodels 모듈을 이용하여 ARMA 모형을 적합한다. 필요한 모듈을 먼저 임포트 한다.
import pandas as pd
import numpy as np
import statsmodels.api as sm
import matplotlib.pyplot as plt
import warnings
warnings.filterwarnings('ignore')
from scipy.stats import norm
from statsmodels.tsa.stattools import acf
from statsmodels.tsa.arima.model import ARIMA4.1 데이터 탐색
먼저 아래의 데이터를 다운받자. 어느 화학 공정의 농도 측정 데이터라고 한다.
먼저 데이터를 읽어보자. 텍스트 파일로 되어 있어서 open을 이용하여 파일을 열고 각 줄을 읽는다. 데이터 값이 문자열이므로 이를 실수형으로 바꿔준다.
data = []
with open('./dataset/concentration.txt', 'r') as f:
lines = f.readlines()
for line in lines:
line = line.split(' ')
line[len(line)-1] = line[len(line)-1].replace('\n', '')
line = line[1:]
data += [x for x in line if x]
data = [float(x) for x in data]
이제 데이터가 어떻게 생겼는지 시각화해보자.
fig = plt.figure(figsize=(10,6))
fig.set_facecolor('white')
plt.plot(range(len(data)), data)
plt.show()
4.2 차수 결정
이제 ARMA 모형의 차수를 결정해보자. 여기서는 ESACF(Extended Sample Autocorrelation Function)과 AIC(Akaike Information Criterion)을 이용하여 차수를 결정해보겠다.
4.2.1 ESACF를 이용한 ARMA 차수 결정
아래 코드는 ESACF 테이블을 만들어 준다.
def esacf(data, ar_max=7, ma_max = 13, alpha=0.05, symbol=True):
sig = norm.ppf(1-alpha/2)
def lag_function(data, lag=1):
res = [np.nan]*lag + list(data[:-lag])
res = np.array(res)
return res
def ar_ols(data, ar_order):
depedent_data = np.array(data[ar_order:])
X = np.empty((0,ar_order))
for i in range(ar_order,len(data)):
temp_row = data[i-ar_order:i][::-1]
X = np.vstack([X, temp_row])
model = sm.OLS(depedent_data,X)
results = model.fit()
return results.params
def reupm(mat, ncol):
k = ncol-1
nrow = mat.shape[0]
# mat2[:] = np.nan
for i in range(k):
i1 = i+1
work = lag_function(mat[:,i])
work[0] = -1
temp = mat[:,i1] - (mat[i1,i1]/mat[i,i])*work
temp[i1] = 0
if i ==0:
mat2 = np.expand_dims(temp, axis=1)
else:
mat2 = np.column_stack((mat2, temp))
return mat2
ar_max += 1
ma_max += 1
nar = ar_max-1
nma = ma_max
ncov = nar + nma + 2
nrow = nar + nma + 1
ncol = nrow - 1
def ceascf(m, cov1, nar, ncol, count, ncov, z, zm):
result = [0]*(nar+1)
result[0] = cov1[ncov+count]
for i in range(nar):
temp = np.column_stack((z[i+1:], zm[i+1:,:i+1])).dot([1]+list(-mat2[:i+1, i]))
result[i+1] = acf(temp, nlags=count+1, fft=False)[count+1]
return result
z = data-np.mean(data)
for i in range(nar):
if i == 0:
zm = np.expand_dims(lag_function(z,i+1), axis=1)
else:
zm = np.column_stack((zm,lag_function(z,i+1)))
cov1 = acf(z, nlags=ncov, fft=False)
cov1 = np.array(list(cov1[1:][::-1]) + list(cov1))
ncov += 1
mat = np.zeros((nrow,ncol))
for i in range(ncol):
mat[:i+1,i] = ar_ols(z, ar_order=i+1)
for i in range(nma):
mat2 = reupm(mat, ncol)
ncol = ncol - 1
if i == 0:
eacfm = np.expand_dims(ceascf(mat2, cov1, nar, ncol, i, ncov, z, zm),axis=1)
else:
eacfm = np.column_stack((eacfm,ceascf(mat2, cov1, nar, ncol, i, ncov, z, zm)))
mat = mat2
if symbol:
work = len(z) - np.array(range(1, nar+2))+1
for j in range(nma):
work = work - 1
temp = np.abs(eacfm[:,j]) > sig/np.sqrt(work)
temp = np.array(['X' if t else 'O' for t in temp])
if j == 0:
symbol = np.expand_dims(temp, axis=1)
else:
symbol = np.column_stack((symbol, temp))
return pd.DataFrame(symbol)
else:
return pd.DataFrame(eacfm)
이제 주어진 데이터에 대하여 ESACF 테이블을 그려보자.
esacf(data, ar_max=5, ma_max = 8, alpha=0.01, symbol=True)
위 코드를 실행해보면 아래와 같은 테이블을 얻을 수 있다.
기호 'O'의 패턴을 보면 빨간색 선을 경계로 삼각 패턴이라는 것을 알 수 있다. 이 삼각형에서 좌측 상단 기호 'O'에 해당하는 MA와 AR 차수는 모두 1이므로 $ARMA(1,1)$이 적당하다는 것을 알 수 있다.
4.2.2 AIC를 이용하여 ARMA 차수 정하기
이번엔 AIC를 이용하여 ARMA 차수를 정해보자. 방법은 간단하다. 여러 차수 조합을 선정해놓고 가장 작은 AIC 값을 찾으면 된다. 나는 $0\leq p \leq 4, 0\leq q \leq 4$를 차수 후보로 정하였다. 각 후보 차수에 대하여 모형 적합은 ARIMA 클래스를 이용한다. 이 클래스는 필수적으로 모형을 적합시킬 데이터와 차수를 인자로 줘야 한다. 차수는 길이가 3인 튜플을 이용하는 데 이는 (AR 차수, 차분 차수, MA 차수)를 의미한다.
ARMA를 적합할 것이므로 가운데 값은 0으로 고정한다. 그런 다음 fit 메서드를 통해서 ARMA 모형을 적합하고 이를 model 변수에 저장한다. AIC 값은 aic 필드를 이용하여 얻을 수 있다. 이제 각 차수 조합별로 AIC 값을 그려본다.
ar_max = 4
ma_max = 4
aic_val = []
order_list = []
for ar in range(ar_max+1):
for ma in range(ma_max+1):
order = (ar, 0, ma)
order_list.append(str(order))
model = ARIMA(data, order=order).fit()
aic_val.append(model.aic)
fig = plt.figure(figsize=(10,6))
fig.set_facecolor('white')
plt.plot(range(len(aic_val)), aic_val)
plt.xticks(range(len(aic_val)), order_list, rotation=90)
plt.show()
위 그림을 확인해보니 $p=1, q=1$일 때 AIC 값이 가장 작았다. 해당 값을 프린트해보았다.
order_aic = dict(zip(order_list, aic_val))
final_order, final_aic = sorted(order_aic.items(), key=lambda x:x[1])[0]
print(f'Order : {final_order}, AIC : {final_aic}')
$p=1, q=1$일 때 AIC 값은 109.50 인 것을 알 수 있다. AIC 값을 이용한 결과도 $ARMA(1,1)$이 적절하다는 것을 알 수 있다.
4.2.3 ESACF와 AIC 결과가 일치하지 않는 경우
만약 ESACF와 AIC의 결과가 일치 하지 않는 경우에는 어떻게 해야할까?
다음의 해결방안을 생각할 수 있다.
1) 모수 절약 원칙에 따라 모수가 적은 모형을 선택한다.
2) 오차 제곱합 등의 방법을 이용하여 적합성을 측정해보고 적합도가 좋은 모형을 선택한다.
어느 방안이 더 좋다고 말할 수는 없고 상황에 맞게 사용하면 될 것 같다.
4.3 모형 추정
차수를 정했으니 최종 모형 $ARMA(p,q)$를 적합해보자. ARIMA 클래스를 사용하며 이에 대한 사용방법은 앞에서 살펴보았으니 생략한다. 적합 결과는 summary 메서드를 이용하여 알아볼 수 있다.
final_model = ARIMA(data, order=(1,0,1)).fit()
final_model.summary()
여기서 ar.~, ma.~ 에 해당하는 값은 각각 AR, MA 계수의 추정치이다. 위 결과를 토대로 최종 모형은 다음과 같음을 알 수 있다.
$$X_t = 17.0652+0.9086X_{t-1}+Z_t-0.5757Z_{t-1}$$
4.4 예측
이제 예측을 해보자. 난 미래 10 시점에 대한 예측값(out-of-sample)과 예측 구간을 시각화해보았다. 여기에는 우리가 적합한 모형을 사용하여 적합값(in-sample prediction)도 추정하였다. 적합값은 fittedvalues 필드를 이용하여 얻을 수 있다. 예측은 get_forecast 메서드를 이용한다. 미래 시점에 대한 예측값은 predicted_mean 필드를 이용하여 얻을 수 있고 예측 구간은 conf_int 메서드를 이용하여 얻을 수 있다.
fig = plt.figure(figsize=(10,6))
fig.set_facecolor('white')
fitted_values = final_model.fittedvalues ## in-sample prediction
## prediction
k = 10
forecast = final_model.get_forecast(k) ## out-of-sample prediction
predicted_values = forecast.predicted_mean
prediction_intervals = forecast.conf_int(alpha=0.5)
upper_limit = [x[1] for x in prediction_intervals]
lower_limit = [x[0] for x in prediction_intervals]
pred_x = range(len(data)+1, len(data)+k+1)
plt.plot(range(len(data)), data, label = 'data')
plt.plot(range(len(data)), fitted_values, label='fitted')
plt.plot(pred_x,predicted_values,color='red',label='prediction')
plt.plot(pred_x,lower_limit,color='red',linestyle='--',linewidth=1)
plt.plot(pred_x,upper_limit,color='red',linestyle='--',linewidth=1)
plt.fill_between(pred_x,lower_limit,upper_limit,color='red',alpha=0.2)
plt.legend()
plt.show()
실제로 적합이 잘된 것을 보려면 잔차를 분석해봐야 하는데 시각적으로 봤을 때 나쁘진 않은 것(?) 같다 ㅎㅎ
참고 자료
Hamilton, J.D. - Time Series Analysis
Tsay, R.S. - Analysis of Financial Time Series
댓글