library(tidyverse)
library(fpp3)
library(ggthemes)
theme_set(theme_calc())
theme_update(axis.title = element_blank())Homework 4
Read data
add_rec_shade <- function(st_date, ed_date, shade_color = "darkgray") {
library(fredr)
library(ecm)
library(ggplot2)
# fredr_set_key("insert_your_api_key")
recession <- fredr(series_id = "USRECD", observation_start = as.Date(st_date), observation_end = as.Date(ed_date))
recession$diff <- recession$value - lagpad(recession$value, k = 1)
recession <- recession[!is.na(recession$diff), ]
recession.start <- recession[recession$diff == 1, ]$date
recession.end <- recession[recession$diff == (-1), ]$date
if (length(recession.start) > length(recession.end)) {
recession.end <- c(recession.end, Sys.Date())
}
if (length(recession.end) > length(recession.start)) {
recession.start <- c(min(recession$date), recession.start)
}
recs <- as.data.frame(cbind(recession.start, recession.end))
recs$recession.start <- as.Date(as.numeric(recs$recession.start), origin = as.Date("1970-01-01"))
recs$recession.end <- as.Date(recs$recession.end, origin = as.Date("1970-01-01"))
if (nrow(recs) > 0) {
rec_shade <- geom_rect(
data = recs, inherit.aes = F,
aes(xmin = recession.start, xmax = recession.end, ymin = -Inf, ymax = +Inf),
fill = shade_color, alpha = 0.5
)
return(rec_shade)
}
}Advanced retail sales: clothing stores
library(tidyverse)
library(fpp3)
library(fredr)
library(ggthemes)
clothing_sales <- fredr(
series_id = "RSCCASN",
observation_start = as.Date("1992-01-01"),
observation_end = as.Date("2023-10-01"),
frequency = "m",
units = "lin"
) |>
select(date, value) |>
mutate(date = tsibble::yearmonth(date)) |>
as_tsibble()Analyze seasonal periods
clothing_sales |>
model(stl = STL(value)) |>
components() |>
autoplot()
Increasing seasonal pattern indicates a need for a box-cox transformation
lambda <- clothing_sales |>
features(value, features = guerrero) |>
pull(lambda_guerrero)clothing_sales |>
autoplot(box_cox(value, lambda)) +
labs(
title = "Monthly seasonally unadjusted retail clothing sales, Box-Cox transformed"
) +
theme_calc() +
theme(axis.title = element_blank())
Training/test splits
train <- clothing_sales |>
filter_index("1992 Jan" ~ "2017 May")
test <- clothing_sales |>
filter_index("2017 June" ~ .)Fit training model
training_fit <- train |>
model(
add_ets = ETS(log(value) ~ error("A") + trend("A") + season("A")),
mult_ets = ETS(log(value) ~ error("M") + trend("A") + season("M")),
auto_arima = ARIMA(log(value), stepwise = FALSE, approximation = FALSE),
dhr = ARIMA(log(value) ~ fourier(K = 6))
)training_fit |> report()# A tibble: 4 × 11
.model sigma2 log_lik AIC AICc BIC MSE AMSE MAE
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 add_ets 0.000914 203. -372. -370. -309. 0.000866 0.000914 0.0228
2 mult_ets 0.0000102 198. -362. -360. -299. 0.000892 0.000945 0.00238
3 auto_arima 0.000848 621. -1231. -1231. -1209. NA NA NA
4 dhr 0.000890 639. -1247. -1245. -1187. NA NA NA
# ℹ 2 more variables: ar_roots <list>, ma_roots <list>Forecast training model
training_fc <- training_fit |>
forecast(h = 12)Check against test set
All models are performing well in terms of forecasting seasonal period (December)
training_fc |>
autoplot(
clothing_sales |>
filter_index("2014 Jan" ~ "2018 May")
)
Accuracy metrics: Additive ETS has lowest RMSE, Multiplicative ETS has lowest MAPE
accuracy(training_fc, clothing_sales)# 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 add_ets Test -30.4 537. 437. -0.219 1.99 0.647 0.636 -0.521
2 auto_arima Test -337. 688. 545. -1.60 2.50 0.806 0.815 -0.427
3 dhr Test -131. 755. 568. -0.547 2.43 0.840 0.895 -0.371
4 mult_ets Test 73.5 553. 421. 0.283 1.86 0.624 0.656 -0.567Forecst on test set
test_fit <- test |>
model(
add_ets = ETS(log(value) ~ error("A") + trend("A") + season("A")),
mult_ets = ETS(log(value) ~ error("M") + trend("A") + season("M")),
auto_arima = ARIMA(log(value), stepwise = FALSE, approximation = FALSE),
dhr = ARIMA(log(value) ~ fourier(K = 6))
) |>
mutate(ensemble = (add_ets + mult_ets + auto_arima + dhr) / 4)test_fit |>
select(!ensemble) |>
glance() |>
arrange(AIC)# A tibble: 4 × 11
.model sigma2 log_lik AIC AICc BIC MSE AMSE MAE ar_roots
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <list>
1 dhr 0.0409 20.8 -11.6 -3.63 23.3 NA NA NA <cpl>
2 auto_arima 0.0651 -4.40 20.8 22.0 34.8 NA NA NA <cpl>
3 add_ets 0.0649 -51.7 137. 148. 177. 0.0512 0.112 0.115 <NULL>
4 mult_ets 0.00103 -69.3 173. 183. 212. 0.0784 0.0803 0.0146 <NULL>
# ℹ 1 more variable: ma_roots <list>Auto ARIMA has the lowest RMSE on the test data
training_fit |>
refit(test) |>
accuracy()# 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 add_ets Training 22547. 23260. 22547. 99.9 99.9 7.09 4.14 0.336
2 mult_ets Training 22547. 23260. 22547. 99.9 99.9 7.09 4.14 0.336
3 auto_arima Training -71.7 4514. 2381. -7.56 18.3 0.749 0.803 0.550
4 dhr Training -229. 4744. 2513. -7.60 18.5 0.791 0.844 0.5611-year ahead forecast
test_fc <- test_fit |>
forecast(h = 12)test_fc |>
autoplot(
clothing_sales |>
filter_index("2018 Jan" ~.),
leves = c(90, 95)
) +
scale_x_yearmonth(
breaks = scales::breaks_width("6 months"),
labels = scales::label_date_short()
) +
scale_y_continuous(
limit = c(0, 100000)
) +
labs(
title = "Forecasted monthly sales, clothing retailers"
) +
add_rec_shade("2020-01-01", "2023-10-01")
Plot each models residuals
models <- names(test_fit)
reg_plots <- purrr::map(models, ~ {
test_fit |>
select(.x) |>
gg_tsresiduals()
})
Overall, the best performing model was auto ARIMA, based on RMSE.