[Time Series] ARMA parameter의 MLE 추정

자 지금까지는 어떤 시계열이 주어졌을 때 그게 어떤 process이고 그 process를 규정하는 parameter들(\(\phi\)같은 애들)의 참값이 어떤값인지 아는 상태에서 covariance가 어떻게 되는지, 예측값은 어떻게 되는지 이런것들을 계산했습니다. 그런데, 실제 세상에서는 이런 parameter들의 참값을 알수는 없겠죠? 대개 주어진 시계열이 어떤 process를 따르는지 조차 모르는 상태에서 모델을 가정하고 parameter값을 통계적으로 추정하게 됩니다. 그리고 이렇게 추정하는 대표적인 방법이 Maximum Likelihood Estimation(MLE)입니다. MLE는 선형회귀분석쪽에서도 제가 다룬 적이 있었죠? 이번 포스팅에서는 ARMA에서의 MLE에 집중해서 살펴보겠습니다. 

 

ARMA 모델의 일반 형태는 다음과 같습니다:

\[ Y_t = c + \sum_{i=1}^p \phi_i Y_{t-i} + \epsilon_t + \sum_{j=1}^q \theta_j \epsilon_{t-j}, \]

• \(Y_t\): 시계열 데이터.
• \(\phi_i\): AR(autoregressive) 계수.
• \(\theta_j\): MA(moving average) 계수.
• \(\epsilon_t\): 평균이 0이고 분산이 \(\sigma^2\)인 white noise.

 

우리는 MLE를 활용해서 주어진 데이터 \((Y_1, Y_2, \ldots, Y_T)\)에 대해 파라미터 \(\theta = \{c, \phi_1, \ldots, \phi_p, \theta_1, \ldots, \theta_q, \sigma^2\}\)를 추정합니다. 보통 우리가 MLE를 사용할 때 그 시작점은 오차항의 분포로부터 시작합니다. 오차항이 어떤 분포를 따른다고 가정했으니, 오차항에 대해서 식을 정리하고 관측치에 대한 likelihood(우도)를 생각할 수 있는 것이죠. 그리고 이 likelihood 식안에 변수로서 parameter값들이 들어있으니 이제 우리가 평소에 자주 접하는 문제(함수를 최대화하는 값을 찾아라)로 바뀌는 겁니다.

ARMA에서도 MLE의 출발은 오차항. 오차항이 Gaussian white noise임을 가정합니다. 

 

AR(1) process의 MLE추정

AR(1)의 형태는 다음과 같습니다:

\[ Y_t = c + \phi Y_{t-1} + \epsilon_t, \quad \epsilon_t \sim N(0, \sigma^2) \]

일단 AR(1) Process의 첫번째 관측치를  \(Y_{1}\)이라고 할 때, 그 mean과 variance가 아래와 같았죠?

\[ E(Y_1) = \mu = \frac{c}{(1 - \phi)} \]

\[ \text{Var}(Y_1) = E((Y_1 - \mu)^2) = \frac{\sigma^2}{(1 - \phi^2)}. \]

AR(1) 식의 정의된 꼴을 보면 오차항이 Gaussian white noise일 경우 \(Y_{1}\)의 probability distribution도 Gaussian이어야 합니다. 그러면 첫번째 관측치에 대한 density는 아래와 같이 쓸 수 있겠죠? (회귀분석에서도 오차항을 Gaussian이라 가정하고 유사하게 likelihood를 계산했었습니다)

\[ f_{Y_1}(y_1; \theta) = \frac{1}{\sqrt{2 \pi \sigma^2 / (1 - \phi^2)}} \exp \left[ - \frac{(y_1 - \mu)^2}{2 \sigma^2 / (1 - \phi^2)} \right] \]

이제 그 다음 \(Y_{2}\) 는 어떨까요? 이제 \(Y_{1}\)이 관측된 상태에서 \(Y_{2}\)의 propbability distribution을 따져야 합니다. 조건부 확률분포죠. 이제 \(Y_{1}\)의 실제 관측된 실현(realization)값이 뭔지 알고 있고 그거를 \(y_{1}\)이라고 표기해보겠습니다. 이제는 \(y_{1}\)는 이미 결정된 값이기 때문에 확률변수가 아닌 확정된 상수(deterministic constant)로 취급할 수 있습니다. 그러면 \( Y_2 \)란 변수는 \(c + \phi y_1\)라는 상수에 Gaussian인 \( \epsilon_2\)를 더한 값이 기에 역시 Gaussian입니다. 

\[ (Y_2 | Y_1 = y_1) \sim N((c + \phi y_1), \sigma^2), \]

\[ f_{Y_2|Y_1}(y_2|y_1; \theta) = \frac{1}{\sqrt{2 \pi \sigma^2}} \exp \left[ - \frac{(y_2 - c - \phi y_1)^2}{2 \sigma^2} \right] \]

이후 관측되는  \(Y_{3}\), \(Y_{4}\), ... 모두 같은 방식으로 생각할 수 있고 전체 데이터셋의 joint likelihood는 다음과 같습니다:

\[ L(\theta) = f(Y_1) \prod_{t=2}^T f(Y_t | Y_{t-1}), \] \[ \log L(\theta) = -\frac{T}{2} \log(2\pi) - \frac{1}{2} \log\left(\frac{\sigma^2}{1-\phi^2}\right) - \sum_{t=2}^T \frac{(Y_t - c - \phi Y_{t-1})^2}{2\sigma^2} \]

 

MA(1) process의 MLE추정

MA(1)의 형태는 다음과 같습니다

\[ Y_t = \mu + \epsilon_t + \theta \epsilon_{t-1}, \quad \epsilon_t \sim N(0, \sigma^2) \]

AR(1) 프로세스와 구체적인 꼴이 다를 뿐, 발상은 비슷합니다. 만일 우리가 \(\epsilon_0=0\)이라고 가정한다면:

\[ (Y_{1}|\epsilon_0=0) \sim N(\mu, \sigma^2) \]

\(\epsilon_{t-1}\)이 관측된 상황에서는 역시 상수로 취급되므로

\[ (Y_{t}|\epsilon_{t-1}) \sim N(\mu + \theta \epsilon_{t-1}, \sigma^2) \]

따라서, MA(1)의 전체 데이터셋의 joint likelihood는 아래와 같이 나타낼 수 있습니다:

\[ \log L(\theta) = -\frac{T}{2} \log(2\pi) - \frac{T}{2} \log(\sigma^2) - \frac{1}{2} \sum_{t=1}^T \frac{\epsilon_t^2}{\sigma^2}. \]

 

ARMA 모델의 Likelihood는 이런식으로 계산할 수 있고, 우도를 최대화하는 과정에서는 어려가지 수치적(Numerical) 방법론이 동원될 수 있습니다. 사실 대부분의 경우에 이런 작업들은 통계소프트웨어 혹은 프로그래밍 package들이 해주는 작업이기에 거기까지 여기서 다루지는 않겠습니다.