1 박스-젠킨스 방법론 ¹

자기회귀이동평균(ARMA) 모형을 이해하고, 자기회귀이동평균 모형을 식별하는 도구 자기상관함수(ACF), 편자기상관함수(PACF)를 살펴본다. 그리고 박스-젠킨스 국제항공 여객데이터(1949 - 1960) 천명 단위로 기록된 시계열 데이터를 모형 식별 및 예측을 위해서 위한 박스-젠킨스 방법론을 적용한다.

박스-젠킨스 방법론은 지난 과거 시계열 데이터를 가장 잘 접합시키는 방법론을 제안했다. 박스-젠킨스 시계열 모형 적합 방법론은 총 3단계로 구성된다.

모형 식별(Model Identification): 시계열 데이터의 정상성을 확인하고 계절변동이 있는지도 확인한다. 특히 편자기상관함수(PACF, Partial Autocorrelation Function), 자기상관함수(ACF, AutoCorrelation Function)를 사용해서 자기회귀이동평균 모형 \(p, q\) 차수를 결정한다.
모수 추정(Parameter Estimation): 시계열 모형에 대한 추정은 회귀분석의 최소제곱방법과 유사하지만 가우스-뉴튼 아이디어에 기초해서 컴퓨터를 활요한 수치해석방법을 적용한다.
모형검정(Model Diagnostics): 자기회귀이동평균 모형을 적용시키고 남은 잔차의 정상성을 확인하는데 중점을 두는데, 특히 잔차가 서로 독립이고 시간에 따라 평균과 분산이 일정한지 검증한다. 시계열 데이터의 자기상관을 검정하는 융-박스(Ljung-Box) 통계량, 평균과 분산이 일정한지, 자기상관함수와 편자기상관함수를 사용하여 추가적으로 모형에 누락된 것이 없는지 검정한다.

추가적으로 시계열 데이터에 대한 최적 모형을 정보이론 (AIC, BIC)에 따라 적절한 차수를 선정한 자기회귀이동평균이 채택되면 예측(forecast)함수를 통해 예측을 한다.

1.1 정상성(Stationarity)

존즌앤존슨 분기 수익률 데이터가 전형적인 추세를 갖고 분산도 커져 변동성이 증가하는 시계열 데이터의 전형이다. 이를 위해 등분산성을 맞추도록 로그 변환을 하고 나서, 추세를 제거하기 위해 차분을 한다. 로그변환과 차분을 통해 정상성을 확보한 후에 자기회귀 모형을 적합시킨다.

AirPassengers 데이터는 박스-젠킨스 시계열 방법론에 소개된 국제항공 여객데이터(1949 - 1960)는 시계열의 붓꽃(iris) 데이터에 해당되는데 크게 3가지 특성이 있다.

추세(trend)를 갖는다; 1949년 이후 국제항공 이용자수가 지속적으로 증가.
변동성(heteroscedasticity)이 커지고 있다; 연도가 지날수록 국제항공 이용자수의 변동성이 커지고 있다.
계절성(seasonality)을 갖는다; 12개월마다 반복되는 계절성이 뚜렷이 관측된다.

autoplot() 함수를 통해서 정상성을 확보하는 차분과 로그변환에 대해서 시각화해보자.

# 0. 환경설정 -----
library(tidyquant)
library(timetk)
library(sweep)
library(extrafont)
library(forecast)
library(ggfortify)
library(cowplot)
loadfonts()
options(family="NanumGothic")

# 1. 데이터 -----
## 1.1. ts 데이터 가져오기
data(AirPassengers)

## 1.2. 시계열 데이터 정상성
ap_ts <- autoplot(AirPassengers)
ap_diff_ts <- autoplot(diff(AirPassengers))
ap_log_ts <- autoplot(log(AirPassengers))
ap_log_diff_ts <- autoplot(diff(log(AirPassengers)))

plot_grid(ap_ts, ap_diff_ts, ap_log_ts, ap_log_diff_ts, nrow=2,
          labels = c("원 시계열", "차분 시계열", "로그변환 시계열", "로그 차분 시계열"))

2 자기회귀이동평균(ARMA) 모형

자기회귀이동평균 모형을 짚고 넘어가기 전에 선형회귀모형을 먼저 살펴보자.

\(Y_i = \beta X_i + \epsilon_i\)

선형회귀모형은 수학적으로 표현하면 위와 같고 오차항 \(\epsilon_i\)는 \(i.i.d\) 가정을 한다. 즉, 동일한 분산을 갖는 서로 독립적인 정규분포다.

자기회귀(AR) 모형 : \(Y_t = \phi Y_{t-1} + \epsilon_t\), 여기서 \(\epsilon_t\)는 백색잡음
이동평균(MA) 모형 : \(\epsilon_t= W_t + \theta W_{t-1}\), 여기서 \(W_t\)는 백색잡음
자기회귀이동평균(ARMA) 모형 : \(X_t = \phi X_{t-1} + W_t + \theta W_{t-1}\)

즉, 정상성 조건이 충족된다는 가정하에 현재 관측점은 과거 \(p\)차까지 자기회귀를 하고 남은 잔차의 상관관계를 필터링하여 백색잡음을 만드는 모형이다.

자기회귀모형(AR)은 세가지 모수를 갖는다.

평균: \(\mu\)
자기회귀계수: \(\phi\)
백색잡음 분산: \(\sigma_{\epsilon}^{2}\)

이동평균모형(MA)도 세가지 모수를 갖는다.

평균: \(\mu\)
이동평균회귀계수: \(\theta\)
백색잡음 분산: \(\sigma_{\epsilon}^{2}\)

2.1 자기회귀이동평균(ARMA) 모형 식별

로그변환(log transformation)과 차분(difference)을 통해 정상성과 더불어 계절변동성을 제거한 후에 편자기상관함수(PACF, Partial Autocorrelation Function), 자기상관함수(ACF, AutoCorrelation Function)를 도식화하고 이를 통해 자기회귀이동평균 차수를 선정한다.

	자기회귀: AR(\(p\))	이동평균: MA(\(p\))	자기회귀이동평균: ARMA(\(p,q\))
자기상관함수(ACF)	지수적 감소, 축퇴하는 사인 형태(tail off)	\(q+1\) 차항부터 절단모양	\(q+1\) 차항부터 지수적 감소 혹은 축퇴하는 사인형태(tail off)
편자기상관함수(PACF)	\(p+1\) 차항부터 절단모양	지수적 감소, 축퇴하는 사인 형태(tail off)	\(p+1\) 차항부터 지수적 감소 혹은 축퇴하는 사인형태(tail off)

3 국제여객 이용객 시계열 데이터 ²

국제여객 이용객 시계열 데이터가 dataset 팩키지에 들어있는데 자료형이 ts 시계열 자료형으로 되어 있어 데이터프레임으로 변화하여 작업하는 것이 수월하다. 이를 위해서 timetk 팩키지 tk_tbl() 함수를 사용하여 데이터프레임으로 변환하고 나서 ggplot()으로 시각화한다.

## 1.3. ts 데이터프레임 변환
ap_df <- AirPassengers %>% 
  tk_tbl(timetk_idx = FALSE) %>% 
  mutate(date = as.Date(as.yearmon(index))) %>% 
  rename(passenger = value)

ap_df %>% 
  ggplot(aes(x=date, y=passenger)) +
    geom_line() +
    labs(x="", title="국제 항공승객수", subtitle="박스-젠킨스 항공여객데이터(1949 - 1960),  단위:천명", y="") +
    theme_minimal(base_family = "NanumGothic") +
    geom_vline(xintercept=lubridate::ymd("1959-01-01"), linetype=2, color="blue")

3.1 훈련/시험 데이터 분할

전체 데이터가 1949년부터 1960년 약 20년정도 되는데 최근 2년을 시계열 예측모형을 시험하는 기간으로 설정하여 1959년 이후를 시험기간으로 두고, 앞선기간을 시계열 예측모형 훈련기간으로 둔다.

## 1.3. 훈련/시험 데이터 분할 -----
ap_train_df <- ap_df %>% filter(date <= "1958-12-31")
ap_test_df  <- ap_df %>% filter(date >= "1959-01-01")

3.2 ARIMA 모형적합

forecast 팩키지 auto.arima() 함수를 사용해서 예측모형을 생성하고 모형에서 예측한 값을 forecast() 함수로 도출하고 이를 ggfortify 팩키지 fortify() 함수를 사용해서 데이터프레임으로 변환시키는 방법을 많이 사용했으나, 최근에 시계열 broom 버젼인 sweep 팩키지가 출시되어서 sw_sweep() 함수로 데이터프레임 변환이 수월해졌다.

# 2. 모형적합 -----
## 2.1. 데이터프레임을 ts 객체 변환 -----
ap_train_ts <- tk_ts(ap_train_df$passenger, start=1949, freq =12, silent = TRUE)

## 2.2. ARIMA 모형 적합 -----
ap_fit <- auto.arima(ap_train_ts)

## 2.3. ARIMA 모형 예측 -----
ap_fcst <- forecast(ap_fit, h=24, level = c(95, 80))

## 2.4. 데이터프레임 변환 -----
ap_fcst_df <- sw_sweep(ap_fcst, timekit_idx = TRUE, rename_index = "date") %>% 
  mutate(date = as.Date(as.yearmon(date)))

## 2.5. 표로 살펴보기 -----
ap_fcst_df %>% 
  arrange(desc(date)) %>% 
  DT::datatable() %>% 
  DT::formatRound(c(3:7), digits = 0)

3.3 ARIMA 모형적합 시각화

원 시계열 데이터에 추가하여 예측 데이터도 모두 정리된 데이터프레임이 준비되었으니 이를 ggplot으로 시각화하여 원데이터, 점예측값, 신뢰구간도 함께 시각화한다.

## 2.6. 그래프 시각화 -----
ap_fcst_df %>% 
  ggplot(aes(x=date, y=value, color=key)) +
    geom_line(size=1) +
    geom_ribbon(mapping = aes_string(x = 'date', ymin = '`lo.95`', ymax = '`hi.95`'), alpha = 0.5, fill="gray") +
    geom_line(aes(y=hi.95), color="pink") +
    geom_line(aes(y=lo.95), color="pink") +
    labs(x="", title="국제 항공승객수", subtitle="박스-젠킨스 항공여객데이터(1949 - 1960),  단위:천명", y="") +
    theme_minimal(base_family = "NanumGothic") +
    geom_line(data=ap_test_df, aes(x=date + lubridate::dyears(0), y=passenger), color="black") +
    scale_x_date(date_labels = "%Y") +
    geom_vline(xintercept=lubridate::ymd("1958-12-01"), linetype=2, color="blue") +
    scale_color_manual(values=c("black", "blue"))

4 추억의 ARIMA

arima.sim() 함수를 통해 AR, MA, ARMA, ARIMA 확률과정을 따르는 데이터를 생성할 수 있다. 정답을 알고 있는 상황에서 ARIMA모형을 식별할 수 있다.

4.1 자기회귀 AR(1) 모형

자기회귀 평균 \(\mu=7\)을 갖고 상관계수 \(\phi=0.7\)을 갖는 모형은 다음과 같다.

\(Y_t = 7 + 0.7(Y_{t-1} -7) + W_t\)

상기 자기회귀모형을 적합시키면 sarima 함수를 통해 적합시킨 결과가 거의 유사함을 확인된다. 또한, acf2 함수를 통해 자기상관함수는 지수적 감소하고, 편자기상관함수는 2차항부터 절단모양이 관찰된다.

# AR(1) 모의시험 모형------------------------------------------------ 
ar_yt <- arima.sim(list(order = c(1, 0, 0), ar = 0.7), n = 250) + 7

# ACF, PACF
astsa::acf2(ar_yt)

        ACF  PACF
 [1,]  0.71  0.71
 [2,]  0.50 -0.01
 [3,]  0.34 -0.03
 [4,]  0.26  0.06
 [5,]  0.20  0.01
 [6,]  0.16  0.00
 [7,]  0.10 -0.04
 [8,]  0.03 -0.07
 [9,]  0.02  0.07
[10,]  0.01 -0.02
[11,]  0.02  0.02
[12,]  0.02  0.02
[13,]  0.01 -0.03
[14,] -0.03 -0.04
[15,] -0.03  0.02
[16,]  0.01  0.06
[17,]  0.04  0.04
[18,]  0.04 -0.05
[19,]  0.00 -0.05
[20,] -0.01  0.04
[21,] -0.06 -0.12
[22,] -0.06  0.02
[23,] -0.09 -0.06
[24,] -0.11 -0.04
[25,] -0.12  0.02
[26,] -0.08  0.05

ar_yt_fit <- astsa::sarima(ar_yt, p = 1, d = 0, q = 0)

initial  value 0.348793 
iter   2 value -0.006755
iter   3 value -0.006790
iter   4 value -0.006810
iter   5 value -0.006811
iter   6 value -0.006818
iter   7 value -0.006818
iter   8 value -0.006818
iter   9 value -0.006818
iter  10 value -0.006818
iter  11 value -0.006818
iter  11 value -0.006818
final  value -0.006818 
converged
initial  value -0.006970 
iter   2 value -0.006975
iter   3 value -0.006979
iter   4 value -0.006980
iter   5 value -0.006984
iter   6 value -0.006984
iter   7 value -0.006984
iter   8 value -0.006984
iter   8 value -0.006984
iter   8 value -0.006984
final  value -0.006984 
converged

ar_yt_fit$ttable

      Estimate     SE t.value p.value
ar1     0.7148 0.0443 16.1447       0
xmean   6.7624 0.2178 31.0447       0

4.2 이동평균 MA(1) 모형

\(\theta=0.7\) 모수를 갖는 시계열을 sarima() 함수 이동평균 모형에 적합시키게 되면 모수의 유사함이 확인된다.

\(Y_t = W_t + 0.7W_{t-1}\)

또한, acf2 함수를 통해 자기상관함수는 2차항부터 절단모향이 관찰되고, 편자기상관함수는 지수적 감소, 축퇴하는 사인형태를 보인다.

# MA(1) 모의시험 모형------------------------------------------------ 
ma_yt <- arima.sim(list(order = c(0, 0, 1), ma = 0.7), n = 250)

# ACF, PACF
astsa::acf2(ma_yt)

        ACF  PACF
 [1,]  0.45  0.45
 [2,] -0.10 -0.38
 [3,] -0.11  0.17
 [4,] -0.07 -0.18
 [5,] -0.13 -0.05
 [6,] -0.09  0.00
 [7,]  0.02  0.02
 [8,]  0.05 -0.01
 [9,]  0.05  0.06
[10,]  0.05 -0.03
[11,]  0.05  0.07
[12,]  0.04 -0.01
[13,]  0.03  0.05
[14,]  0.00 -0.04
[15,] -0.07 -0.04
[16,] -0.07 -0.01
[17,] -0.05 -0.04
[18,]  0.01  0.06
[19,]  0.01 -0.07
[20,] -0.06 -0.07
[21,] -0.08 -0.03
[22,] -0.02  0.00
[23,]  0.03  0.01
[24,]  0.05  0.04
[25,]  0.10  0.08
[26,]  0.04 -0.08

ma_yt_fit <- astsa::sarima(ma_yt, p = 0, d = 0, q = 1)

initial  value 0.247306 
iter   2 value 0.079966
iter   3 value 0.047366
iter   4 value 0.041614
iter   5 value 0.038363
iter   6 value 0.037646
iter   7 value 0.037619
iter   8 value 0.037604
iter   9 value 0.037604
iter  10 value 0.037604
iter  10 value 0.037604
iter  10 value 0.037604
final  value 0.037604 
converged
initial  value 0.038915 
iter   2 value 0.038912
iter   3 value 0.038910
iter   4 value 0.038910
iter   5 value 0.038910
iter   5 value 0.038910
iter   5 value 0.038910
final  value 0.038910 
converged

ma_yt_fit$ttable

      Estimate     SE t.value p.value
ma1     0.7011 0.0437 16.0405  0.0000
xmean  -0.0101 0.1115 -0.0903  0.9281

4.3 자기회귀이동평균 ARMA(1,1) 모형

\(\theta=0.7\), \(\phi=0.7\) 모수를 갖는 시계열을 sarima() 함수에 ARMA 모형을 적합시키게 되면 모수의 유사함이 확인된다.

\(Y_t = 0.7Y_{t-1} + W_t + 0.7W_{t-1}\)

또한, acf2 함수를 통해 자기상관함수와 편자기상관함수는 모두 지수적 감소, 축퇴하는 사인형태를 보인다.

# ARMA(1,1) 모의시험 모형------------------------------------------------ 
arma_yt <- arima.sim(list(order = c(1, 0, 1), ar=0.7, ma = 0.7), n = 250)

# ACF, PACF
astsa::acf2(arma_yt)

       ACF  PACF
 [1,] 0.84  0.84
 [2,] 0.58 -0.47
 [3,] 0.39  0.31
 [4,] 0.29 -0.08
 [5,] 0.24  0.08
 [6,] 0.21 -0.03
 [7,] 0.19  0.10
 [8,] 0.19 -0.01
 [9,] 0.21  0.15
[10,] 0.24  0.00
[11,] 0.27  0.09
[12,] 0.28  0.01
[13,] 0.25 -0.09
[14,] 0.19  0.03
[15,] 0.16  0.10
[16,] 0.17 -0.02
[17,] 0.17  0.02
[18,] 0.17  0.01
[19,] 0.14 -0.07
[20,] 0.09 -0.05
[21,] 0.05  0.02
[22,] 0.04  0.00
[23,] 0.06  0.07
[24,] 0.10  0.05
[25,] 0.13 -0.02
[26,] 0.14  0.04

arma_yt_fit <- astsa::sarima(arma_yt, p = 1, d = 0, q = 1)

initial  value 0.773040 
iter   2 value 0.010183
iter   3 value 0.010047
iter   4 value -0.024920
iter   5 value -0.041487
iter   6 value -0.042491
iter   7 value -0.042766
iter   8 value -0.042771
iter   9 value -0.042771
iter  10 value -0.042771
iter  11 value -0.042771
iter  12 value -0.042771
iter  13 value -0.042771
iter  14 value -0.042771
iter  14 value -0.042771
iter  14 value -0.042771
final  value -0.042771 
converged
initial  value -0.039603 
iter   2 value -0.039612
iter   3 value -0.039619
iter   4 value -0.039620
iter   5 value -0.039628
iter   6 value -0.039628
iter   7 value -0.039628
iter   7 value -0.039628
iter   7 value -0.039628
final  value -0.039628 
converged

arma_yt_fit$ttable

      Estimate     SE t.value p.value
ar1     0.7082 0.0470 15.0675  0.0000
ma1     0.7183 0.0487 14.7394  0.0000
xmean  -0.4983 0.3524 -1.4142  0.1586

4.4 정보이론 기반 자기회귀이동평균(ARMA) 모형 선정

AIC와 BIC를 활용하여 최적차수를 갖는 자기회귀이동평균 모형을 선택한다. AIC와 BIC는 모수가 추가됨에 따라 벌점을 부과하는데 AIC는 \(k=2\), BIC는 \(k=log(n)\)이 선택된다. 물론, AIC와 BIC는 가장 작은 값을 선택하고 이에 해당되는 모형이 ARMA(\(p,q\)) 최적 모형이 된다.

\(\mathrm{AIC} = 2(p+q) - 2\ln(\hat L)\)
\(\mathrm{BIC} = {\ln(n)(p+q) - 2\ln({\hat L})}\)

\(i.i.d\) 가정을 만족하는 정규과정인 경우 오차제곱(\(\widehat{\sigma_e^2}\))으로 BIC를 재표현할 수 있다.

\(\mathrm{BIC} = n \cdot \ln(\widehat{\sigma_e^2}) + (p+q) \cdot \ln(n)\)

4.5 잔차 분석

sarima() 함수를 통해 잔차분석을 효과적으로 수행된다.

표준화잔차(Standardized residuals)
잔차 자기상관함수(ACF)
정규분포 적합 QQ-플롯
융-박스(Ljung-Box) Q-통계량 p-값

4.6 예측(forecast)

시계열 모형에 갖는 이유는 시간에 따른 미래 시점에 대한 예측에 관심을 두기 때문이다.

데이터 과학 – 금융(Finance)

자기회귀이동평균(ARMA) 모형 - airpassenger

1 박스-젠킨스 방법론 ¹

1.1 정상성(Stationarity)

2 자기회귀이동평균(ARMA) 모형

2.1 자기회귀이동평균(ARMA) 모형 식별

3 국제여객 이용객 시계열 데이터 ²

3.1 훈련/시험 데이터 분할

3.2 ARIMA 모형적합

3.3 ARIMA 모형적합 시각화

4 추억의 ARIMA

4.1 자기회귀 AR(1) 모형

4.2 이동평균 MA(1) 모형

4.3 자기회귀이동평균 ARMA(1,1) 모형

4.4 정보이론 기반 자기회귀이동평균(ARMA) 모형 선정

4.5 잔차 분석

4.6 예측(forecast)

데이터 과학 – 금융(Finance)

자기회귀이동평균(ARMA) 모형 - airpassenger

1 박스-젠킨스 방법론 1

1.1 정상성(Stationarity)

2 자기회귀이동평균(ARMA) 모형

2.1 자기회귀이동평균(ARMA) 모형 식별

3 국제여객 이용객 시계열 데이터 2

3.1 훈련/시험 데이터 분할

3.2 ARIMA 모형적합

3.3 ARIMA 모형적합 시각화

4 추억의 ARIMA

4.1 자기회귀 AR(1) 모형

4.2 이동평균 MA(1) 모형

4.3 자기회귀이동평균 ARMA(1,1) 모형

4.4 정보이론 기반 자기회귀이동평균(ARMA) 모형 선정

4.5 잔차 분석

4.6 예측(forecast)

1 박스-젠킨스 방법론 ¹

3 국제여객 이용객 시계열 데이터 ²