본 포스팅에서는 수식을 포함하고 있습니다.
티스토리 피드에서는 수식이 제대로 표시되지 않을 수 있으니
PC 웹 브라우저 또는 모바일 웹 브라우저에서 보시기 바랍니다.
이번 포스팅에서는 ARMA 모형의 차수를 결정하는데 도움이 되는 Extended Sample Autocorrelation Function(ESACF)을 소개하는 논문을 리뷰하고자 한다.
여기서 다루는 내용은 다음과 같다.
3. Extended Sample Autocorrelation Functions
4. Tentative Model Identification
5. Properties of the Iterated AR Estimates
7. Discussions and Comparisons
1. Introduction
이 논문은 일치성(Consistency)을 갖는 AR 모형의 최소 제곱 추정량과 잠정적인 ARMA 모형의 차수를 결정하는 과정을 제안한다.
먼저 $Z_t$가 $ARMA(p, q)$ 모형을 따른다고 해보자.
$$\Phi(B)Z_t = C + \theta(B)a_t\tag{1.1}$$
여기서 $\Phi(B) = U(B)\phi(B) = 1-\Phi_1B-\cdots-\Phi_pB^p$, $U(B) = 1 - U_1B - \cdots - U_dB^d$, $\phi(B) = 1-\phi_1B-\cdots -\phi_{p-d}B^{p-d}$ 그리고 $\theta(B) = 1 - \theta_1B - \cdots - \theta_qB^q$이다. 그리고 $B$는 Backshift Operator이다. $C$는 상수이고 $a_t~N(0, \sigma_a^2)$인 백색 잡음이다.
자기 회귀 모형 AR의 추정은 (최대 우도 추정 방법도 있지만) 주로 최소 제곱 법을 이용하여 추정하게 된다. 하지만 ARMA 모형을 최소 제곱법으로 추정했을 때 얻어진 추정량은 일치성을 갖지 않는다고 한다. 본 논문에서는 일치성을 갖는 추정량을 얻을 수 있는 Iterated Regression Procedure(IRP)를 소개한다. 만약 $Z_t$가 정상성을 갖는 경우 즉 $U(B)=1$인 경우 IRP로 얻어진 AR 파라미터 추정치는 점근적으로 estimated Generalized Yule-Walker 방정식의 솔루션과 동일하다고 한다. 하지만 $Z_t$가 정상성이 아닌 경우에는 이러한 동등성이 만족되지 않는다고 한다.
ARMA 차수를 결정하는 문제는 보통 Sample Autocorrelation Function(SACF), Akaike Information Criterion(AIC), R-array, S-array, Corner 방법 등이 있다고 한다. SACF는 AR 또는 MA 모형의 차수를 결정하는 데는 도움이 되지만 ARMA의 차수를 즉, AR과 MA가 혼합되어 있는 모형의 차수를 결정하는 데에는 도움이 되지 않는다. AIC를 이용한 방법은 일치성을 갖지 않는다고 한다. AIC는 실제 파라미터들보다 더 많은 파라미터를 갖는 모형을 선택한다고 한다. R-array, S-array 방법은 사용하기에 다소 복잡하며 통계적인 성질이 많이 알려져 있지 않다고 한다. Corner 방법 또한 통계적인 성질이 알려져 있지 않다고 한다.
정상성이 아닌 시계열을 정상성을 갖도록 하는 변환으로 차분을 생각한다. 하지만 차분을 얼마나 해야 하는가에 대한 질문이 남아 있다. 이 논문에서는 정상성을 가지는 시계열과 그렇지 않은 시계열에 대하여 모두 적용할 수 있는 ARMA 모형 식별에 대한 방법론을 제시한다. 특히 이 논문에서 제시하는 ESACF를 이용하는 방법은 차분을 신경 쓸 필요가 없다는 이점이 있다.
2. Iterated Regression
먼저 $n$개의 데이터 $Z_1, \ldots, Z_n$가 다음의 $ARMA(p,q)$ 모형으로부터 생성되었다고 가정하자.
$$Z_t = \sum_{l=1}^p\Phi_lZ_{t-l} - \sum_{j=1}^q\theta_ja_{t-j}+a_t\tag{2.1}$$
우리의 첫 목표는 일치성을 갖는 최소 제곱 AR 파라미터를 얻는 것이다. 먼저 주어진 데이터를 이용하여 $AR(p)$ 모형을 최소 제곱법을 이용하여 적합한다.
$$Z_t = \sum_{l=1}^p\Phi_{l(p)}^{(0)}Z_{t-l} + e_{p,t}^{(0)}, t=p+1, \ldots, n\tag{2.2}$$
여기에서 (0)의 의미는 AR 파라미터 추정 시 초기값은 최소 제곱법으로 추정하겠다는 것이다. 그리고 $p$는 AR 모형의 차수이다. 그리고 $e_{p,t}^{(0)}$은 오차항이다.
만약 $Z_t$가 $AR(p)$ 또는 순수 비정상 과정을 따른다면($\phi(B) = 1$) 위에서 구한 최소 제곱 추정량은 일치성을 만족한다고 한다. 하지만 이외의 모형이라면 일치성을 만족하지 않는다고 한다.
식 (2.2)에서 모수 추정치 $\hat{\Phi}_{l(p)}^{(0)}$이 일치성을 만족하지 않는다고 해보자. 식 (2.2)에서 계산한 잔차 $\hat{e}_{p, t}^{(0)} = Z_t - \sum_{l=1}^p\hat{\Phi}_{l(p)}^{(0)}Z_{t-l}$은 데이터 개수 $n$이 크다 하더라도 백색잡음이 안된다고 한다. 따라서 이전 시점들의 값 $\hat{e}_{p, t-j}^{(0)}, j>0$는 $Z_t$를 설명하는데 필요한 정보를 담고 있다. 따라서 아래와 같은 첫 번째 반복 최소 제곱법을 생각하게 된다.
$$Z_t = \sum_{l=1}^p\Phi_{l(p)}^{(1)}Z_{t-l}+\beta_{1(p)}^{(0)}\hat{e}_{p, t-1}^{(0)}+e_{p, t}^{(1)}, t=p+2, \ldots, n\tag{2.3}$$
여기서 (1)은 첫 번째 반복 단계를 의미하며 $e_{p,t}^{(1)}$은 해당 오차항이다. 이때 식 (2.3)에서 추정한 $\hat{\Phi}_{l(p)}^{(1)}$은 $q\leq 1$ 또는 $\phi(B)=1$인 경우에 일치성을 갖는다고 한다. 만약 일치성을 갖지 않는 경우에는 다음 두 번째 반복 단계로 들어간다. 이를 위해 첫 번째 반복 단계로부터 얻어진 잔차를
$$\hat{e}_{p,t}^{(1)} = Z_t - \sum_{l=1}^p\hat{\Phi}_{l(p)}^{(2)}Z_{t-l}-\hat{\beta}_{1(p)}^{(1)}\hat{e}_{p,t-1}^{(0)}$$
라 하면 이 잔차도 $Z_t$를 설명하는데 필요한 정보가 남아 있어서 더 개선할 수 있다. 두 번째 단계에서는 다음과 같은 모형을 최소 제곱법으로 적합한다.
$$Z_t = \sum_{l=1}^p\Phi_{l(p)}^{(2)}Z_{t-l}+\beta_{1(p)}^{(2)}\hat{e}_{p,t-1}^{(0)} + \beta_{2(p)}^{(2)}\hat{e}_{p,t-2}^{(0)}+e_{p,t}^{(2)}, t=p+3, \ldots, n \tag{2.4}$$
이때 $q\leq 2$ 또는 $\phi(B) =1$이라면 식 (2.4)에서 얻어진 AR 파라미터 추정량은 일치성을 갖는다고 한다.
일반적으로 모형을 식별하는 단계에서 ARMA의 실제 차수 $p, q$는 알 수 없다. 이 경우 여러 차수 조합에 대해서 모형을 적합해본다. 이때 $j$ 번째 반복 단계에서 $AR(k)$ 회귀 모형을 다음과 같이 정의하자.
$$Z_t = \sum_{l=1}^k\Phi_{l(k)}^{(j)}Z_{t-l}+\sum_{i=1}^j\beta_{i(k)}^{(j)}\hat{e}_{k, t-i}^{(j-i)}+e_{k, t}^{(j)}, \\ t=p+k+1, \ldots, n, j= 0, \ldots, k= 1,2, \ldots, \tag{2.5}$$
여기서
$$\hat{e}_{k, t}^{(i)} = Z_t-\sum_{l=1}^k\hat{\Phi}_{l(k)}^{(i)}Z_{t-l}-\sum_{h=1}^i\hat{\beta}_{h(k)}^{(i)}\hat{e}_{k, t-h}^{(i-h)}\tag{2.6}$$
는 $i$ 번째 반복 단계에서 최소 제곱법을 이용하여 구해진 잔차이며 $\hat{\Phi}_{l(k)}^{(i)}$, $\hat{\beta}_{h(k)}^{(i)}$는 최소 제곱법을 이용하여 구한 회귀계수들이다.
AR 파라미터의 추정치들은 다음의 점화식을 만족한다고 한다.
$$\hat{\Phi}_{l(k)}^{(j)} = \hat{\Phi}_{l(k+1)}^{(j-1)}-\hat{\Phi}_{l-1(k)}^{(j-1)}\hat{\Phi}_{k+1(k+1)}^{(j-1)}/\hat{\Phi}_{k(k)}^{(j-1)}\tag{2.7}$$
여기서 $\hat{\Phi}_{0(k)}^{(j-1)}=-1, l=1, \ldots, k, k\geq 1, j \geq 1$이다.
정리 2.1은 패스하고 정리 2.2는 특정 조건하에서 반복 단계로 얻어진 AR 파라미터들은 일치성을 갖는다는 것을 말해준다.
3. Extended Sample Autocorrelation Functions
지금까지는 일치성을 갖는 추정량을 구하는 방법에 대해서 알아보았다면 이제는 잠정적으로 ARMA 모형의 차수를 결정하는 방법에 대해서 알아보자. 정리 2.2는 특정 조건하에서 반복 단계로 얻어진 AR 파라미터들은 일치성을 갖는다는 것을 말해준다. 이러한 일치성을 갖는 추정치를 이용하여 우리는 Extended Sample Autocorrelation Functions(ESACF) $r_{j(k)}$을 다음과 같은 방식으로 정의한다.
1 단계) $p=0$이면 $Z_t$는 $MA(q)$를 따르고 SACF는 차수 $q$에서 cut off가 점근적으로 일어난다. 즉,
$$r_{j(0)} \approx 0 if j>q, p=0\tag{3.1}$$
다시 말하면 ESACF $r_{j(0)}$는 우리가 잘 알고 있는 SACF가 된다.
2 단계) $p=1$이면 $Z_t$는 $ARMA(1, q)$를 따른다. 이때 이전 섹션에서 알아본 방법으로 $\Phi_1$의 일치 추정량 $\hat{\Phi}_1$을 얻을 수 있고 $W_t = Z_t-\hat{\Phi}_1$는 점근적으로 $MA(q)$를 따른다고 한다. 하지만 $\hat{\Phi}_1$을 얻으려면 반복 회귀를 적용해야 하는데 반복 횟수는 $q$에 의존한다고 한다. 하지만 $q$를 모르는 것이 문제이다. 따라서 다음의 사항을 고려한다.
만약 $q=0$이면 정리 2.2에 의해 모든 $j\geq 0$에 대하여$\hat{\Phi}_{1(1)}^{(j)}$는 $\Phi_1$으로 확률 수렴한다. 이때 시계열 $W_{1,t}^{(0)}=Z_t-\hat{\Phi}_{1(1)}^{(0)}Z_{t-1}$은 점근적으로 백색 잡음이 된다고 한다. 따라서
$$r_s(W_{1,t}^{(0)}) \approx 0 \text{ for } s > 0$$
이때 $r_s(X_t)$는 시계열 $X_t$의 $s$차 SACF이다.
만약 $q=1$이면 정리 2.2에 의해 모든 $j\geq 1$에 대하여 $\hat{\Phi}_{1(1)}^{(j)}$는 $\Phi_1$으로 확률 수렴한다. 이때 시계열 $W_{1,t}^{(1)}=Z_t-\hat{\Phi}_{1(1)}^{(1)}Z_{t-1}$은 점근적으로 $MA(1)$을 따른다. 따라서
$$r_s(W_{1,t}^{(1)}) \approx 0 \text{ for } s > 1$$
이다.
일반적으로 $j\geq q$에 대하여 $\hat{\Phi}_{1(1)}^{(j)}$은 $\Phi_1$로 확률수렴하고 $W_{1,t}^{(q)}=Z_t-\hat{\Phi}_{1(1)}^{(q)}Z_{t-1}$은 점근적으로 $MA(q)$을 따른다. 따라서
$$r_s(W_{1,t}^{(q)}) \approx 0 \text{ for } s > q$$
이다.
이제 $Z_t$의 1st ESACF를 다음과 같이 정의할 수 있다.
$$r_{j(1)} = r_j(W_{1,t}^{(j)})\tag{3.2}$$
여기서 subscript (1)은 1st ESACF를 의미한다.
$r_{j(1)}$는 다음과 같은 성질이 있다.
$$\begin{align} r_{j(1)} & \approx 0 & \text{ for } j > q, p=1 \\ &\neq 0 & \text{ for } j=q, p=1 \end{align}\tag{3.3}$$
일반적으로 임의의 양의 정수 $k$에 대하여 우리는 $k$th ESACF를 정의할 수 있다.
$$r_{j(k)} = r_j(W_{k,t}^{(j)})\tag{3.4}$$
여기서 $W_{k,t}^{(j)} = Z_t-\sum_{l=1}^k \hat{\Phi}_{l(k)}^{(j)}Z_{t-l}$ 이다. 만약 $Z_t$가 $ARMA(p, q)$를 따른다면 $W_{p,t}^{(j)}$은 $j\geq q$에 대하여 점근적으로 $MA(q)$를 따른다.
$$\begin{align} r_{j(k)} & \approx 0 & \text{ for } j > q, k=p \\ &\neq 0 & \text{ for } j=q, k=p \end{align}\tag{3.5}$$
(3.1)과 (3.5)를 비교해보면 $p$th ESACF도 cutting-off 성질이 있는 것을 알 수 있다. 그리고 ESACF도 결국 변환된 시계열의 SACF이므로 ESACF의 샘플링 성질들도 SACF의 성질들을 이용하여 알아낼 수 있다. 하지만 $k=0$인 경우를 제외하고 ESACF는 $Z_t$의 어떠한 변환의 SACF가 아니라고 한다. 이러한 특징의 효과는 뒤 챕터에서 논의한다고 한다.
- Overfitting -
$ARMA(p,q)$ 모형에서 $p$를 알고 있다고 하자. 그렇다면 (3.5)의 cutting-off 성질은 MA 차수 q를 결정하는데 아주 도움이 된다. 이러한 성질은 정리 2.2에서 소개한 $\hat{\Phi}_{l(k)}^(j)$의 일치성 성질에 의존한다.
실제 문제에서는 $p, q$는 알지 못한다. 특히 AR 차수 $k$가 실제 AR 차수 $p$보다 큰 경우 그리고 반복 회귀의 반복 횟수 $j$가 실제 MA 차수 $q$보다 큰 경우 문제가 되는데 이러한 문제를 Overfitting 이라고 한다. 섹션 6에서는 AR 차수 $k$가 실제 AR 차수 $p$보다 큰 경우의 Overfitting 효과가 MA의 차수를 실제보다 더 증가시키는 것을 증명한다고 한다.
ESACF에 대하여 저자는 다음을 증명한다고 한다. $ARMA(p,q)$ 모형에 대하여
$$\begin{align} r_{j(k)} & \approx c(k-p, j-q) & \text{ for } 0 \leq j-q \leq k-p \\ & = 0 & \text{ for } j-q>k-p\leq 0 \end{align}\tag{3.6}$$
여기서 $c(k-p, j-q)$는 0이 아닌 -1과 1사이의 값을 갖는 상수 혹은 확률 변수이다.
4. Tentative Model Identification
본 섹션에서는 ESACF를 이용하여 ARMA의 차수를 결정하는 방법을 제공한다. 이를 위하여 앞에서 다룬 $r_{j(k)}$을 2-way 테이블로 나타낸다.
또한 실제 모형이 $ARMA(1,2)$인 경우, 이론적인 테이블 형태는 다음과 같다.
이때 * 는 -1과 1 사이의 값을 갖는 상수이며 기호 'X'는 ESACF가 통계적으로 0이 아님을 의미하고 기호 '0'는 통계적으로 0임을 의미한다. 식 (3.6)을 이용하면 실제 모형이 $ARMA(1,2)$ 인 경우 $r_{j(1)} \approx 0, j\geq 3$, $r_{j(k)} \approx 0, j=k+2$이며 이는 위 테이블에서 나타난 결과와 같다.
저자는 ESACF의 테이블을 이용하여 0값이 삼각 경계를 이루는 경우 좌측 상단 꼭지점이 0이 되는 지점을 ARMA의 차수로 잠정적으로 정할 수 있다고 한다. 즉, $k=c_1\geq 1, j-k=c_2\geq 0$이라 한다면 $ARMA(c_1,c_2)$를 잠정 모형을 정한다. 실제 데이터 분석 시에는 유한 샘플을 이용하여 분석하기 때문에 $r_{j(k)}$는 정확히 0은 아닐 것이다. 이때에는 $r_{j(k)}$의 점근 분산을 이용한다. 즉, X는 $r_{j(k)}$이 구간 $\pm 2\sqrt{ \text{asymp.}Var(r_{j(k)})}$ 밖에 있다는 뜻이고 0은 그 사이에 있는 경우를 나타낸다. 점근 분산 $\text{asymp.}Var(r_{j(k)})$은 Bartlett 공식을 사용하여 구하거나 식 (3.4)에서 $W_{k,t}^{(j)}$이 백색 잡음이라는 가정 하에 간단하게 $(n-k-j)^{-1}$로 계산할 수도 있다.
아래의 Example은 저자의 방법과 다른 방법들을 비교하여 저자의 방법이 우월하다는 것을 설명한다.
Example 1.
Box, Jenkins의 온도 데이터를 이용하여 ESACF 테이블을 완성하고 ARMA 모형의 차수를 결정하는 예제이다. 그 결과 $AR(2), q=0$인 것으로 판단된다. 이렇게 차수를 결정한 뒤 모형의 계수를 구한다. 저자는 다른 방법에 비하여 ESACF를 이용하는 방법이 차분을 이용하지 않고도 ARMA 차수를 결정할 수 있다고 한다.
Example 2.
다음으로 저자는 정상성이 아닌 $ARMA(4,1)$ 모형을 다음과 같이 설정하여 시뮬레이션 스터디를 수행한다.
$$(1-B)^2(1-\sqrt{2}B+B^2)Z_t = (1-0.5B)a_t, \sigma^2_a = 1\tag{4.1}$$
이에 대한 ESACF 테이블을 채운 것이 아래 그림이다.
예상대로 5번째 행과 두 번째 열을 좌측 상단 꼭짓점으로 해서 삼각형 형태의 0에 가까운 값을 보이고 있다. 저자는 반복 회귀를 4번 적용하여 각 단계별로 AR 계수를 계산하였다. 이에 대한 결과는 아래 테이블로 나타내었다.
각 단계별로 각 계수들이 실제 계수 값에 가까워지는 즉, 일치성을 보이는 것을 알 수 있다. 저자는 다른 방법은 ARMA 차수를 결정하기 위해 여러 스텝이 존재하지만 저자의 방법은 비교적 직접적으로 차수를 정할 수 있다고 한다.
Example 3.
저자는 Box, Jenkins의 독서 집중력에 대한 데이터를 세 번째 예제로 소개한다. Box, Jenkins는 아래의 $ARMA(p, q)$ 모형을 적합했다고 한다.
$$(1-0.92B)Z_t = 1.45+(1-0.58B)a_t, \hat{\sigma}_a^2 = 0.097\tag{4.2}$$
이에 대한 ESACF 테이블은 아래와 같다.
Indicator Symbol을 보면 두 번째 행과 두 번째 칼럼을 좌측 상단 꼭지점으로 하는 삼각형 패턴을 볼 수 있는데 이때 7번째 칼럼에서 3개의 값이 X로 나타나서 약간 애매하다. 하지만 이 값들은 X의 상한선인 0.15보다 약간 큰 값이라서 이 정도는 X가 아닌 0이라고 봐도 무방하다고 한다. 즉, 저자는 약간의 차이가 있지만 이 예제에서도 ESACF 테이블이 $ARMA(p,q)$를 식별하는데 유효하다는 것을 보여주었다.
Example 4.
여기에서는 카페인 섭취에 대한 데이터를 이용하였다. 저자는 계절성을 나타내는 차수도 ESACF 테이블을 이용하여 식별할 수 있다고 한다.
5. Properties of the Iterated AR Estimates
5.1 A Recursion
여기에서는 식 (2.7)의 점화식이 성립한다는 여러 Lemma를 이용하여 증명한다.
5.2 Relationship to Solutions of Moment Equations
여기에서는 저자가 소개한 반복 회귀 방법을 이용한 AR 계수 추정치와 Yule-Walker 방정식을 이용하여 구한 AR 계수의 추정치의 관계가 점근적으로 동일함을 증명하였다.
5.3 Consistency Properties of Iterated AR Estimates
여기에서는 정상(Stationary) ARMA와 비정상(non Stationary) ARMA인 경우 반복 회귀 방법을 이용한 AR 추정치의 일치성 성질을 소개한다.
6. Properties of the ESACF
여기에서는 ESACF의 점근 성질인 (3.5)와 (3.6)에 대한 증명을 소개한다.
7. Discussions and Comparisons
여기에서는 ESACF의 추가적인 성질과 기존 ARMA 차수를 식별하는 방법들과의 비교 결과를 소개한다.
ESACF는 ARMA 모형 차수의 최대값을 결정한다고 한다. 그 이유는 ESACF가 $Z_t$의 어떤 특정 변환의 SACF가 아니라고 한다(이 부분 잘 모르겠다).
저자는 ESACF 테이블을 이용하여 잠정적으로 ARMA 차수를 결정하고 그다음 반복 회귀를 이용하여 계산된 AR 다항식 $\hat{\Phi}_p^{(q)}(B)$와 $W_{p,t}^{(q)}=\hat{\Phi}_p^{(q)}(B)Z_t$의 SACF의 패턴을 실제 모형을 추정하기 전에 좀 더 조사해보라고 한다.
식 (1.1)에서 $a_t$를 정규분포라 가정했는데 이러한 정규성 가정은 AR 계수 추정치의 일치성에는 영향을 미치지 않는다고 한다.
8. Implementation with Python
8.1 Implementation
여기서는 ESACF 테이블을 만드는 함수를 파이썬으로 구현한다. 원래는 논문에서 제시한 알고리즘을 이해하고 이를 바탕으로 구현하려고 했으나 내가 잘 이해를 못한 것인지 논문에 나온 결과와 일치하지 않았다.
R에서 TSA라는 패키지에 eacf 함수가 ESACF 테이블을 만들어주는데 나도 이 함수를 참고했다. 먼저 필요한 모듈을 임포트 해준다.
import pandas as pd
import numpy as np
import statsmodels.api as sm
import warnings
warnings.filterwarnings('ignore')
from statsmodels.tsa.stattools import acf
아래는 ESACF 테이블을 만들어주는 함수이다. ESACF 테이블을 만드는 과정을 간략히 설명하면 다음과 같다. 먼저 ESACF 테이블의 칼럼은 MA 차수를 의미한다. 또한 이 칼럼은 반복 회귀의 반복 횟수를 결정하게 된다. 1) 반복 회귀에서 얻어진 AR 계수는 식 (2.7)을 이용하여 계산하고 2) 이를 이용하여 식 (3.4)의 $W_{k,t}^{(j)}$를 계산한다. 3) 마지막으로 $W_{k,t}^{(j)}$의 0부터 정해진 AR 차수까지 SACF를 계산하면 하나의 칼럼이 완성되며 이를 MA 차수까지 반복한다.
아래 코드는 ESACF 테이블을 만드는 esacf 함수를 나타낸 것이다. 인자는 시계열 데이터와 테이블에 표시할 최대 AR, MA 차수 그리고 기호를 사용할지 말지를 결정하는 인자가 있다. 기호를 사용하지 않으면 ESACF의 값을 계산하고 기호를 사용할 경우 ESACF값이 0이라고 판단되면 'O' 아니면 'X'로 한다.
def esacf(data, ar_max=7, ma_max = 13, symbol=True):
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]) > 2/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 함수 내부에는 테이블을 만드는 데 필요한 내부 함수가 정의되어있다. 이에 대한 간략한 설명은 다음과 같다.
line 5~8
배열을 한 시점 이전으로 이동시킨다. 이때 배열의 길이는 입력 배열과 똑같이 유지하는데 비는 원소는 Nan으로 처리한다.
line 10~19
OLS를 이용하여 자기 회귀 모형을 적합한 뒤 추정된 계수를 출력한다.
line 21~35
최소 제곱 회귀를 반복적으로 적용하여 AR 계수의 추정치를 행렬로 표시한다. 이때 재귀식 (2.7)을 이용한다.
line 45~51
ESACF 테이블의 각 칼럼을 출력한다.
8.2 Test
이제 구현한 것이 잘 작동하는지 확인해보자. 여기서 사용할 데이터는 Example 1에서 소개한 온도 데이터이다. 해당 데이터는 아래에 첨부해놓았다.
먼저 데이터를 불러온 뒤 숫자형으로 바꿔준다.
data = []
with open('./temperature.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 += line
data = [float(x) for x in data]
이제 ESACF 테이블을 만들어보자. 난 ESACF 테이블을 ESACF 값과 기호로 둘 다 나타내었다.
esacf(data, ar_max=5, ma_max = 8, symbol=False) # ESACF 값 출력
esacf(data, ar_max=5, ma_max = 8, symbol=True)
결과를 보니 해당 데이터는 $AR(2)$ 모형이라고 판단되며 이는 논문에 나온 결과와 일치하였다. 아래 그림은 위 결과와 논문에 나온 결과를 비교한 것이다.
이번 포스팅에서는 Extended Sample Autocorrelation Function에 대해서 알아보았다. 예전부터 ARMA의 차수는 어떻게 결정할 수 있을지에 대하여 궁금했는데 이번 기회에 ESACF를 이용하여 차수를 정할 수 있다는 것을 알게 되어 뿌듯하다. 하지만 알고리즘에 대해선 아직 완벽하게 이해한 것은 아니어서 나중에 시간 되면 한번 더 봐야겠다.
댓글