시계열 분석 시리즈 IV (Time Series Toolbox)

 

시계열 예측 도구

과거 데이터로 올해 매출량을 예측한다고 해보자. 몇 가지 원초적인 방법들이 있을 수 있다. 단순하게 지금까지 평균 매출량 정도 팔릴 것으로 예측하는 방법, 아니면 작년의 매출량 보다는 조금 더 팔릴 것으로 예측하는 방법 등이 있을 것이다. 이런 방법은 어려운 시계열 데이터 분석을 배우지 않아도 실제 사업하는 분들이 현장에서 하고 있는 방법이다. 과연 데이터 분석가인 우리가 하는 정교한 시계열 데이터 예측 방법이 감으로 하는 상인들의 예측보다 더 정확할까? 정확하다면 얼마나 더 정확할까?
이번 시리즈에서는 시계열 예측을 수행하는 절차와 예측하는 데 사용하는 도구, 예측 결과를 확인하는 방법 등에 대해 다룬다.

Tidy한 예측 절차

시계열 예측은 데이터를 tidy 한 포맷으로 정리하고, 결측치 제거, log 변환 등 예측에 적합하도록 데이터를 준비하는 것으로 시작한다. 적용하는 모델에 따라 데이터 요구사항이 다르기 때문에 모델에 맞도록 데이터를 준비해야 한다.

# Data 준비
gdppc <- global_economy %>% 
  mutate(GDP_per_capita = GDP / Population)   # 1인당 GDP 계산 
#  A tsibble: 15,150 x 10 [1Y] , Key: Country [263]

데이터 시각화 단계는 준비된 데이터를 가시화시켜 패턴과 특성을 파악하는 과정이다. 시각화 단계는 본 시리즈 II에서 소개된 각종 플롯을 활용하여 데이터를 파악한다. 아래는 국가별 경제지표에 대한 데이터에서 한국 데이터만 발췌하여 가시화한 소스이다.

# Data Visualization 
gdppc %>% 
  filter(Country == "Korea, Rep.") |>
  autoplot(GDP_per_capita) +
  labs(y = "$US", title = "GDP per capita for Korea")

예측 모델에는 여러 가지가 있으며 fable 패키지 포함되어 있다. 우리는 데이터의 특성에 따라 적합한 모델을 지정하여 학습 시키면 된다. 예를들어 시계열 회귀를 모델로 지정하려면 TSLM(Time Series Linear Model)을 사용해 모델을 정의한다. “TSLM(GDP_per_capita ~ trend())”은 GDP_per_capita를 trend()로 선형회귀 하라고 지정한 것이다.

TSLM(GDP_per_capita ~ trend())  # GDP_per_capita를 trend()로 시계열 선형회귀 

<TSLM model definition>

통상은 모델 지정(Specify)과 학습(적합) 과정을 model() 함수로 통합하여 한꺼번에 진행한다. 아래 예제는 위에서 준비한 gdppc 데이터를 TSLM 모델로 지정하여 model()함수로 학습 시킨것이다.

# Train the model (estimate)  : 적합(학습)시키기 
fit <- gdppc %>%                 # 학습된 모델명은 fit로 저장  
  model(trend_model = TSLM(GDP_per_capita ~ trend()))  # model로 TSLM 지정 

fit  %>% print(n=5)  # A mable: 263 x 2,    Key:   Country [263]

# A mable: 263 x 2
# Key:     Country [263]
  Country        trend_model
  <fct>              <model>
1 Afghanistan         <TSLM>
2 Albania             <TSLM>
3 Algeria             <TSLM>
4 American Samoa      <TSLM>
5 Andorra             <TSLM>
# ℹ 258 more rows

학습 결과인 fit를 보면 Key 변수인 국가(Country)별 데이터를 적용한 회귀모델 263개가 생성되었음을 알 수 있다. 이 학습결과를 model table이란 의미에서 mable 객체라고 한다.

학습이 끝나면 학습된 모델이 얼마나 잘 학습되었는지를 평가지표에 따라 평가(Evaluate)한다. 평가 지표는 예측값과 실제값이 얼마나 차이가 나는지를 확인하는 지표이다. 회귀모델에 대한 다양한 평가지표에 대해서는 본 글의 “예측정확도 평가” 챕터에서 자세하게 다룬다.

평가 결과를 수용할 만한 수준이 되면, 학습시킨 모델로 예측을 수행한다. 예측은 학습된 모델(mable)을 forecast()함수에 파이프라이닝 시키면서 예측할 기간을 지정해서 예측한다. 예측된 결과 모델을 forecast model이란 의미에서 fable 객체라고 한다.
아래 소스는 위에서 적합 시킨 mable객체 fit를 forecast() 함수에 파이프라이닝 시겨 3년간 예측하고 이중 “Korea, Rep.”만 발췌하여 가시화 시킨 것이다.

# Produce forecasts : 예측수행하기  

fit %>% 
  forecast(h = 3 ) %>%  #  3년간 예측해서 fable로 리턴 
# A fable: 789 x 5 [1Y],  Key:     Country, .model [263] --> 263개, 각 3개년 
  filter(Country == "Korea, Rep.") %>%  
  autoplot(gdppc) +
  labs(title = "한국의 개인당 GDP ", y = "$US")

몇 가지 단순한 예측 기법

복잡한 모델을 통해 정밀하게 예측하는 경우도 있지만, 때론 아주 단순하게 예측하는 방법도 있다. 단순 예측기법에는 지금까지의 평균값으로 예측하는 방법, 마지막 값이 계속될 것으로 나이브하게 생각한다고 하여 이름 붙여진 나이브 기법, 나이브하게 예측하되 이전 계절의 값으로 예측하는 계절성 나이브(SNAIVE) 등의 방법이 있다. 이런 방법들은 실제 예측에 사용되기보다는 예측 모델의 참조모델로 예측 결과를 비교하기 위한 벤치마킹 모델로 사용되곤 한다.

아래는 aus_production 데이터를 필터링하여 1992년 1분기부터 2006년 4분기까지 60개 분기의 데이터를 학습용 데이터로 만들고, 이후 14개 분기의 데이터로 호주의 맥주 생산량을 예측한 예제이다. 예측 방법은 단순한 기법인 평균, 나이브, 계절 나이브의 세가지 방법으로 예측하였다.

# Set training data from 1992 to 2006
train <- aus_production %>%
  filter_index("1992 Q1" ~ "2006 Q4")

# Fit the models   ---> beer_fit 모델은 mable 
beer_fit <- train %>%          # train 데이터를 세가지로 학습해서 beer_fit 로 
  model(                       # 세가지 모델(MEAN, NAIVE, SNAIVE)로 학습
    평균 = MEAN(Beer),         # 여러 변수중 맥주 생산량을 평균으로 예측 
    나이브 = NAIVE(Beer),
    계절나이브 = SNAIVE(Beer)
  )    # beer_fit 모델에는 "평균", "나이브", "계절나이브" 세가지 모델이 들어있음 

# Generate forecas ts for 14 quarters
beer_fc <- beer_fit %>% forecast(h = 14)    # 학습된 모델(beer_fit)로 14개분기 예측 

# Plot forecasts against actual values
beer_fc %>%
  autoplot(train, level = NULL) +   
  autolayer(.vars = Beer, #  2007년 Q1 이후 데이터 출력 
    filter_index(aus_production, "2007 Q1" ~ .),    
    colour = "black"
  ) +
  labs(
    y = "Megalitres",
    title = "분기별 맥주 생산량 예측 "
  ) +
  guides(colour = guide_legend(title = "Forecast"))    

예측 결과를 보면, 평균, 나이브, 계절나이브로 예측한 것이 원래 데이터와 오버랩되어 도시되었다.

적합값과 잔차

참고 : 적합값이라고 사용하는 이유는 적합값은 학습된(fitted) 모델이 학습한 값이란 의미이고, 예측값은 학습된 모델이 forcast 함수를 이용하여 예측한 값을 의미하기 때문에 둘을 구분하기 위해서 적합값이란 용어를 사용하였다.

적합값 \(\hat{y}_{t|t-1}\) 은 \(y_1, y_2.., y_{t-1}\) 값으로 \(y_{t|t-1}\)를 예측한 값이며 이를 적합값(fitted value)라고 한다.

평균값으로 예측한 경우에는 \(\hat{y}_t = \hat{c}\) 로 표현할 수 있다.

잔차(Resudial)는 관측값과 적합값의 차이를 \(e = y_t - \hat{y}_t\) 의미한다. 모델의 적합값과 잔차는 augment() 함수를 써서 확인 할 수 있다.

augment(beer_fit) 

# A tsibble: 180 x 6 [1Q]
# Key:       .model [3]
   .model Quarter  Beer .fitted .resid .innov
   <chr>    <qtr> <dbl>   <dbl>  <dbl>  <dbl>
 1 평균   1992 Q1   443    436.   6.55   6.55
 2 평균   1992 Q2   410    436. -26.4  -26.4 
 3 평균   1992 Q3   420    436. -16.4  -16.4 
 4 평균   1992 Q4   532    436.  95.6   95.6 
 5 평균   1993 Q1   433    436.  -3.45  -3.45
 6 평균   1993 Q2   421    436. -15.4  -15.4 
 7 평균   1993 Q3   410    436. -26.4  -26.4 
 8 평균   1993 Q4   512    436.  75.6   75.6 
 9 평균   1994 Q1   449    436.  12.6   12.6 
10 평균   1994 Q2   381    436. -55.4  -55.4 
# ℹ 170 more rows

augment() 함수는 적합값(.fitted), 잔차(.resid) 및 혁신잔차(.innov)를 제시한다.
참고 혁신잔차(innovation residuals) : 혁신잔차는 필요에 따라서 데이터를 변환(예 로그 변환)했을 때 변환된 스케일의 잔차를 의미한다.

교재에서 제시한 다음 예제는 시계열 데이터 분석 및 예측에 대한 프로세스를 잘 정리해 놓았다. 다음 문제를 읽고 문제를 어떻게 해결할지를 생각하면서 예제를 따라가 보자.

문제 : gafa_stock 데이터는 구글, 애플, 페이스북, 아마존에 대한 주식가격이 들어있는 시계열 데이터이다. 구글사의 종가 기준 주식 시세를 예측해 보고자 한다. 다음 순서대로 문제를 해결하시오.

1. gafa_stock 데이터의 구조를 파악하시오.
2. 2015년 이후 구글사 데이터만을 발췌하시오.
3. 발췌한 데이터의 인덱스가 불규칙한 경우 이를 개선하시오.
4. 데이터를 분할하되, 2015년 데이터를 train_set으로 하고, 가시화하시오.
5. 2016년 1월 데이터를 test_set으로 하시오.
6. train_set으로 다음 세 가지 모델을 적합(fit)/학습시킨 모델을 생성하시오.

• 평균으로 예측, 단순하게(NAIVE) 에측, drift 예측
• drift 예측은 (dirft 모델 : https://otexts.com/fpp3/simple-methods.html 랜덤워크 모델 참조 : https://otexts.com/fpp3/stationarity.html )

7. test_set으로 예측을 수행하시오.
8. 예측 결과를 테이블 형태와 신뢰구간을 포함하여 가시화하시오.

#  1. gafa_stock 데이터의 구조를 파악하시오. 
#  2. 2015년 이후 구글사 데이터만을 발췌하시오. 
#  3. 발췌한 데이터의 인덱스가 불규칙한 경우 이를 개선하시오.

google_stock <- gafa_stock %>%
  filter(Symbol == "GOOG" & year(Date) >= 2015) %>%   
  mutate(Day = row_number()) %>%   # 규칙성 개선 
  update_tsibble(index = Day, regular = TRUE) 
 
#  4. 2015년 데이터를 google_2015(train_set)으로 하고, 가시화 시키시오. 
google_2015 <- google_stock %>%
  filter(year(Date) == 2015) 
google_2015 %>% autoplot(Close)

# 5. 2016년 1월 데이터를  test_set으로 하시오.
test_set <- google_stock %>%
  filter(yearmonth(Date) == yearmonth("2016. 1"))  # 이렇게 비교해도 되고.

test_set <- google_stock %>%
  filter(as.character(yearmonth(Date)) == "2016 1")  # 요렇게 해도 되겠지...

# 6. google_2015으로 다음 세가지 모델을 적합(fit)/학습 시킨 모델을 생성하시오.

fitted_model <- model(google_2015, mean = MEAN(Close),
                                 naive = NAIVE(Close),
                                 drift = NAIVE(Close ~ drift()))  
                 # drift 모델 --> 교재 chapter 5.2 혹은 9.1 참조 
# 또는 tidyverse 하게..
fitted_model <- google_2015 %>%
   model(mean = MEAN(Close),
         naive = NAIVE(Close),
         drift = NAIVE(Close ~ drift()))

#  7. test_set으로 예측을 수행하시오.
forecast_result<- fitted_model %>%
  forecast(new_data = test_set) 
 
# 8. 예측결과 출력 및 가사화  
forecast_result  # fable

# A fable: 57 x 11 [1]
# Key:     Symbol, .model [3]
   Symbol .model   Day        Close .mean Date        Open  High   Low Adj_Close
   <chr>  <chr>  <int>       <dist> <dbl> <date>     <dbl> <dbl> <dbl>     <dbl>
 1 GOOG   mean     253 N(602, 6766)  602. 2016-01-04  743   744.  731.      742.
 2 GOOG   mean     254 N(602, 6766)  602. 2016-01-05  746.  752   739.      743.
 3 GOOG   mean     255 N(602, 6766)  602. 2016-01-06  730   747.  729.      744.
 4 GOOG   mean     256 N(602, 6766)  602. 2016-01-07  730.  738.  719.      726.
 5 GOOG   mean     257 N(602, 6766)  602. 2016-01-08  731.  733.  713       714.
 6 GOOG   mean     258 N(602, 6766)  602. 2016-01-11  717.  719.  704.      716.
 7 GOOG   mean     259 N(602, 6766)  602. 2016-01-12  722.  729.  717.      726.
 8 GOOG   mean     260 N(602, 6766)  602. 2016-01-13  731.  735.  699.      701.
 9 GOOG   mean     261 N(602, 6766)  602. 2016-01-14  705.  722.  689.      715.
10 GOOG   mean     262 N(602, 6766)  602. 2016-01-15  692.  707.  685.      694.
# ℹ 47 more rows
# ℹ 1 more variable: Volume <dbl>

# 예측결과 가시화 
forecast_result %>%
  autoplot(google_2015) + # 2015년 데이터와 예측값이 도시 
  autolayer(test_set, color ="red")+ # 실제 2016,1 월 데이터 
  labs(title = "구글 주식 예측 ")

Plot variable not specified, automatically selected `.vars = Open`

잔차의 진단

잔차(residual)의 진단은 어떤 모델이 적절하게 예측을 수행했는지를 판단할 때 유용하다. 제대로 예측한 모델의 잔차는 다음과 같은 특징을 갖는다.

• 잔차(residual는 상관관계가 없다. 잔차 사이에 상관관계(correlation)가 있다면, 잔차에 예측값을 계산할 때 사용한 정보가 남아 있다는 의미이다.
• 잔차의 평균은 0이다. 잔차의 평균이 0이 아니라면, 예측값이 편향(bias)되어 있다는 의미이다.
• 잔차의 분산이 상수이다. (등분산이다)
• 잔차는 정규분포를 따른다.

잔차의 진단은 적합시킨 모델의 잔차가 위 가정을 만족하는지를 확인하는 과정이다. 앞에서 예측한 예제인 google_2015의 innovation 잔차로 확인해 보자.

aug <- google_2015 %>%    
  model(NAIVE(Close)) %>%  
  augment()    #  NAIVE(Close) 모델 fitting 결과들의 결과 출력 
  #   Symbol, .model,  Day, Close,  .fitted,  .resid,  .innov 값 

autoplot(aug, .innov) +
  labs(y = "$US",
       title = "Naïve 방법으로 예측했을떄의 잔차")

140일경에 잔차가 매우 특이하게 튄 것을 볼 수 있고 이후 2차례 정도 잔차의 변화가 큰 것을 식별할 수 있다. 이번에는 히스토그램을 통해 확인해보자.

aug %>%
  ggplot(aes(x = .innov)) +
  geom_histogram() +
  labs(title = "잔차 히스토그램 ")

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

히스토그램을 보면 평균이 0인 정규분포의 형상을 가지고 있지만 오른쪽으로 길게 늘어진 형태가 보인다. 이번에는 ACF 도표를 확인해 보자.

aug %>%
  ACF(.innov) %>% 
  autoplot() + 
  labs(title = "innvation 잔차 ACF 도표 ")

ACF 도표상으로 대부분의 lag에서 자기상관이 없는 것으로 확인된다. 이상의 잔차 진단을 종합해 보면 Naive 방법으로 모델을 적합 시켰을 때, innovation 잔차는 위 가정들을 만족한다고 볼 수 있다. 잔차를 진단하기 위해 여러 그래프를 한꺼번에 그려주는 함수인 gg_tsresiduals() 함수를 이용하면 전체를 쉽게 확인 할 수 있다. gg_tsresiduals() 함수는 tsibble 데이터를 적합 시킨 mable 객체에 적용시켜야 한다.

google_2015 %>%
  model(Naive = NAIVE(Close)) %>%
  gg_tsresiduals() +
  labs(title = "gg_tsresidual() 함수를 이용한 잔차 진단용 그래프 ")

검정을 통해 잔차진단을 엄격하게 하는 방법도 있다. 잔차를 검정하는 여러 검정 방법을 묶어 Portmanteau(가방이라는 뜻) 검정이라고 한다. 여기에는 Box-Pierce 검정, Ljung-Box 검정 등이 있다.

aug %>% 
  features(.innov, box_pierce, lag = 10)   # aug의 .innov를 box_pierce 로  

# A tibble: 1 × 4
  Symbol .model       bp_stat bp_pvalue
  <chr>  <chr>          <dbl>     <dbl>
1 GOOG   NAIVE(Close)    7.74     0.654

box_pierec 검정을 실시한 결과 p-value가 0.654로 상관관계가 0이라는 귀무설을 기각하지 못한다.

aug %>% 
  features(.innov, ljung_box, lag = 10)    # aug의 .innov를 ljung_box 로 

# A tibble: 1 × 4
  Symbol .model       lb_stat lb_pvalue
  <chr>  <chr>          <dbl>     <dbl>
1 GOOG   NAIVE(Close)    7.91     0.637

또한 ljung_box 검정 실시 결과 역시 p-value가 0.637로 자기상관이 0이라는 귀무가설을 기각하지 못한다.
따라서 잔차 진단을 수행해 본 결과, naive 모델로 적합 시킨 모델은 타당하다고 결론을 내릴 수 있다.

분해 예측(Decomposition)

시계열 예측 시리즈 II에서 논의한 시계열 데이터의 분해된 요소들로 예측을 수행할 수 있다. 시계열 데이터를 STL 함수로 분해하는 모델을 만들어 요소별 컴포넌트를 확인하면, 시계열 데이터의 원 value(아래 예제 소스에서 Employed에 해당)가 분해되어 추세성분(trend), 계절성분(season_year), 잔여성분(remainder)으로 분해되고 계절성분이 제거된 계절조정(season_adjust) 성분이 제시된다. 이를 가법모델 수식으로 표현하면, 계절 조성 성분은 \(A_t\)는 \(A_t = T_t + R_t\)이고, \(Y_t = \hat{T}_t + \hat{A}_t\) 이된다. 실제 계절성분은 계절별로 변화의 폭이 크지 않다. 아래 소스를 살펴보자.

us_retail_employment <- us_employment %>%
  filter(year(Month) >= 1990, Title == "Retail Trade")

dcmp <- us_retail_employment %>%
  model(STL(Employed ~ trend(window = 7), robust = TRUE)) %>%  # Employed를 STL로 성분 분해한 모델  
  components() %>%          # complonet별로 정리 --> dable이라고 함 
  select(-.model)           # 모델을 제거하고 나면 tsibble 객체가 됨 
  # dcmp를 dable에서 tsibble로 만들어 tidyverse하게 처리하기 위해 ???                               

성분 분해가 된 시계열 데이터(tsibble 객체)를 계절조정 요소(season_adjust)를 NAIVE() 방법으로 적합시킨 모델을 만들어 예측한다.

dcmp %>%                            # 요소별로 분해된 tsibble 
  model(NAIVE(season_adjust)) %>%   # season_adjust 값을 NAIVE로 적합하는 모델 
  forecast() %>%                    # 적합된 모델로 예측하면, 디폴트는 2년(h=24) 
  autoplot(dcmp) +                  # 원데이터와 forecast 결과를 도시 
  labs(y = "소매관련 고용인수 ",
       title = "미국 소매관련 고용인 예측 ")

이 과정은 아래와 같이 decomposition_model() 함수를 써서 분해와 적합을 동시에 수행하여 단순화 시킬 수 있다. 예측은 디폴트로 NAIVE 방법을 사용한다.

us_retail_employment  %>%
  model(stlf = decomposition_model(
    STL(Employed ~ trend(window = 7), robust = TRUE),
    NAIVE(season_adjust)
  ))   %>%
  forecast() |>   #SNAIVE()
  autoplot(us_retail_employment)+
  labs(y = "소매관련 고용인수",
       title = "미국 소매관련 고용인 예측")

예측 정확도 평가

잔차의 크기는 실제 예측 오류가 얼마나 클지 나타내는 신뢰할 만한 지표는 아니다. 예측의 정확성은 모델을 피팅할 때 사용되지 않은 새로운 데이터에 대해 모델이 얼마나 잘 예측하는 지를 고려해야 한다. 모델을 적합할 때는 훈련 데이터를 사용하고 정확도 평가할때는 테스트 데이터를 사용해야 한다. 테스트 데이터는 적어도 예측하려는 구간의 크기만큼은 되어야 한다.

데이터를 분할할때 사용하는 함수는 filter() 함수, filter_index()함수 또는 slice() 함수를 사용한다.

aus_production %>% filter(year(Quarter) >= 1995)  # 1995년 이후 데이터만 

aus_production %>% filter_index("1995 Q1" ~ .)    # 1995년 이후 데이터만

aus_production %>% slice(n()-19:0)  #  aus_production %>% tail(n=20)와 동일 

aus_retail %>% 
  group_by(State, Industry) %>%     # key 변수별로 그룹핑해서 
  slice(1:12)                       # 각각 1년치 (12개월)씩 slicing 한다. 

예측 정확도 지표

예측 오차를 표현하는 정확도 지표에는 다음과 같은 지표가 있다.

• MAE : Mean Absolute Error
• RMSE : Root Mean Square Error
• MAPE : Mean Absolute Percent Error

recent_production <- aus_production %>%
  filter(year(Quarter) >= 1992)
beer_train <- recent_production %>%
  filter(year(Quarter) <= 2007)       #1992년부터 2007년까지 데이터가 train 

beer_fit <- beer_train %>%            # 네가지 모델로 적합시킨 모델을 beer_fit
  model(                              # 네가지 모델로 적합 
    Mean = MEAN(Beer),                # 평균으로 적합 
    `Naïve` = NAIVE(Beer),            # NAIVE
    `Seasonal naïve` = SNAIVE(Beer),  # Seasonal NAIVE
    Drift = RW(Beer ~ drift())        # Drift 
  )

beer_fc <- beer_fit %>%
  forecast(h = 10)     # aus_procuction은 2010년 Q2까지고, train은 2007년 까지
                       # 따라서 h = 10이면  2008년,2009년, 2010 Q1, Q2 10분기임 

beer_fc %>%
  autoplot(
    aus_production %>% filter(year(Quarter) >= 1992),
    level = NULL
  ) +
  labs(
    y = "백만 리터",
    title = "호주의 분기별 맥주 생산량 "
  ) +
  guides(colour = guide_legend(title = "Forecast"))

정확도 지표를 보기 위해서는 예측을 위한 10개 분기 데이터를 적용했을 때, 정확도를 계산해야 한다. 정확도는 accuracy() 함수를 이용한다. accuracy(f, x)는 첫 번째 인자로 예측모델 객체(f), 두 번째 인자로 평가할 데이터의 수치형 벡터(x)를 제공해야 한다. 따라서 recent_production 으로 데이터를 주면, 예측용 10개 분기 데이터에 대한 여러 가지 정확도를 알아서 계산해서 제시해 준다.

accuracy(beer_fc, recent_production)   # 예측모델에,  recent_production 이후 데이터 

# A tibble: 4 × 10
  .model         .type    ME  RMSE   MAE    MPE  MAPE  MASE RMSSE    ACF1
  <chr>          <chr> <dbl> <dbl> <dbl>  <dbl> <dbl> <dbl> <dbl>   <dbl>
1 Drift          Test  -54.0  64.9  58.9 -13.6  14.6  4.12  3.87  -0.0741
2 Mean           Test  -13.8  38.4  34.8  -3.97  8.28 2.44  2.29  -0.0691
3 Naïve          Test  -51.4  62.7  57.4 -13.0  14.2  4.01  3.74  -0.0691
4 Seasonal naïve Test    5.2  14.3  13.4   1.15  3.17 0.937 0.853  0.132 

정확도 지표를 보면, Seasonal naïve 방법이 가장 정확한 것으로 나타났으며, 이는 직관적인 판단과 일치한다.

시계열 교차검증

시계열 교차 검증을 위해서 테스트 데이터와 예측용 데이터는 위 그림과 같이 학습용 데이터와 예측용 데이터를 반복적으로 생성하여 각각을 예측하고 예측 결과의 정확도는 각 예측 정확도를 평균하여 평가한다. 훈련데이터를 반복적으로 생성하는 함수는 stretch_tsibble() 함수이다. 아래 예를 참조하자.

# Time series cross-validation accuracy
google_2015_tr <- google_2015 %>%
  stretch_tsibble(.init = 3, .step = 1) %>%  # 최초 3개로 시작해서 하나씩 늘림 
  relocate(Date, Symbol, .id)   # id는 자동으로 생성된 몇번째 교차검증인지 표시 

google_2015_tr

# A tsibble: 31,875 x 10 [1]
# Key:       Symbol, .id [250]
   Date       Symbol   .id  Open  High   Low Close Adj_Close  Volume   Day
   <date>     <chr>  <int> <dbl> <dbl> <dbl> <dbl>     <dbl>   <dbl> <int>
 1 2015-01-02 GOOG       1  526.  528.  521.  522.      522. 1447600     1
 2 2015-01-05 GOOG       1  520.  521.  510.  511.      511. 2059800     2
 3 2015-01-06 GOOG       1  512.  513.  498.  499.      499. 2899900     3
 4 2015-01-02 GOOG       2  526.  528.  521.  522.      522. 1447600     1
 5 2015-01-05 GOOG       2  520.  521.  510.  511.      511. 2059800     2
 6 2015-01-06 GOOG       2  512.  513.  498.  499.      499. 2899900     3
 7 2015-01-07 GOOG       2  504.  504.  497.  498.      498. 2065100     4
 8 2015-01-02 GOOG       3  526.  528.  521.  522.      522. 1447600     1
 9 2015-01-05 GOOG       3  520.  521.  510.  511.      511. 2059800     2
10 2015-01-06 GOOG       3  512.  513.  498.  499.      499. 2899900     3
# ℹ 31,865 more rows

생성된 교차검증용 학습 데이터를 보면 최초 3개의 데이터로부터 시작하여 매번 1개 씩 추가해 가며 데이터를 생성해 나가는 것을 확인 할 수 있다. id는 몇 번째 교차검증용 데이터 인지를 나타낸다. 그럼 생성된 교차검증용 데이터로 학습을 시키고 예측을 수행해 보자.

# TSCV accuracy
google_2015_tr %>%                # 위에서 생성한 교차검증용 학습 데이터  
  model(RW = RW(Close ~ drift())) %>%  # Drift 모델(랜덤워크 모델)로 적합 
  forecast(h = 1) %>%             # h = 1로 예측
  accuracy(google_2015)           # 위에 google_2015_tr에는 예측의 기준 값은 없음 

# A tibble: 1 × 11
  .model Symbol .type    ME  RMSE   MAE   MPE  MAPE  MASE RMSSE   ACF1
  <chr>  <chr>  <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>  <dbl>
1 RW     GOOG   Test  0.726  11.3  7.26 0.112  1.19  1.02  1.01 0.0985

                                  # 예측하는 값은 google_2015에 들어있다 !!!!     
# Training set accuracy
google_2015 %>%                   # 비교를 위해서 그냥 데이터로 적합해서 
  model(RW(Close ~ drift())) %>%
  accuracy()

# A tibble: 1 × 11
  Symbol .model     .type        ME  RMSE   MAE     MPE  MAPE  MASE RMSSE   ACF1
  <chr>  <chr>      <chr>     <dbl> <dbl> <dbl>   <dbl> <dbl> <dbl> <dbl>  <dbl>
1 GOOG   RW(Close … Trai… -2.97e-14  11.1  7.16 -0.0267  1.18  1.00 0.996 0.0976

교차검증한 정확도 보다 전체 데이터로 적합 시킨 모델의 정확도가 더 높은 것을 알 수 있다. 교차검증은 일부 데이터만 가지고 예측한 것이기 때문에 전체 데이터로 적합시킨 정확도 보다 떨어지는 것은 당연하다. 교차 검증을 하는 이유는 정확도를 높이기 위해서가 아니라 과적합 문제를 해결하기 위해서 사용하는 기법이다.

아래는 예측 구간(h 참고:책에서는 h를 horizon이라고 표현)의 길이에 따라 정확도 지표인 RMSE값이 어떻게 달라지는지를 보여준다.

fc <- google_2015_tr %>%
  model(RW(Close ~ drift())) %>%
  forecast(h = 8) %>%           # 예측 구간을 8개 까지 
  group_by(.id) %>%             # 8개씩 예측했으니, id로 그룹화화
  mutate(h = row_number()) %>%  # 각 id별로 1부터 8까지 번호 부여 
  ungroup() %>%
  as_fable(response = "Close", distribution = Close) # fable로 변환 

fc %>%
  accuracy(google_2015, by = c("h", ".model")) %>%   
  ggplot(aes(x = h, y = RMSE)) +
  geom_point()

예측구간(h)가 늘어 날 수록 RMSE가 증하는 것을 확인 할 수 있다.