ETS

library(tidyverse)
library(fpp3)

1 Data

1.1 Data Import

We will download data from FRED regarding the 10-year real interest rate (REAINTRATREARAT10Y), and the employment in Private Education and Health Services (CEU6500000001).

datos_tbl <- tidyquant::tq_get(
  x    = c("REAINTRATREARAT10Y", "CEU6500000001"), 
  get  = "economic.data",
  from = "1982-01-01",
  to   = "2023-02-01"
)

Registered S3 method overwritten by 'quantmod':
  method            from
  as.zoo.data.frame zoo 

datos_tbl

1.2 Data Wrangle

Convert to tsibble.

datos_tsbl <- datos_tbl |> 
  mutate(
    date = yearmonth(date),
    symbol = if_else(symbol == "REAINTRATREARAT10Y",
                     "int_rate", "employment")
  ) |> 
  as_tsibble(
    index = date,
    key   = symbol
  )

datos_tsbl

2 Train & Test Splits

We will set the last two years as test data.

datos_train <- datos_tsbl |> 
  filter_index(. ~ "2021 feb.")

3 Visualization

datos_train |> 
  autoplot(price) +
  facet_wrap(~ symbol, scales = "free_y", ncol = 1) +
  theme(legend.position = "none")

datos_train |> 
  gg_season(price)

We noticed that the employment time series has both seasonality and an upward trend, whereas the interest rate time series has only a downward trend. Thus, we will need different models for each.

datos_wide <- datos_train |> 
  pivot_wider(
    names_from  = symbol,
    values_from = price
  )
datos_wide

datos_tsbl_wide <- datos_tsbl |> 
  pivot_wider(
    names_from  = symbol,
    values_from = price
  )

4 Employment

4.1 Model

employment_fit <- datos_wide |> 
  model(
    drift    = RW(employment ~ drift()),
    snaive   = SNAIVE(employment),
    ets_aada = ETS(employment ~ error("A") + trend("Ad") + season("A")),
    ets_madm = ETS(employment ~ error("M") + trend("Ad") + season("M"))
  )
employment_fit

4.2 Residual Diagnostics

employment_fit |> 
  select(ets_aada) |> 
  gg_tsresiduals() +
  ggtitle("Residual Diagnostics for the ETS(A,Ad,A) model")

employment_fit |> 
  select(ets_madm) |> 
  gg_tsresiduals() +
  ggtitle("Residual Diagnostics for the ETS(M,Ad,M) model")

Both ETS models show still some patterns on their residuals. Hence, we will try to improve them.

We will focus on the MAE error metric for this analysis.

accuracy(employment_fit) |> 
  select(.model, .type, MAE, RMSE, MAPE, MASE) |> 
  arrange(MAE)

lambda <- datos_wide |> 
  features(employment, guerrero) |> 
  pull(lambda_guerrero)

employment_fit2 <- datos_wide |> 
  model(
    ets_bc_aada = ETS(box_cox(employment, lambda) ~ error("A") + trend("Ad") + season("A")),
    ets_bc_madm = ETS(box_cox(employment, lambda) ~ error("M") + trend("Ad") + season("M"))
  )

employment_fit2 |> 
  select(ets_bc_aada) |> 
  gg_tsresiduals() +
  ggtitle("Residual Diagnostics for the ETS(A,Ad,A) model using a Box-Cox transformation")

employment_fit2 |> 
  select(ets_bc_madm) |> 
  gg_tsresiduals() +
  ggtitle("Residual Diagnostics for the ETS(M,Ad,M) model using a Box-Cox transformation")

It doesn’t seem to improve the fit.

Let’s try using a decomposition model.

datos_wide |> 
  model(
    stl = STL(employment, robust = TRUE)
  ) |> 
  components() |> 
  autoplot()

employment_fit3 <- datos_wide |> 
  model(
    stl_drift_snaive = decomposition_model(
      STL(employment, robust = TRUE),
      SNAIVE(season_year),
      RW(season_adjust ~ drift())
    ),
    stl_ets = decomposition_model(
      STL(employment, robust = TRUE),
      ETS(season_year ~ error("M") + trend("N") + season("M")),
      ETS(season_adjust ~ error("A") + trend("A") + season("N"))
    )
  )

employment_fit3

accuracy(employment_fit3) |> 
  select(.model, .type, MAE, RMSE, MAPE, MASE) |> 
  arrange(MAE)

4.3 Forecasting on the test set

employment_fit_full <- employment_fit |> 
  cross_join(employment_fit2) |> 
  cross_join(employment_fit3)

employment_fc <- employment_fit_full |> 
  forecast(h = "2 years")

employment_fc

employment_fc |> 
  autoplot(datos_tsbl_wide |> filter_index("2010 jan." ~ .), level = NULL)

employment_fc |> 
  filter(.model %in% c("drift", "stl_drift_snaive", "stl_ets")) |> 
  autoplot(datos_tsbl_wide |> filter_index("2010 jan." ~ .), level = NULL)

employment_fc |> 
  accuracy(datos_tsbl_wide) |> 
  select(.model, .type, MAE, RMSE, MAPE, MASE) |> 
  arrange(MAE)

employment_fit_combination <- employment_fit_full |> 
  mutate(combination = (drift + stl_ets + stl_drift_snaive)/3) |> 
  select(contains("drift"), contains("stl"), combination) 

employment_fc2 <- employment_fit_combination |> forecast(h = "2 years")

employment_fc2 |> 
  autoplot(datos_tsbl_wide |> filter_index("2015 jan." ~ .), level = NULL)

employment_fc2 |> 
  accuracy(datos_tsbl_wide) |> 
  select(.model, .type, MAE, RMSE, MAPE, MASE) |> 
  arrange(MAE)

4.4 Forecasting the future

employment_refit <- datos_tsbl_wide |> 
  drop_na() |> 
  model(
    combination = combination_model(
      RW(employment ~ drift()),
      decomposition_model(
      STL(employment, robust = TRUE),
      SNAIVE(season_year),
      RW(season_adjust ~ drift())
    ),
    decomposition_model(
      STL(employment, robust = TRUE),
      ETS(season_year ~ error("M") + trend("N") + season("M")),
      ETS(season_adjust ~ error("A") + trend("A") + season("N"))
    )
    )
  )
employment_refit

employment_refit |> 
  forecast(h = "2 years") |> 
  autoplot(datos_tsbl_wide |> filter_index("2015 jan." ~ .))