Discussion 9

Author

Ethan Wright

1.

library(fpp3)

Warning: package 'fpp3' was built under R version 4.5.2

Registered S3 method overwritten by 'tsibble':
  method               from 
  as_tibble.grouped_df dplyr

── Attaching packages ──────────────────────────────────────────── fpp3 1.0.2 ──

✔ tibble      3.3.0     ✔ tsibble     1.1.6
✔ dplyr       1.1.4     ✔ tsibbledata 0.4.1
✔ tidyr       1.3.1     ✔ feasts      0.4.2
✔ lubridate   1.9.4     ✔ fable       0.5.0
✔ ggplot2     3.5.2

Warning: package 'tsibble' was built under R version 4.5.2

Warning: package 'tsibbledata' was built under R version 4.5.2

Warning: package 'feasts' was built under R version 4.5.2

Warning: package 'fabletools' was built under R version 4.5.2

Warning: package 'fable' was built under R version 4.5.2

── Conflicts ───────────────────────────────────────────────── fpp3_conflicts ──
✖ lubridate::date()    masks base::date()
✖ dplyr::filter()      masks stats::filter()
✖ tsibble::intersect() masks base::intersect()
✖ tsibble::interval()  masks lubridate::interval()
✖ dplyr::lag()         masks stats::lag()
✖ tsibble::setdiff()   masks base::setdiff()
✖ tsibble::union()     masks base::union()

library(fredr)

Warning: package 'fredr' was built under R version 4.5.2

library(tidyverse)

Warning: package 'stringr' was built under R version 4.5.2

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ forcats 1.0.1     ✔ readr   2.1.5
✔ purrr   1.1.0     ✔ stringr 1.6.0

── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter()     masks stats::filter()
✖ tsibble::interval() masks lubridate::interval()
✖ dplyr::lag()        masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(patchwork)
library(feasts)
library(fabletools)
library(quantmod)

Loading required package: xts
Loading required package: zoo

Attaching package: 'zoo'

The following object is masked from 'package:tsibble':

    index

The following objects are masked from 'package:base':

    as.Date, as.Date.numeric


######################### Warning from 'xts' package ##########################
#                                                                             #
# The dplyr lag() function breaks how base R's lag() function is supposed to  #
# work, which breaks lag(my_xts). Calls to lag(my_xts) that you type or       #
# source() into this session won't work correctly.                            #
#                                                                             #
# Use stats::lag() to make sure you're not using dplyr::lag(), or you can add #
# conflictRules('dplyr', exclude = 'lag') to your .Rprofile to stop           #
# dplyr from breaking base R's lag() function.                                #
#                                                                             #
# Code in packages is not affected. It's protected by R's namespace mechanism #
# Set `options(xts.warn_dplyr_breaks_lag = FALSE)` to suppress this warning.  #
#                                                                             #
###############################################################################

Attaching package: 'xts'

The following objects are masked from 'package:dplyr':

    first, last

Loading required package: TTR
Registered S3 method overwritten by 'quantmod':
  method            from
  as.zoo.data.frame zoo

library(tseries)

fredr_set_key("18f0dfa970a5782481171e8fa3980e41")


cpi_raw <- fredr(series_id = "CPIAUCSL",
                 observation_start = as.Date("1990-01-01"),
                 observation_end   = as.Date("2023-12-01"))

oil_raw <- fredr(series_id = "MCOILWTICO",
                 observation_start = as.Date("1990-01-01"),
                 observation_end   = as.Date("2023-12-01"))

cpi <- cpi_raw |>
  mutate(Month = yearmonth(date)) |>
  as_tsibble(index = Month) |>
  select(Month, CPI = value)

oil <- oil_raw |>
  mutate(Month = yearmonth(date)) |>
  as_tsibble(index = Month) |>
  select(Month, Oil = value)

both <- cpi |>
  left_join(oil, by = "Month") |>
  drop_na()

both |> autoplot(CPI) +
  labs(title = "US CPI", y = "Index", x = NULL)

both |> autoplot(Oil) +
  labs(title = "WTI Crude Oil Price", y = "USD per barrel", x = NULL)

n_train <- floor(nrow(both) * 0.8)

train <- both |> slice(1:n_train)
test  <- both |> slice((n_train + 1):nrow(both))

fit <- train |>
  model(
    ETS   = ETS(CPI),
    ARIMA = ARIMA(CPI)
  )

glance(fit) |> select(.model, AIC, AICc, BIC)

# A tibble: 2 × 4
  .model   AIC  AICc   BIC
  <chr>  <dbl> <dbl> <dbl>
1 ETS    1418. 1418. 1440.
2 ARIMA   421.  422.  444.

tidy(fit)

# A tibble: 10 × 6
   .model term     estimate std.error statistic   p.value
   <chr>  <chr>       <dbl>     <dbl>     <dbl>     <dbl>
 1 ETS    alpha      1.000    NA          NA    NA       
 2 ETS    beta       0.112    NA          NA    NA       
 3 ETS    phi        0.980    NA          NA    NA       
 4 ETS    l[0]     127.       NA          NA    NA       
 5 ETS    b[0]       0.660    NA          NA    NA       
 6 ARIMA  ma1        0.547     0.0554      9.87  2.95e-20
 7 ARIMA  ma2        0.0627    0.0625      1.00  3.16e- 1
 8 ARIMA  ma3       -0.116     0.0564     -2.07  3.97e- 2
 9 ARIMA  sma1      -0.194     0.0543     -3.58  3.97e- 4
10 ARIMA  constant   0.355     0.0305     11.6   2.15e-26

fc <- fit |>
  forecast(h = nrow(test))

fc |>
  accuracy(both) |>
  select(.model, RMSE, MAE, MAPE, MASE)

# A tibble: 2 × 5
  .model  RMSE   MAE  MAPE  MASE
  <chr>  <dbl> <dbl> <dbl> <dbl>
1 ARIMA   17.0  11.0  3.74  2.52
2 ETS     23.0  15.3  5.27  3.53

fc |>
  autoplot(both, level = 95) +
  facet_wrap(~ .model, ncol = 1) +
  labs(title    = "ETS vs ARIMA — 80/20 Test Split",
       subtitle  = "Vertical line = train/test cutoff",
       x = NULL, y = "CPI Index") +
  geom_vline(xintercept = as.numeric(yearmonth("2017 Feb")),
             linetype = "dashed", colour = "black")

No, although both models do no not forecast the positive shock, ARIMA performs better across all accuracy metrics and has a tighter confidence interval.

2.

fit_all <- train |>
  model(
    ETS = ETS(CPI),
    ARIMA  = ARIMA(CPI),
    Dynamic = ARIMA(CPI ~ Oil)
  )

glance(fit_all) |> select(.model, AIC, AICc, BIC)

# A tibble: 3 × 4
  .model    AIC  AICc   BIC
  <chr>   <dbl> <dbl> <dbl>
1 ETS     1418. 1418. 1440.
2 ARIMA    421.  422.  444.
3 Dynamic  340.  341.  359.

tidy(fit_all)

# A tibble: 14 × 6
   .model  term      estimate std.error statistic   p.value
   <chr>   <chr>        <dbl>     <dbl>     <dbl>     <dbl>
 1 ETS     alpha       1.000   NA           NA    NA       
 2 ETS     beta        0.112   NA           NA    NA       
 3 ETS     phi         0.980   NA           NA    NA       
 4 ETS     l[0]      127.      NA           NA    NA       
 5 ETS     b[0]        0.660   NA           NA    NA       
 6 ARIMA   ma1         0.547    0.0554       9.87  2.95e-20
 7 ARIMA   ma2         0.0627   0.0625       1.00  3.16e- 1
 8 ARIMA   ma3        -0.116    0.0564      -2.07  3.97e- 2
 9 ARIMA   sma1       -0.194    0.0543      -3.58  3.97e- 4
10 ARIMA   constant    0.355    0.0305      11.6   2.15e-26
11 Dynamic ma1         0.287    0.0608       4.71  3.60e- 6
12 Dynamic sma1       -0.175    0.0533      -3.29  1.13e- 3
13 Dynamic Oil         0.0559   0.00582      9.60  2.19e-19
14 Dynamic intercept   0.351    0.0238      14.7   6.10e-38

fc_all <- fit_all |>
  forecast(new_data = test)

fc_all |>
  accuracy(both) |>
  select(.model, RMSE, MAE, MAPE, MASE)

# A tibble: 3 × 5
  .model   RMSE   MAE  MAPE  MASE
  <chr>   <dbl> <dbl> <dbl> <dbl>
1 ARIMA    17.0  11.0  3.74  2.52
2 Dynamic  16.0  10.1  3.45  2.33
3 ETS      23.0  15.3  5.27  3.53

fc_all |>
  autoplot(both, level = 95) +
  facet_wrap(~ .model, ncol = 1) +
  geom_vline(xintercept = as.numeric(yearmonth("2017 Feb")),
             linetype = "dashed", colour = "black") +
  labs(title    = "ETS vs ARIMA vs Dynamic Regression — 80/20 Test Split",
       subtitle  = "Dashed line = train/test cutoff",
       x = NULL, y = "CPI Index")

Out of the 3 models (ETS, ARIMA, ARIMA with Oil regressor), the dynamic regression (CPI~Oil) ARIMA(0,1,1)(0,1,1)[12] wins out on every metric, and also has a tighter confidence interval than the basic ARIMA model. So, oil price, which can push inflation as oil is a raw input in many goods, is a decent regressor.