Nonstationary process, ARIMA, SARIMA

일반적인 시계열 데이터는 stationary 하지 않습니다. 대부분 추세를 가지고 있고 따라서 stationary의 가정과 같이 평균이 상수로 일정하지 않습니다. 이전 포스트에서 알아봤듯이 추세를 제거해서 stationary 한 신호를 만들어 분석을 수행하게 됩니다. 

$X_t = X_{t-1} + a_t$의 데이터는 어떻게 될까요? 노이즈가 껴있긴 하지만 선형으로 증가하는 추세를 가지게 될겁니다. 그렇다면 이를 1차 difference ($\nabla X_t = X_t - X_{t-1}$) 를 수행하면 평균이 일정하게 되므로 stationary 하게 만들 수 있습니다.

추세가 quadratic하게 증가하는 시계열 데이터는 어떻게 다룰 수 있을까요? 이런 상황에는 2차 difference ($\nabla X_t^2$=$\nabla\nabla X_t$)를 수행하면 추세를 없애고 stationary 하게 만들 수 있습니다. 만약 시계열 데이터가 증가는 하는데, quadratic하게 세게 증가하지 않고 선형적으로 증가하지 않으면서 시간이 갈수록 variablity가 커지는 경우가 있습니다. 이러한 경우에는 box-cox transformation 을 통해 데이터에 $\log$를 취하게 됩니다.

우리는 이러한 데이터에 대해 difference를 수행한 이후에 ARMA 모델을 적용할 것입니다. 즉, ARIMA(p,d,q) 모델이 되는데, 가운데 d는 difference의 차원으로 d=1일 경우 1차 difference를 한 이후에 ARMA(p,q)를 적용하는 것입니다. (ARIMA의 I는 Integrated의 약자로 difference를 합쳤다고 생각하면 됩니다)

ARIMA

ARIMA(p,d,q) 를 수식으로 표현하면 다음과 같습니다.

$\Phi_p(B)\nabla^d X_t=\Psi_q(B)a_t + \mu$

$\Phi_p(B)$는 AR(p) 파트의 다항식이고 $\nabla^d$는 추세를 제거하기 위한 difference, $\Psi_q(B)$는 MA(q)의 다항식이니다. d는 일반적으로 AIC, BIC 등으로 찾아냅니다. 만약 차분 d를 지나치게 높인다면 필요없는 AR 파트도 생성하게 되어 불필요한 파라미터 예측 과정이 수행되게 됩니다. 

Prediction

ARIMA(1,1,1)의 예를 한 번 살펴보죠. ARIMA(1,1,1)은 $(1-\phi B)(1-B)X_t=(1+\psi B)a_t+\mu$로 표현할 수 있습니다. 이를 $X_t$에 대해 표현하면 $X_t=(\phi+1)X_{t-1}-\phi X_{t-2}+a_t+\psi a_{t-1}+\mu$가 될 것입니다. 

이를 이용해서 t+1 시점의 예측치를 추정하면, $E[X_{t+1}|x_1,...,x_t]$=$(\phi+1)x_t-\phi x_{t-1}+\tilde{a}_t+\mu$ 가 됩니다. t+2 시점의 예측치는 $(\phi+1)\hat{X}_{t+1}-\phi x_t+\mu$ 이 됩니다. ARIMA에서는 ARMA와는 달리 예측치가 특정 값으로 수렴하지 않습니다.

예측치의 에러 바운드는 어떻게 될까요? 마찬가지로 $X_t$를 노이즈의 causal한 선형결합으로 표현합니다. ($X_t$=$a_t+\phi_1 a_{t-1} + \phi_2 a_{t-2} + ...$) $\hat{X}_{t+1}, \hat{X}_{t+2}, ... $를 계산해보면 미래시점의 white noise는 평균을 0이라 생각하므로 ARMA와 마찬가지로 $Var(X_{t+1}-\hat{X}_{t+1})=\sigma_a^2$, $Var(X_{t+2}-\hat{X}_{t+2})=(1+\psi_1^2)\sigma_a^2$ 임을 알 수 있습니다. 하지만 ARIMA의 예측치 에러 바운드는 ARMA 와 달리 수렴하지 않고 시간이 갈수록 계속 커지게 됩니다.

간단한 $X_t=X_{t-1}+a_t+\mu$ 예에 대해 생각해 보겠습니다. $Var(X_t)$=$Var(X_{t-1})+\sigma_a^2$ 형태로 되어 예측치의 에러 바운드가 선형으로 증가하게 됩니다. 즉, $Var(X_{t+1}-\hat{X}_{t+1})=\sigma_a^2$가 되고 $Var(X_{t+2}-\hat{X}_{t+2})=2\sigma_a^2$ 의 형태로 예측치의 분산이 증가하게 됩니다. 따라서 ARIMA 또한 먼 미래의 값을 예측하는 데 적합하지 않습니다. 

 

Seasonal ARIMA

Seasonality는 짧은 주기성을 말하는 것으로 요일, 분기, 계절 등이 포함됩니다. SARIMA는 autocorrelation function이 0으로 감소하지 않고 주기성 (periodically)을 띄게 되었을 때 사용하게 됩니다. 밑의 그림에서는 seasonality가 12임을 알 수 있겠죠.

Seasonality를 갖고 있는 시계열 데이터에 ARIMA 를 적용해도 되지만 너무 많은 difference가 필요하므로 보통 적합하지 않습니다. Seasonality 6임을 찾기 위해서는 6번의 difference를 하여야 래그가 6인 시계열 데이터를 연결시킬 수 있고 이는 위에서 언급했듯이 불필요한 파라미터를 포함하게 되게 됩니다. 

단순한 $X_t=a_t - \theta a_{t-12}$를 생각해 보겠습니다. $Cov(X_t, X_{t-1})=0$이지만 $Cov(X_t, X_{t-12})$=$-\theta \sigma_a^2$ 임을 알 수 있습니다.

Seasonal MA(Q)는 MA$(Q)_S$로 표기되며 이는 수식으로 $X_t=a_t-\Theta_1 a_{t-s}-\Theta_2 a_{t-2s}-...$ 로 표현할 수 있습니다. Seasonal AR(P)는 AR$(P)_S$로 표기되며 $X_t=\Phi_1 X_{t-s}+\Phi_2 X_{t-2s} + ... \Phi_p X_{t-ps} + a_t$ 로 표현할 수 있습니다. Seasonal AR이 stationary 하기 위한 조건은 AR에서와 같습니다. 특히 AR$(1)_{12}$의 경우에는 $|\Phi|<1$ 이어야 stationary 조건을 만족하고 Yule-Walker equation은 $\rho_k=\Phi \rho_{k-12}$가 되어 $\rho_{12k}=\Phi^k$가 됩니다. 

이를 확장한 SARIMA(P,D,Q) 모델의 일반적인 수식은 다음과 같습니다.

$\Phi_P (B^S)\nabla_S^D X_t=\Phi_P (B^S)(1-B^S)^D X_t=\Theta_Q( B^S)a_t$

Seasonality를 고려하여 원래 시계열 데이터에서 Seasonality 부분을 제거한 나머지 부분에 대해서 ARIMA(p,d,q) 를 적용하면 됩니다. 즉, seasonality에 대한 ARIMA와 seasonality가 아닌 부분에 대한 ARIMA를 동시에 고려하겟다는 것이죠. 이를 합쳐서 더 일반적인 식으로 표현할 수 있습니다. 이를 $(p,d,q)\times (P,D,Q)_S$로 표기합니다.

$\Phi_P (B^S)\phi_p(B)\nabla_S^D \nabla^d X_t=\Phi_P (B^S)\phi_p(B)(1-B^S)^D (1-B)^d X_t=\Theta_Q(B^S)\theta_q(B)a_t$

예를 들어 $(0,1,1)\times (0,1,1)_{12}$은 $(1-B)(1-B^{12})X_t=(1-\theta B)(1-\Theta B^{12})a_t$ 모델이 됩니다.

일반적으로 시계열 데이터 모델링은

1) ACF, PACF 가 stationary 하도록 전처리를 합니다. (difference, log 등)
2) seasonality가 있다면 seasonality S를 정하고 seasonality ARIMA P,D,Q를 정하고
3) seasonality 제외한 래그에 대해 ACF, PACF를 살펴보고 이를 통해 ARIMA p,d,q 를 정합니다.

모델의 선택은 cross-validation으로 정확하게 정할 수 있지만 시간이 오래 걸리므로 보통 AIC, BIC 등으로 빠르게 판단합니다.