이번 포스팅에서는 이동 평균 모형(Moving Average : MA)에 대해서 알아보고 파이썬으로 구현해보고자 한다. 또한 statsmodels를 이용하여 이동 평균 모형을 적합하는 방법도 알아보겠다.
1. 이동 평균 모형이란 무엇인가?
정상성(Stationary)을 갖는 시계열 데이터 $X_t$와 백색 잡음(White Noise) $Z_t$에 대하여 $q$차 이동 평균 모형 $MA(q)$는 다음과 같이 정의한다(백색 잡음의 정의는 여기를 참고하자).
$$X_t = c + Z_t + \theta_1Z_{t-1} + \cdots + \theta_qZ_{t-q}\tag{1-1}$$
이다. 위 식 양변에 기대값을 취해보자.
$$\begin{align} E(X_t) &= c + E(Z_t) + \theta_1E(Z_{t-1}) + \cdots + \theta_qE(Z_{t-q}) \\ &= c \end{align}$$
즉, 식 (1-1) 우변에서 $c$는 $X_t$의 기대값인 것을 알 수 있다. $X_t$는 정상성을 가지므로 $c = \mu$라 하겠다.
그런데 왜 식 (1-1)을 이동 평균이라고 부르는 걸까? '이동'이란 단어는 평균 $\mu$를 중심으로 해당 시계열이 왔다 갔다 이동한다는 의미이고 '평균'이라는 단어는 왔다 갔다 움직이는 정도를 백색 잡음 $\{Z_t\}$의 현재 값과 과거의 값들로 가중 합(weighted sum)을 했다는 것이다.
2. 언제 사용하는가?
이동 평균 모형은 자기 상관 함수(Autocorrelation Function : ACF)를 이용하여 사용 여부를 판단할 수 있다. 이는 다음과 같은 이론 정리에 기반한다.
정리)
시계열 $X_t$의 $l$차 ACF를 $\rho_l$이라 하자. 이때 $l>q$인 모든 $l$에 대하여 $\rho_l=0$이라면 $X_t$는 $MA(q)$를 따른다.
위 정리에 의하여 Sample ACF를 여러 차수에 대해서 그려보고 특정 차수 이후부터 Sample ACF 값이 크게 떨어진다면 이동 평균 모형을 사용하여 주어진 시계열을 모델링하는 것을 고려할 수 있다. 아래 그림을 살펴보자.
위 그림은 Sample ACF의 Correlogram을 그린 것이다. 이때 1차 이후부터 Sample ACF 값이 확 떨어진 것을 알 수 있다. 이 경우 $MA(1)$를 고려해 볼 수 있다.
자기 상관 함수(ACF)와 부분 자기 상관 함수(Partial ACF)에 대한 내용은 아래에 포스팅 해두었으니 참고하기 바란다.
[시계열 분석] 5. 자기 상관 함수(Autocorrelation Function : ACF)과 부분 자기 상관 함수(Partial Autocorrelatio
본 포스팅에서는 수식을 포함하고 있습니다. 티스토리 피드에서는 수식이 제대로 표시되지 않을 수 있으니 PC 웹 브라우저 또는 모바일 웹 브라우저에서 보시기 바랍니다. 이번 포스팅에서는 자
zephyrus1111.tistory.com
3. 모형 추정 방법
모형을 추정하기 위해선 먼저 차수 $q$를 정해야 하는데 이는 앞에서 언급했던 정리에 기반하여 ACF에 대한 Correlogram을 이용하여 정할 수 있다.
이동 평균 모형은 두 가지 방법으로 추정할 수 있다. 조건부 최대 우도 추정법(Conditional Maximum Likelihood Estimation)과 최대 우도 추정법(Exact Maximum Likelihood Estimation)이 그것이다. 이에 대해서 알아보자.
1. 조건부 최대 우도 추정법(Conditional Maximum Likelihood Estimation)
2. 정확한 최대 우도 추정법(Exact Maximum Likelihood Estimation)
1. 조건부 최대 우도 추정법
- $MA(1)$ -
먼저 다음의 $MA(1)$의 모형을 생각해보자.
$$X_t = \mu + Z_t+\theta Z_{t-1} \tag{3.1.1}$$
이때 $Z_t \sim \text{i.i.d.} N(0, \sigma^2) $이고 $\boldsymbol \theta = (\mu, \theta, \sigma^2)^t$를 추정해야 할 모수 벡터라고 하자. 여기에서 $Z_{t-1} = z_{t-1}$로 주어졌다면 우리는 다음과 같음을 알 수 있다.
$$X_t|Z_{t-1} = z_{t-1} \sim N(\mu+\theta z_{t-1}, \sigma^2)$$
$z_0=0$이라고 가정하자(이 가정이 조건부 최대 우도 추정법의 핵심이다). 만약 $x_1$의 값이 주어졌다면 식 (3.1.1)에 의하여 $z_1$의 값을 알 수 있다.
$$z_1 = x_1-\mu$$
이때 $X_1=x_1, Z_0=0$으로 주어진 경우 $X_2$의 조건부 확률분포는 다음과 같다.
$$\begin{align} f_{X_2|X_1, Z_0}(x_2|x_1, z_0; \boldsymbol{\theta}) &= \frac{1}{\sqrt{2\pi\sigma^2}}\exp \left [\frac{-(x_2-\mu-\theta z_1)^2}{2\sigma^2} \right ] \\ &= \frac{1}{\sqrt{2\pi\sigma^2}}\exp \left [-\frac{z_2^2}{2\sigma^2} \right ] \end{align}$$
다음으로 $Z_1=z_1$과 $X_2=x_2$가 주어진 경우 마찬가지로 식 (3.1.1)로부터 $z_2$의 값을 알 수 있다.
$$z_2 = x_2-\mu - \theta z_1$$
이러한 방식으로 계속 진행하면 다음을 알 수 있다.
$$z_t = x_t-\mu-\theta z_{t-1}\tag{3.1.2}$$
그리고 조건부 확률 밀도 함수 또한 다음과 같음을 알 수 있다.
$$\begin{align} f_{X_t|X_{t-1}, \ldots, X_1, Z_0}(x_t|x_{t-1}, \ldots, x_1, z_0=0; \boldsymbol{\theta}) &= f_{X_t|Z_{t-1}}(x_t|z_{t-1};\boldsymbol{\theta}) \\ &= \frac{1}{\sqrt{2\pi\sigma^2}}\exp \left [ -\frac{z_t^2}{2\sigma^2} \right ] \end{align}$$
이제 시계열 데이터 $X_T , X_{T-1}, \ldots, X_1$의 우도 함수는 다음과 같이 계산할 수 있다.
$$\begin{align} f_{X_T, X_{T-1}, \ldots, X_1|Z_0}(x_t, x_{t-1}, \ldots, x_1|z_0=0; \boldsymbol{\theta})\\= f_{X_1|Z_0}(x_1|z_0 ; \boldsymbol{\theta}) \prod_{t=2}^Tf_{X_t|X_{t-1}, \ldots, X_1, Z_0}(x_t|x_{t-1}, \ldots, x_1, z_0=0;\boldsymbol{\theta} )\end{align}$$
이제 위 식 양변에 로그를 취하고 정리해주면 조건부 로그 우도(conditional log likelihood) 함수를 다음과 같이 얻을 수 있다.
$$L(\boldsymbol \theta) = -\frac{T}{2}\log (2\pi) - \frac{T}{2}\log (\sigma^2) - \sum_{t=1}^T\frac{z_t^2}{2\sigma^2}\tag{3.1.3}$$
식 (3.1.3)는 모수에 대해서 쉽게 풀리지 않는다. 따라서 수치적인 방법을 이용해서 구해야 한다.
식 (3.1.2)를 이용하면 임의의 $z_0$ 값에 대하여
$$z_t = (y_t-\mu) - \theta (x_{t-1}-\mu) + \theta^2(x_{t-2}-\mu) - \cdots +(-1)^{t-1}\theta^{t-1}(x_1-\mu)+ (-1)^t\theta^tz_0$$
만약 $|\theta| < 1$이라면 $z_0$에 대한 영향은 지수적(exponential)으로 줄어들게 되며 조건부 우도 함수는 정확한 우도 함수(Exact Likelihood Function)를 잘 근사하게 된다. 하지만 $|\theta|>1$이라면 조건부 우도 함수를 이용한 모형 추정방법은 그리 좋지 못하다고 한다. 따라서 최종적으로 추정된 $\hat{\theta}$값이 1보다 크다면 $\hat{\theta}$의 역수를 초기값으로 두고 추정해야 한다고 한다.
- $MA(q)$ -
지금까지 $MA(1)$ 모형에 대한 조건부 우도 함수를 계산해보았는데 $MA(q)$ 모형에 대한 조건부 우도 함수를 계산해보자. 먼저 다음의 $MA(q)$ 모형을 고려해보자.
$$X_t = \mu + Z_t + \theta_1 Z_{t-1} + \cdots + \theta_q Z_{t-q}$$
이때 다음을 가정한다.
$$z_0=z_{-1}=\cdots=z_{-q+1}=0$$
이러한 가정을 이용하면 우리는 $MA(1)$의 조건부 우도 함수를 유도했던 것과 비슷한 점화식(3.1.2)을 계산할 수 있다.
$$z_t = x_t-\mu - \theta_1 z_{t-1} - \cdots - \theta_q z_{t-q}\tag{3.1.*}$$
$MA(1)$에서 비슷한 방법으로 조건부 로그 우도 함수를 다음과 같이 계산할 수 있다.
$$L(\boldsymbol \theta) = -\frac{T}{2}\log (2\pi) - \frac{T}{2}\log (\sigma^2) - \sum_{t=1}^T\frac{z_t^2}{2\sigma^2}\tag{3.1.4}$$
여기서 $\boldsymbol \theta = (\mu, \theta_1, \ldots, \theta_q, \sigma^2)^t$이다. 식 (3.1.4) 역시 모수에 대하여 꽤 복잡한 비선형함수이다. 따라서 수치적인 방법을 이용하여 추정해야 한다. 이에 대한 방법은 아래에서 다루겠다.
2. 정확한 최대 우도 추정법
- $MA(1)$ -
조건부 우도 추정법에서 $z_0=0$으로 가정했지만 이러한 가정을 하지 않고 최대 우도 추정을 할 수 있다. 이를 위해서 $X_T, X_{T-1}, \ldots, X_1$의 우도 함수를 구해야 한다. 먼저 $MA(1)$부터 살펴보자. 여기에서도 $Z_t \sim \text{i.i.d. } N(0, \sigma^2)$이라 가정한다. 이 경우 $\mathbf{X} = (X_1, X_2, \ldots, X_T)^t$은 $T$차원 정규 분포를 따르며 확률 밀도 함수는 다음과 같다.
$$f_{\mathbf{X}}(\mathbf{x};\boldsymbol \theta) = (2\pi)^{-T/2}|\boldsymbol \Omega|^{-1/2}\exp \left [-\frac{1}{2} (\mathbf{x}-\boldsymbol \mu)^t\boldsymbol \Omega^{-1} (\mathbf{x} - \boldsymbol \mu)\right] \tag{3.2.1}$$
여기서 $\boldsymbol \mu = (\mu, \mu, \ldots, \mu)^t$, $T\times 1$ 벡터이고 $\boldsymbol \Omega$는 $T\times T$ 행렬로 아래와 같다.
$$\boldsymbol \Omega =
\sigma^2\begin{pmatrix}
(1+\theta^2) & \theta & 0 & \cdots & 0 \\
\theta & (1+\theta^2) & \theta & \cdots & 0 \\ 0 & \theta & (1+\theta^2) &\cdots & 0 \\ \vdots & \vdots &\vdots & \ddots & \vdots \\
0 & 0 & 0& \cdots & (1+\theta^2)
\end{pmatrix}$$
계산의 편의성을 위해 $\boldsymbol \Omega$를 다음과 같이 분해한다.
$$\boldsymbol \Omega = ADA^t$$
여기서
$$A =
\begin{pmatrix}
1& 0 & \cdots &0 & 0 \\
\frac{\theta}{1+\theta^2}& 1& \cdots & 0 & 0 \\
\vdots & \vdots & \ddots & \vdots &\vdots \\
0 & 0 & \cdots & \frac{\theta (1+\theta^2 + \cdots + \theta^{2(T-1)})}{1+\theta^2 + \cdots + \theta^{2(T-2)}}& 1
\end{pmatrix}$$
$$D =
\begin{pmatrix}
1+\theta^2& 0 & \cdots & 0 \\
0 & \frac{1+\theta^2 + \theta^4}{1+\theta^2} & \cdots & 0 \\
\vdots & \vdots & \ddots &\vdots \\
0 & 0 & \cdots & \frac{1+\theta^2 + \cdots + \theta^{2T}}{1+\theta^2 + \cdots + \theta^{2(T-1)}}
\end{pmatrix}$$
이다.
이를 이용하여 식 (3.2.1)을 다시 쓰면 다음과 같다.
$$f_{\mathbf{X}}(\mathbf{x};\boldsymbol \theta) = (2\pi)^{-T/2}|ADA^t|^{-1/2}\exp \left [-\frac{1}{2} (\mathbf{x}-\boldsymbol \mu)^t (A^t)^{-1}D^{-1}A^{-1} (\mathbf{x} - \boldsymbol \mu)\right]$$
$A$의 대각 원소는 모두 1이므로 $|A| = |A^t| = 1$이다. 따라서 $|ADA^t| = |A||D||A^t| = |D|$이다. 그리고 $\bar{\mathbf{x}} = A^{-1}(\mathbf{x}-\boldsymbol \mu)$라 하면 식 (3.2.1)은
$$f_{\mathbf{X}}(\mathbf{x};\boldsymbol \theta) = (2\pi)^{-T/2}|D|^{-1/2}\exp \left [-\frac{1}{2} \bar{\mathbf{x}}^tD^{-1} \bar{\mathbf{x}}\right]$$
$A\bar{\mathbf{x}} = \mathbf{y}-\boldsymbol \mu$이므로 이를 통하여 $\bar{\mathbf{x}}$를 계산할 수 있다. $D$의 $t$번째 대각 원소를 $d_{tt}$라 하면 로그 우도 함수는 다음과 같이 계산할 수 있다.
$$L(\boldsymbol \theta) = -\frac{T}{2} \log (2\pi) - \frac{1}{2} \sum_{t=1}^T\log (d_{tt}) - \frac{1}{2}\sum_{t=1}^T\frac{\bar{\mathbf{x}}^2}{d_{tt}}\tag{3.2.2}$$
주어진 $\mu, \theta, \sigma^2$에 대하여 $\bar{\mathbf{x}}$는 계산할 수 있다. 식 (3.2.2)를 최대화하는 문제도 수치적인 방법을 이용하여 구한다. 조건부 우도 추정과는 다르게 $\theta$의 값과 상관없이 최대 우도 추정량을 구할 수 있다.
- $MA(q)$ -
다음으로 $MA(q)$에 대해서 살펴보자. 앞서 살펴본 것과 같이 $(X_1, X_2, \ldots, X_T)$는 $T$ 차원 정규분포를 따르며 확률 밀도 함수는 (3.2.1)과 같다. 다만 몇 가지 $\boldsymbol \Omega$가 달라진다.
여기서 $\gamma_k$는 $k$차 자기 공분산 함수로서 다음과 같이 계산할 수 있다.
$$\gamma_k = \sigma^2\rho_k = \begin{cases} \sigma^2(\theta_k +\theta_{k+1}\theta_1+\cdots+\theta_q\theta_{q-k}) & \text{ for } k=1, \ldots, q \\ 0 & \text{ for } k>q \end{cases}$$
$MA(1)$에서 살펴본 것과 같이 $\boldsymbol \Omega$를 분해하고 (3.2.2)의 형태로 만들 수 있다. 모형 추정은 수치적인 방법을 통하여 계산한다.
3. 모형 추정 알고리즘
로그 우도 함수를 구했으니 이제 이를 최대화하는 모수 $\boldsymbol \theta$를 찾아야 한다. 여기에서는 Newton-Rhapson 방법을 이용하려고 한다. 여기에서 소개하는 알고리즘은 조건부 최대 우도 추정이다. 정확한 최대 우도 추정은 구현해봤는데 잘 안돼서 잘되면 정리해보겠다. 로그 우도 함수를 $L(\boldsymbol \theta)$라 하자. 그리고
$$g(\boldsymbol \theta) = \frac{\partial L}{\partial \boldsymbol \theta} \\ H(\boldsymbol \theta) = -\frac{\partial ^2 L}{\partial \boldsymbol \theta \partial \boldsymbol \theta^t} $$
라 하자.
반복 단계로 업데이트하는 모수는 $\theta_1, \ldots, \theta_q, \mu$이며 $\tilde{\boldsymbol \theta} = \theta_1, \ldots, \theta_q, \mu$라 하자. 이 경우 모형 추정 알고리즘은 다음과 같다.
알고리즘
1 단계) 초기값 $\boldsymbol \theta^{(0)}$과 오차범위 $\epsilon$를 설정한다.
2 단계) $m$번째 스텝에서 추정치$\tilde{\boldsymbol \theta}^{(m)}$과 $\hat{\sigma}^2_m$주어졌을 때 $m+1$번째 추정량은 다음과 같이 업데이트한다.
$$\tilde{\boldsymbol \theta}^{(m+1)} = \tilde{\boldsymbol \theta}^{(m)} + sH(\tilde{\boldsymbol \theta}^{(m)})^{-1}g(\tilde{\boldsymbol \theta}^{(m)})$$
이때 $s>0$는 여러 가지 후보 값들 중에서 $g$값을 최대로 하는 값을 선택한다.
3 단계) $\tilde{\boldsymbol \theta}^{(m+1)}$을 이용하여 $\hat{\sigma}^2_{m+1}$ 추정치를 계산한다. 계산식은 다음과 같다.
$$\hat{\sigma}^2_{m+1} = \frac{1}{T}\sum_{t=1}^Th(\tilde{\boldsymbol \theta}^{(m+1)})$$
이때 식 (3.1.4)에서 $h(\tilde{\boldsymbol \theta})=z^2_t$이다.
4 단계) $\| \tilde{\boldsymbol \theta}^{(m+1)} -\tilde{\boldsymbol \theta}^{(m)} \| < \epsilon$이면 알고리즘을 종료하고 아닌 경우 2 단계)로 넘어간다.
4. 파이썬 구현
이제 이동 평균 모형을 적합하는 코드를 만들어보려고 한다. 여기에서는 조건부 최대 우도 추정 방법으로 모형을 적합한다. 먼저 이번 포스팅에서 필요한 모듈을 임포트 한다.
import pandas as pd
import numpy as np
np.set_printoptions(suppress=True) ## scientific notation mode off
import matplotlib.pyplot as plt
import statsmodels.api as sm
from scipy.stats import norm
from scipy.linalg import lu
from tqdm import tqdm
from statsmodels.tsa.arima.model import ARIMA
이동 평균 모형을 적합하는데 필요한 클래스를 정의한다. 이 클래스는 차수 $q$와 데이터를 인자로 준다.
class movingAverageModel():
def __init__(self, order, x):
self.x = x
self.order = order
self.theta = None
self.mu = None
self.sigma2 = None
다음으로 모형을 실질적으로 적합하는 fit 함수를 정의한다. fit 함수는 method 인자에 따라서 조건부 최대 우도 추정 또는 정확 최대 우도 추정법을 이용하여 파라미터를 추정한다. 원래는 조건부뿐만 아니라 정확한 최대 우도 추정법까지 모두 아우르는 fit 함수를 구현하고 싶었으나 정확 최대 우도 추정은 원하는 결과가 잘 안 나와서 생략한다. ㅠ.ㅠ
class movingAverageModel():
'''생략'''
def fit(self, step_size=0.01, nr_epsilon=0.001, max_cnt=100, epsilon=0.0001, method='conditional'):
assert method in ['conditional', 'exact']
x = self.x
if method=='conditional':
self._fit_cond(x, step_size, nr_epsilon, max_cnt, epsilon)
else:
pass
## to do
# self._fit_exact()
다음으로 조건부 최대 우도 추정법으로 파라미터를 추정하는 코드이다. 최대 반복 횟수 max_cnt를 설정하고 Newton-Rhapson 방법을 이용하여 파라미터를 업데이트한다. 이때 파라미터 현재값과 다음값의 차이가 오차범위 nr_epsilon 내에 있다면 알고리즘을 종료한다.
class movingAverageModel():
'''생략'''
def _fit_cond(self, x, step_size=0.1, nr_epsilon=0.0001, max_cnt=20, epsilon=0.0001):
## initialization
parameter = [0]*order + [np.mean(x)]
parameter = np.array(parameter)
sigma2 = np.var(x)
cnt = 1
while cnt <= max_cnt:
print(cnt, ' ', parameter)
first_deriv = self.first_deriv_likelihood(x, parameter, sigma2, order, epsilon=epsilon)
hess_matrix = self.get_hessian(x, parameter, sigma2, order, epsilon=epsilon)
next_parameter = parameter + step_size*np.linalg.inv(hess_matrix).dot(first_deriv)
next_sigma2 = np.mean(np.square(self.get_z(x,next_parameter, order)))
if np.linalg.norm(next_parameter-parameter) < nr_epsilon:
break
cnt += 1
parameter = next_parameter
sigma2 = next_sigma2
self.theta = parameter[:-1]
self.mu = parameter[-1]
self.sigma2 = sigma2
아래 함수들은 조건부 최대 우도 추정을 위해 필요한 함수들이다. 차례대로 (3.1.4)에서 $z$를 계산하는 함수, 로그 우도 함수의 근사(approximation) 1차 미분 벡터, 헤시안 행렬의 근사값을 계산하는 함수이다.
class movingAverageModel():
'''생략'''
def get_z(self, x, parameter, order):
z = [0]*len(x)
theta = parameter[:-1]
mu = parameter[-1]
z_init = [0]*order
for i in range(len(x)):
z_total = z_init + z
target_z = z_total[i:i+order]
target_z = target_z[::-1]
z[i] = x[i] - mu - np.dot(theta, target_z)
return np.array(z)
def first_deriv_likelihood(self, x, parameter, sigma2, order, epsilon=0.0001):
deriv_vector = []
p = len(parameter)
for idx in range(p):
e = np.zeros(p)
e[idx] = epsilon
z1 = self.get_z(x,parameter+e, order)
z2 = self.get_z(x,parameter, order)
deriv = (np.square(z1) - np.square(z2))/(2*sigma2*epsilon)
temp = -np.sum(deriv)
deriv_vector.append(temp)
deriv_vector = np.array(deriv_vector)
return deriv_vector
def get_hessian(self, x,parameter, sigma2, order, epsilon=0.0001):
p = len(parameter)
hess_matrix = np.zeros((p, p))
for i in range(p):
for j in range(i+1):
e_i = np.zeros(p)
e_i[i] = epsilon
e_j = np.zeros(p)
e_j[j] = epsilon
z1 = self.get_z(x,parameter+e_i+e_j, order)
z2 = self.get_z(x,parameter+e_i-e_j, order)
z3 = self.get_z(x,parameter-e_i+e_j, order)
z4 = self.get_z(x,parameter-e_i-e_j, order)
temp_vector = np.square(z1) - np.square(z2) - (np.square(z3)-np.square(z4))
temp_vector = 0.25*temp_vector/(2*sigma2*np.square(epsilon))
temp = np.sum(temp_vector)
hess_matrix[i][j] = temp
hess_matrix = hess_matrix + hess_matrix.T - np.diag(np.diagonal(hess_matrix))
return hess_matrix
4. 예측(Forecasting)
이제 $MA(q)$ 모형을 이용하여 미래 시점의 값을 예측해보자.
1. 예측값
먼저 $I_n = (X_1, X_2, \ldots, X_n)$이라 하자. $n$ 시점까지의 시계열 데이터가 주어졌을 때 $n+k$ 시점의 예측값 $X_{n+k, n}$은 다음과 같이 예측 오차 제곱을 최소화하는 $f$로 한다.
$$E(X_{n+k}-f(I_n))^2$$
계산해보면 $E(X_{n+k}|I_n)$인 것을 알 수 있다. 즉, $X_{n+k, n} = E(X_{n+k}|I_n)$ 이다.
먼저 $X_{n+1, n}$을 계산해보자.
$$\begin{align} X_{n+1, n} &= E(X_{n+1}|I_n) \\ &= \mu + E(Z_{n+1}|I_n) + \theta_1 E(Z_n|I_n) + \cdots + \theta_q E(Z_{n-q+1}|I_n) \\ &= \mu + \theta_1 Z_{n} + \cdots + \theta_q Z_{n-q+1} \end{align}$$
두 번째 등식에서 $E(Z_{n+1}|I_n) = 0$이고 세 번째 등식에서 $k=n-p+1, \ldots, n$에 대하여 $E(X_k|I_n) = X_k$이다. 이를 바탕으로 $X_{n+k, n}$을 계산하면 다음과 같다.
$$X_{n+k, n}=\begin{cases} \mu +\theta_kZ_{n} + \cdots + \theta_qZ_{n-q+k} & \text{ for } k=1, \ldots, q \\ \mu & \text{ for } k>q \end{cases}$$
실제로 예측할 때에는 $\mu, \theta_1, \ldots, \theta_q$ 대신 추정치를 사용하며 $Z_t$의 경우 (3.1.*)의 관계식을 이용하여 계산한다.
2. 예측 구간
예측 구간을 구하기 위하여 $n+k$ 시점의 예측 오차 분산 $var (e_{n+k, n})$을 구해야 한다. 여기서 $e_{n+k, n} = X_{n+k} - X_{n+k, n} = X_{n+k} - E(X_{n+k}|I_n)$이다. 먼저 $var (e_{n+1, n})$을 계산해보자.
$$\begin{align} var (e_{n+1, n}) &= E(X_{n+1}-X_{n+1, n})^2 \\ &= E(\mu + Z_{n+1} +\theta_1 Z_n +\cdots +\theta_qZ_{n-q+1} - \mu - \theta_1 Z_n + \cdots + \theta_q Z_{n-q+1})^2 \\ &= E(Z_{n+1}^2) = \sigma^2 \end{align}$$
위에서 살펴본 것을 바탕으로 $var(e_{n+k, n})$은 다음과 같다.
$$var(e_{n+k, n}) = \begin{cases} (1+\theta_1^2 + \cdots + \theta_{k-1}^2)\sigma^2 & \text{ for } k=1, \ldots, q \\ (1+\theta_1^2+\cdots +\theta_q^2)\sigma^2 & \text{ for } k>q \end{cases}$$
예측 구간을 계산할 때에는 정규분포의 꼬리 값을 이용한다. 따라서 $n+k$시점에서의 ($1-\alpha$)% 예측 구간은 다음과 같다.
$$E(X_{n+k}|I_n) \pm z_{\alpha/2}\sqrt{var (e_{n+k, n})}$$
내용
3. 파이썬 코드
이번엔 예측을 수행하는 코드를 만들어보려고 한다. 나는 기본적으로 예측하려는 최대 시점(max_step)을 정한 뒤 1 시점부터 최대 시점 까지의 예측값을 계산한다. 예측 오차 분산과 예측 구간도 마찬가지이다.
아래 코드는 각각 1) 1시점 부터 최대 시점까지의 예측값과 개별 예측값, 2) 최대 시점까지의 예측 오차 분산과 개별 시점 예측 오차 분산 그리고 3) 최대 시점까지의 예측 구간과 개별 시점 예측 구간을 계산하는 함수이다.
class movingAverageModel():
'''생략'''
def predict(self, max_step):
return [self._predict(step) for step in range(1,max_step+1)]
def _predict(self, step):
theta = self.theta
mu = self.mu
q = len(theta)
if step<=q:
order = self.order
parameter = list(theta) + [mu]
x = self.x
z = self.get_z(x, parameter, order)
target_z = z[len(z)-q+step-1:]
target_z = target_z[::-1]
future_val = mu + theta[step-1:].dot(target_z)
else:
future_val = mu
return future_val
def prediction_err_var(self, max_step):
return [self._prediction_err_var(step) for step in range(1, max_step+1)]
def _prediction_err_var(self, step):
theta = self.theta
sigma2 = self.sigma2
q = len(theta)
theta_total = [1] + list(theta)
theta_total = np.array(theta_total)
if step <= q:
pev = np.sum(np.square(theta_total[:step]))*sigma2
else:
pev = np.sum(np.square(theta_total))*sigma2
return pev
def prediction_interval(self, max_step, alpha = 0.05):
return [self._prediction_interval(step, alpha) for step in range(1, max_step+1)]
def _prediction_interval(self, step, alpha = 0.05):
future_val = self._predict(step)
z_alpha = norm.ppf(1-alpha*0.5)
pev = self._prediction_err_var(step)
upper_limit = future_val+z_alpha*np.sqrt(pev)
lower_limit = future_val-z_alpha*np.sqrt(pev)
return (lower_limit, upper_limit)
이제 구현은 끝났다. 실제 데이터를 이용하여 모형을 적합시켜 보자~
5. 예제
Ruey S. Tsay 책 'Analysis of Finantial Time Series'에서 다룬 IBM 회사의 1926년 1월부터 2003년 12까지 Equal-weighted 월간 이익을 사용하려고 한다. 먼저 전체 데이터를 다운 받는다. 여기에는 월간 이익, Value-weighted 월간 이익, 그리고 S&P 500 값이 포함되어 있다.
1. 데이터 불러오기 및 시각화
데이터를 다운 받았으면 불러오자. 이 데이터는 정확하게 탭으로 분리된 것이 아닌 띄어쓰기 여러 번으로 분리되어 있다. 이때에는 read_csv 함수의 sep 인자에 '\s+'를 넣어주면 된다. 그리고 날짜 칼럼을 스트링에서 날짜 형식을 변환해주었다.
names=['Date','IBM','VW','EW','SP']
df = pd.read_csv('./dataset/monthly_simple_return_of_IBM.txt', sep='\s+',header=None, names=names)
df['Date'] = pd.to_datetime(df['Date'], format='%Y%m%d')
df.head()
데이터가 어떻게 생겼는지 확인해보자.
fig = plt.figure(figsize=(8,8))
fig.set_facecolor('white')
plt.plot(df['Date'], df['EW'])
plt.show()
모형의 차수를 정하기 위하여 ACF Correlogram을 그려본다. statsmodels 라이브러리에서는 graphics.tsa.plot_acf 함수를 이용하여 그릴 수 있다.
fig = plt.figure(figsize=(8,8))
fig.set_facecolor('white')
ax = fig.add_subplot()
sm.graphics.tsa.plot_acf(df['EW'], ax=ax)
plt.show()
위 그림에서 보면 1차, 3차, 9차 ACF가 유의하다고 나온다(10차 이후는 무시). 나는 그중에서 3차를 선택하였고 $MA(3)$차 모형을 적합하려고 한다.
2. 모형 적합
2.1 구현한 모델 이용
미래 1~5 시점까지는 예측값과 비교해보기 위하여 이는 빼놓은 나머지 데이터를 이용하여 $MA(3)$ 모형을 적합한다.
order = 3
x = df['EW']
train_x = x.iloc[:-5]
ma = movingAverageModel(order,train_x)
ma.fit(step_size=1)
3번째 스텝 만에 모형 추정이 완료되었다. 결과를 출력해보았다.
print('MA parameter : ', ma.theta)
print('Mean : ', ma.mu)
print('WN Variance : ', ma.sigma2)
위 결과를 통해 최종 모형은 다음과 같음을 알 수 있다.
$$X_t = 0.0129 + Z_t + 0.2000Z_{t-1} + 0.0313Z_{t-2} - 0.0936 Z_{t-3}$$
2.2 statsmodels 이용
이번엔 statsmodels를 이용하여 $MA(3)$을 적합해보자. statsmodels에서는 ARIMA라는 클래스를 제공하여 이동 평균 모형뿐만 아니라 다른 여러 가지 모형들도 적합할 수 있게 해 준다. 우리는 3차 이동 평균 모형을 적합할 것이므로 order인자에 (0,0,3)을 넣어준다.
ma_model = ARIMA(train_x, order=(0,0,3)).fit()
이제 summary 함수를 이용하여 결과를 확인해보자.
ma_model.summary()
위 결과를 보면 const 부분이 $\mu$의 추정치이고 ma.L1~L3 까지는 각각 $\theta_1$~$\theta_3$의 추정치이며 마지막으로 sigma2는 백색 잡음의 분산 $\sigma^2$의 추정치이다. 내가 구현한 것과 비교해보니 거의 비슷했다. 아마도 미분 계산 방식이나 알고리즘 구현에 의한 차이인 것 같다.
3. 예측
3.1 구현한 모델 이용
먼저 미래 5 시점까지의 예측값과 예측 구간을 계산했다.
max_step = 5
predict_val_list = ma.predict(max_step)
predict_interval_list = ma.prediction_interval(max_step)
upper_limits = [x[1] for x in predict_interval_list]
lower_limits = [x[0] for x in predict_interval_list]
다음으로 이를 시각화하여 예측값과 실제값을 비교해보고자 한다.
fig = plt.figure(figsize=(6,6))
fig.set_facecolor('white')
plt.plot(df['Date'].iloc[-10:], df['EW'].iloc[-10:], label='Data', color='b')
plt.plot(df['Date'].iloc[-5:], predict_val_list, label='Prediction', color='g')
## prediction, prediction interval
plt.plot(df['Date'].iloc[-5:],upper_limits,color='g',linestyle='--',linewidth=1)
plt.plot(df['Date'].iloc[-5:],lower_limits,color='g',linestyle='--',linewidth=1)
plt.fill_between(df['Date'].iloc[-5:],lower_limits,upper_limits,color='purple',alpha=0.2)
plt.xticks(rotation=90)
plt.plot()
plt.legend()
plt.show()
예측이 많이 빗나가긴 했다 ㅎㅎ;; 그래도 예측 구간 안에는 다 들어온다.
3.2 statsmodels 이용
이번엔 statsmodels를 이용해보자. forecast 함수를 사용하면 미래 시점의 값을 예측할 수 있고 params.sigma2 필드를 이용하여 백색 잡음 추정치를 가져올 수 있다.
alpha = 0.05
predict_val_list = ma_model.forecast(max_step)
std = np.sqrt(ma_model.params.sigma2)
upper_limits = [x+norm.ppf(1-alpha*0.5)*std for x in predict_val_list]
lower_limits = [x-norm.ppf(1-alpha*0.5)*std for x in predict_val_list]
마찬가지로 시각화해보았다.
fig = plt.figure(figsize=(6,6))
fig.set_facecolor('white')
plt.plot(df['Date'].iloc[-10:], df['EW'].iloc[-10:], label='Data', color='b')
plt.plot(df['Date'].iloc[-5:], predict_val_list, label='Prediction', color='g')
## prediction, prediction interval
plt.plot(df['Date'].iloc[-5:],upper_limits,color='g',linestyle='--',linewidth=1)
plt.plot(df['Date'].iloc[-5:],lower_limits,color='g',linestyle='--',linewidth=1)
plt.fill_between(df['Date'].iloc[-5:],lower_limits,upper_limits,color='purple',alpha=0.2)
plt.xticks(rotation=90)
plt.plot()
plt.legend()
plt.show()
앞에서 내가 구현한 것과 비교해보기 위해 앞의 결과도 같이 넣었다. 왼쪽이 statsmodels, 오른쪽에 내가 구현한 모델의 예측 결과이다.
약간의 차이는 있지만 거의 동일하다.
6. 마치며
이번 포스팅에서는 이동 평균 모형에 대해서 알아보았다. 원래는 더 빨리 포스팅할 수 있었으나 정확 최대 우도 추정 알고리즘을 구현하는 데 있어서 추정 값이 잘 수렴하지 않아서 석이 나갔고 이로 인해 포스팅이 좀 늦어졌다. 시계열은 정말 어려운 것 같다. 다음 포스팅에서는 ARMA 모형에 대해서 다루려고 한다.
참고자료
Ruey S. Tsay - Analysis of Financial Time Series
Hamilton, J.D. - Time Series Analysis
댓글