8.1

library(fpp3)
## Registered S3 method overwritten by 'tsibble':
##   method               from 
##   as_tibble.grouped_df dplyr
## ── Attaching packages ──────────────────────────────────────────── fpp3 1.0.2 ──
## ✔ tibble      3.3.1     ✔ tsibble     1.1.6
## ✔ dplyr       1.2.0     ✔ tsibbledata 0.4.1
## ✔ tidyr       1.3.2     ✔ feasts      0.4.2
## ✔ lubridate   1.9.5     ✔ fable       0.5.0
## ✔ ggplot2     4.0.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()
# 1) Data: pigs slaughtered in Victoria
pigs_vic <- aus_livestock |>
  filter(State == "Victoria", Animal == "Pigs") |>
  select(Month, Count)

# Optional quick plot
autoplot(pigs_vic, Count) +
  labs(title = "Pigs slaughtered in Victoria", y = "Count", x = "Month")

# 2) Fit SES using ETS(A,N,N)
fit_ses <- pigs_vic |>
  model(ETS(Count ~ error("A") + trend("N") + season("N")))

# Show model details (includes alpha and initial level l)
report(fit_ses)
## Series: Count 
## Model: ETS(A,N,N) 
##   Smoothing parameters:
##     alpha = 0.3221247 
## 
##   Initial states:
##      l[0]
##  100646.6
## 
##   sigma^2:  87480760
## 
##      AIC     AICc      BIC 
## 13737.10 13737.14 13750.07
# Extract alpha and l0 programmatically (works with fable/ets objects)
fit_tbl <- fit_ses |>
  tidy()

alpha_hat <- fit_tbl |>
  filter(term == "alpha") |>
  pull(estimate)

l0_hat <- fit_tbl |>
  filter(term %in% c("l", "l0")) |>
  pull(estimate)

cat("Estimated alpha:", alpha_hat, "\n")
## Estimated alpha: 0.3221247
cat("Estimated l0:", l0_hat, "\n")
## Estimated l0:
# 3) Forecast next 4 months
fc4 <- fit_ses |>
  forecast(h = 4)

# Print forecasts
fc4
## # A fable: 4 x 4 [1M]
## # Key:     .model [1]
##   .model                                                        Month
##   <chr>                                                         <mth>
## 1 "ETS(Count ~ error(\"A\") + trend(\"N\") + season(\"N\"))" 2019 Jan
## 2 "ETS(Count ~ error(\"A\") + trend(\"N\") + season(\"N\"))" 2019 Feb
## 3 "ETS(Count ~ error(\"A\") + trend(\"N\") + season(\"N\"))" 2019 Mar
## 4 "ETS(Count ~ error(\"A\") + trend(\"N\") + season(\"N\"))" 2019 Apr
## # ℹ 2 more variables: Count <dist>, .mean <dbl>
# Plot forecasts
autoplot(fc4, pigs_vic) +
  labs(title = "SES forecasts (ETS(A,N,N)) — next 4 months",
       y = "Count", x = "Month")

# 4) Residual SD (s)
resid_sd <- fit_ses |>
  augment() |>
  summarise(s = sd(.resid, na.rm = TRUE)) |>
  pull(s)

cat("Residual SD (s):", resid_sd, "\n")
## Residual SD (s): NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
# 5) Manual 95% PI for the first forecast: yhat ± 1.96*s
first_mean <- fc4 |>
  as_tibble() |>
  slice(1) |>
  pull(.mean)

lower_manual <- first_mean - 1.96 * resid_sd
upper_manual <- first_mean + 1.96 * resid_sd

cat("First forecast mean:", first_mean, "\n")
## First forecast mean: 95186.56
cat("Manual 95% PI (first forecast): [", lower_manual, ", ", upper_manual, "]\n")
## Manual 95% PI (first forecast): [ NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA ,  NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA ]
# 6) R-produced 95% interval (via hilo)
pi_r <- fc4 |>
  hilo(level = 95) |>
  as_tibble() |>
  select(Month, `95%`)

pi_r
## # A tibble: 4 × 2
##      Month                  `95%`
##      <mth>                 <hilo>
## 1 2019 Jan [76854.79, 113518.3]95
## 2 2019 Feb [75927.17, 114445.9]95
## 3 2019 Mar [75042.22, 115330.9]95
## 4 2019 Apr [74194.54, 116178.6]95
# If you want to see the numeric lower/upper bounds explicitly:
pi_r_bounds <- fc4 |>
  hilo(level = 95) |>
  unpack_hilo(`95%`) |>
  as_tibble() |>
  select(Month, .mean, `95%_lower`, `95%_upper`)

pi_r_bounds
## # A tibble: 4 × 4
##      Month  .mean `95%_lower` `95%_upper`
##      <mth>  <dbl>       <dbl>       <dbl>
## 1 2019 Jan 95187.      76855.     113518.
## 2 2019 Feb 95187.      75927.     114446.
## 3 2019 Mar 95187.      75042.     115331.
## 4 2019 Apr 95187.      74195.     116179.
# Compare manual PI vs R PI for first forecast
first_r <- pi_r_bounds |>
  slice(1)

cat("\n--- Comparison (first forecast) ---\n")
## 
## --- Comparison (first forecast) ---
cat("Manual 95% PI: [", lower_manual, ", ", upper_manual, "]\n")
## Manual 95% PI: [ NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA ,  NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA ]
cat("R 95% PI:      [", first_r$`95%_lower`, ", ", first_r$`95%_upper`, "]\n")
## R 95% PI:      [ 76854.79 ,  113518.3 ]

8.5

library(fpp3)

# 0) Choose ONE country
country_pick <- "Australia"   # <- change if you want (e.g., "United States", "China", "India")

# 1) Build dataset: annual Exports for selected country
exports_cty <- global_economy |>
  filter(Country == country_pick) |>
  select(Year, Exports)

# Quick check
glimpse(exports_cty)
## Rows: 58
## Columns: 2
## $ Year    <dbl> 1960, 1961, 1962, 1963, 1964, 1965, 1966, 1967, 1968, 1969, 19…
## $ Exports <dbl> 12.99445, 12.40310, 13.94301, 13.00589, 14.93825, 13.22018, 12…
# 2) Plot the Exports series + discuss features (write your discussion under the plot in your report)
autoplot(exports_cty, Exports) +
  labs(
    title = paste("Annual Exports —", country_pick),
    x = "Year",
    y = "Exports"
  )

#  If you want to see growth/volatility a bit clearer, try log scale:
# autoplot(exports_cty, log(Exports)) +
#   labs(title = paste("Log Annual Exports —", country_pick), y = "log(Exports)")

# 3) Fit ETS(A,N,N) and ETS(A,A,N)
fit_ets <- exports_cty |>
  model(
    ETS_ANN = ETS(Exports ~ error("A") + trend("N") + season("N")),
    ETS_AAN = ETS(Exports ~ error("A") + trend("A") + season("N"))
  )

# View estimated parameters / model summary
report(fit_ets)
## Warning in report.mdl_df(fit_ets): Model reporting is only supported for
## individual models, so a glance will be shown. To see the report for a specific
## model, use `select()` and `filter()` to identify a single model.
## # A tibble: 2 × 9
##   .model  sigma2 log_lik   AIC  AICc   BIC   MSE  AMSE   MAE
##   <chr>    <dbl>   <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 ETS_ANN   1.36   -126.  257.  258.  264.  1.32  1.79 0.914
## 2 ETS_AAN   1.34   -124.  258.  259.  269.  1.25  1.59 0.893
# 4) Forecast and plot (choose horizon)
h <- 10  # next 10 years (change if your assignment specifies something else)

fc <- fit_ets |>
  forecast(h = h)

# Plot both sets of forecasts (faceted)
autoplot(fc, exports_cty) +
  labs(
    title = paste("Exports forecasts —", country_pick),
    x = "Year",
    y = "Exports"
  )

# overlay both forecast means in one plot
fc_means <- fc |>
  as_tibble() |>
  select(.model, Year, .mean)

ggplot() +
  geom_line(data = exports_cty, aes(x = Year, y = Exports)) +
  geom_line(data = fc_means, aes(x = Year, y = .mean, linetype = .model)) +
  labs(
    title = paste("Forecast means comparison —", country_pick),
    x = "Year",
    y = "Exports",
    linetype = "Model"
  )

# 5) Compute RMSE on TRAINING data
# accuracy(fit_ets) gives in-sample accuracy for each model
acc_train <- accuracy(fit_ets)

acc_train
## # A tibble: 2 × 10
##   .model  .type              ME  RMSE   MAE    MPE  MAPE  MASE RMSSE   ACF1
##   <chr>   <chr>           <dbl> <dbl> <dbl>  <dbl> <dbl> <dbl> <dbl>  <dbl>
## 1 ETS_ANN Training  0.232        1.15 0.914  1.09   5.41 0.928 0.928 0.0125
## 2 ETS_AAN Training -0.000000746  1.12 0.893 -0.387  5.39 0.907 0.904 0.109
rmse_ann <- acc_train |>
  filter(.model == "ETS_ANN") |>
  pull(RMSE)

rmse_aan <- acc_train |>
  filter(.model == "ETS_AAN") |>
  pull(RMSE)

cat("Training RMSE ETS(A,N,N):", rmse_ann, "\n")
## Training RMSE ETS(A,N,N): 1.146794
cat("Training RMSE ETS(A,A,N):", rmse_aan, "\n")
## Training RMSE ETS(A,A,N): 1.116727
# 6) Compare forecasts from both methods (table)
fc_table <- fc |>
  as_tibble() |>
  select(.model, Year, .mean)

fc_table
## # A tibble: 20 × 3
##    .model   Year .mean
##    <chr>   <dbl> <dbl>
##  1 ETS_ANN  2018  20.6
##  2 ETS_ANN  2019  20.6
##  3 ETS_ANN  2020  20.6
##  4 ETS_ANN  2021  20.6
##  5 ETS_ANN  2022  20.6
##  6 ETS_ANN  2023  20.6
##  7 ETS_ANN  2024  20.6
##  8 ETS_ANN  2025  20.6
##  9 ETS_ANN  2026  20.6
## 10 ETS_ANN  2027  20.6
## 11 ETS_AAN  2018  20.8
## 12 ETS_AAN  2019  21.0
## 13 ETS_AAN  2020  21.1
## 14 ETS_AAN  2021  21.3
## 15 ETS_AAN  2022  21.4
## 16 ETS_AAN  2023  21.5
## 17 ETS_AAN  2024  21.7
## 18 ETS_AAN  2025  21.8
## 19 ETS_AAN  2026  21.9
## 20 ETS_AAN  2027  22.1
# 7) 95% prediction interval for FIRST forecast using RMSE and normal errors
# Manual PI: yhat ± 1.96 * RMSE
first_fc <- fc |>
  as_tibble() |>
  group_by(.model) |>
  slice(1) |>
  ungroup() |>
  select(.model, Year, .mean)

first_fc
## # A tibble: 2 × 3
##   .model   Year .mean
##   <chr>   <dbl> <dbl>
## 1 ETS_AAN  2018  20.8
## 2 ETS_ANN  2018  20.6
# Manual PI for each model
manual_pi <- first_fc |>
  mutate(
    RMSE = if_else(.model == "ETS_ANN", rmse_ann, rmse_aan),
    lower_manual_95 = .mean - 1.96 * RMSE,
    upper_manual_95 = .mean + 1.96 * RMSE
  )

manual_pi
## # A tibble: 2 × 6
##   .model   Year .mean  RMSE lower_manual_95 upper_manual_95
##   <chr>   <dbl> <dbl> <dbl>           <dbl>           <dbl>
## 1 ETS_AAN  2018  20.8  1.12            18.6            23.0
## 2 ETS_ANN  2018  20.6  1.15            18.4            22.9
# 8) R-produced 95% intervals (model-based) using hilo()
r_pi <- fc |>
  hilo(level = 95) |>
  unpack_hilo(`95%`) |>
  as_tibble() |>
  group_by(.model) |>
  slice(1) |>
  ungroup() |>
  select(.model, Year, .mean, r_lower_95 = `95%_lower`, r_upper_95 = `95%_upper`)

r_pi
## # A tibble: 2 × 5
##   .model   Year .mean r_lower_95 r_upper_95
##   <chr>   <dbl> <dbl>      <dbl>      <dbl>
## 1 ETS_AAN  2018  20.8       18.6       23.1
## 2 ETS_ANN  2018  20.6       18.3       22.9
# 9) Side-by-side comparison: Manual PI vs R PI for first forecast
pi_compare <- manual_pi |>
  left_join(r_pi, by = c(".model", "Year", ".mean")) |>
  select(.model, Year, .mean,
         RMSE,
         lower_manual_95, upper_manual_95,
         r_lower_95, r_upper_95)

pi_compare
## # A tibble: 2 × 8
##   .model  Year .mean  RMSE lower_manual_95 upper_manual_95 r_lower_95 r_upper_95
##   <chr>  <dbl> <dbl> <dbl>           <dbl>           <dbl>      <dbl>      <dbl>
## 1 ETS_A…  2018  20.8  1.12            18.6            23.0       18.6       23.1
## 2 ETS_A…  2018  20.6  1.15            18.4            22.9       18.3       22.9
# ---------------------------
#
# - Main features: describe overall level, trend (up/down), volatility changes, any sudden shifts.
# - ETS(A,N,N) assumes no trend; ETS(A,A,N) allows additive trend (one extra state/parameter).
# - Merits:
#    * If data has strong trend, ETS(A,A,N) often fits/forecasts better.
#    * If trend is weak/unstable, ETS(A,N,N) can be more robust (less overfitting).
# - Compare forecast reasonableness: do forecasts follow recent direction without exploding?
# - Prediction intervals:
#    * Manual RMSE-based PI uses constant error size.
#    * R’s ETS intervals come from the model’s estimated error variance/state uncertainty,
#      so they can differ (often widen with horizon, and differ across models).
# ---------------------------

8.6

# 1) Extract Chinese GDP
china_gdp <- global_economy |>
  filter(Country == "China") |>
  select(Year, GDP)

# Plot original series
autoplot(china_gdp, GDP) +
  labs(title = "China GDP", y = "GDP", x = "Year")

# 2) Fit multiple ETS models
fit_china <- china_gdp |>
  model(
    ETS_auto      = ETS(GDP),                              # Automatic selection
    ETS_AAN       = ETS(GDP ~ error("A") + trend("A") + season("N")),
    ETS_damped    = ETS(GDP ~ error("A") + trend("Ad") + season("N")),
    ETS_boxcox    = ETS(box_cox(GDP, lambda = 0.2)),
    ETS_log       = ETS(box_cox(GDP, lambda = 0))          # log transform
  )

report(fit_china)
## Warning in report.mdl_df(fit_china): Model reporting is only supported for
## individual models, so a glance will be shown. To see the report for a specific
## model, use `select()` and `filter()` to identify a single model.
## # A tibble: 5 × 9
##   .model       sigma2 log_lik    AIC   AICc    BIC      MSE     AMSE      MAE
##   <chr>         <dbl>   <dbl>  <dbl>  <dbl>  <dbl>    <dbl>    <dbl>    <dbl>
## 1 ETS_auto   1.08e- 2 -1546.  3102.  3103.  3112.  4.00e+22 1.61e+23 7.93e- 2
## 2 ETS_AAN    3.88e+22 -1624.  3258.  3259.  3268.  3.61e+22 1.31e+23 9.59e+10
## 3 ETS_damped 3.96e+22 -1624.  3260.  3262.  3273.  3.62e+22 1.30e+23 9.49e+10
## 4 ETS_boxcox 3.94e+ 2  -289.   588.   589.   598.  3.67e+ 2 1.32e+ 3 1.58e+ 1
## 5 ETS_log    8.81e- 3    21.5  -33.1  -31.9  -22.8 8.20e- 3 2.30e- 2 7.22e- 2
# 3) Use large forecast horizon to see differences clearly
h <- 30

fc_china <- fit_china |>
  forecast(h = h)

# 4) Plot forecasts (faceted)
autoplot(fc_china, china_gdp) +
  labs(title = "Chinese GDP Forecasts — ETS Variations",
       y = "GDP",
       x = "Year")

# 5) Overlay forecast means in one comparison plot
fc_means <- fc_china |>
  as_tibble() |>
  select(.model, Year, .mean)

ggplot() +
  geom_line(data = china_gdp,
            aes(x = Year, y = GDP)) +
  geom_line(data = fc_means,
            aes(x = Year, y = .mean, color = .model)) +
  labs(title = "Comparison of ETS Forecast Means",
       y = "GDP",
       x = "Year",
       color = "Model")

## 8.7

# ---------------------------
# Gas production — ETS modeling
# ---------------------------

# 1) Extract Gas data (Quarterly)
gas_data <- aus_production |>
  select(Quarter, Gas)

# Plot series
autoplot(gas_data, Gas) +
  labs(title = "Australian Gas Production",
       y = "Gas Production",
       x = "Quarter")

# 2) Fit automatic ETS model
fit_auto <- gas_data |>
  model(ETS(Gas))

report(fit_auto)
## Series: Gas 
## Model: ETS(M,A,M) 
##   Smoothing parameters:
##     alpha = 0.6528545 
##     beta  = 0.1441675 
##     gamma = 0.09784922 
## 
##   Initial states:
##      l[0]       b[0]      s[0]    s[-1]    s[-2]     s[-3]
##  5.945592 0.07062881 0.9309236 1.177883 1.074851 0.8163427
## 
##   sigma^2:  0.0032
## 
##      AIC     AICc      BIC 
## 1680.929 1681.794 1711.389
# 3) Fit specific models for comparison
fit_models <- gas_data |>
  model(
    ETS_MAM = ETS(Gas ~ error("M") + trend("A") + season("M")),
    ETS_MAdM = ETS(Gas ~ error("M") + trend("Ad") + season("M")),
    ETS_AAM = ETS(Gas ~ error("A") + trend("A") + season("M"))
  )

report(fit_models)
## Warning in report.mdl_df(fit_models): Model reporting is only supported for
## individual models, so a glance will be shown. To see the report for a specific
## model, use `select()` and `filter()` to identify a single model.
## # A tibble: 3 × 9
##   .model     sigma2 log_lik   AIC  AICc   BIC   MSE  AMSE    MAE
##   <chr>       <dbl>   <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>  <dbl>
## 1 ETS_MAM   0.00324   -831. 1681. 1682. 1711.  21.1  32.2 0.0413
## 2 ETS_MAdM  0.00329   -832. 1684. 1685. 1718.  21.1  32.0 0.0417
## 3 ETS_AAM  18.2       -899. 1816. 1817. 1847.  17.6  25.1 2.84
# 4) Forecast next 5 years (20 quarters)
h <- 20

fc <- fit_models |>
  forecast(h = h)

# Plot forecasts
autoplot(fc, gas_data) +
  labs(title = "Gas Forecasts — ETS Models",
       y = "Gas Production",
       x = "Quarter")

# 5) Compare accuracy (training RMSE)
accuracy(fit_models)
## # A tibble: 3 × 10
##   .model   .type          ME  RMSE   MAE    MPE  MAPE  MASE RMSSE    ACF1
##   <chr>    <chr>       <dbl> <dbl> <dbl>  <dbl> <dbl> <dbl> <dbl>   <dbl>
## 1 ETS_MAM  Training -0.115    4.60  3.02  0.199  4.08 0.542 0.606 -0.0131
## 2 ETS_MAdM Training -0.00439  4.59  3.03  0.326  4.10 0.544 0.606 -0.0217
## 3 ETS_AAM  Training  0.218    4.19  2.84 -0.920  5.03 0.510 0.553  0.0405

8.8

# ---------------------------
# 0) Recreate "retail time series" from aus_retail
#    (Choose ONE series: set state + industry)
# ---------------------------

state_pick <- "New South Wales"
industry_pick <- "Cafes, restaurants and takeaway food services"

retail_data <- aus_retail |>
  filter(State == state_pick, Industry == industry_pick) |>
  select(Month, Turnover)

# Plot the series
autoplot(retail_data, Turnover) +
  labs(
    title = paste("Retail turnover:", state_pick, "|", industry_pick),
    y = "Turnover",
    x = "Month"
  )

# ---------------------------
# 1) Why multiplicative seasonality?
#    (Answer in words; code to help you see it)
# ---------------------------

# Seasonal plot: if seasonal amplitude grows with level -> multiplicative seasonality makes sense
retail_data |>
  gg_season(Turnover) +
  labs(title = "Seasonal plot of turnover")
## Warning: `gg_season()` was deprecated in feasts 0.4.2.
## ℹ Please use `ggtime::gg_season()` instead.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

# ---------------------------
# 2) Holt-Winters multiplicative method (ETS(M,A,M))
#    + experiment with damped trend ETS(M,Ad,M)
# ---------------------------

fit_hw <- retail_data |>
  model(
    HW_mult        = ETS(Turnover ~ error("M") + trend("A")  + season("M")),
    HW_mult_damped = ETS(Turnover ~ error("M") + trend("Ad") + season("M"))
  )

# Model details
report(fit_hw |> select(HW_mult))
## Series: Turnover 
## Model: ETS(M,A,M) 
##   Smoothing parameters:
##     alpha = 0.7729142 
##     beta  = 0.004491741 
##     gamma = 0.0001004778 
## 
##   Initial states:
##      l[0]      b[0]     s[0]     s[-1]    s[-2]    s[-3]    s[-4]   s[-5]
##  142.3738 0.4409842 0.994291 0.9301837 1.022985 1.140497 1.023531 1.01942
##     s[-6]     s[-7]     s[-8]    s[-9]    s[-10]    s[-11]
##  0.987387 0.9847622 0.9841327 0.944313 0.9872559 0.9812413
## 
##   sigma^2:  0.0015
## 
##      AIC     AICc      BIC 
## 5253.739 5255.186 5323.253
report(fit_hw |> select(HW_mult_damped))
## Series: Turnover 
## Model: ETS(M,Ad,M) 
##   Smoothing parameters:
##     alpha = 0.7681617 
##     beta  = 0.009440331 
##     gamma = 0.0001011594 
##     phi   = 0.9799997 
## 
##   Initial states:
##      l[0]     b[0]      s[0]     s[-1]   s[-2]    s[-3]    s[-4]    s[-5]
##  142.2247 1.499298 0.9910721 0.9273323 1.02064 1.138143 1.022578 1.022682
##      s[-6]     s[-7]     s[-8]     s[-9]    s[-10]    s[-11]
##  0.9925157 0.9886548 0.9869216 0.9442324 0.9864236 0.9788046
## 
##   sigma^2:  0.0015
## 
##      AIC     AICc      BIC 
## 5259.854 5261.475 5333.457
# ---------------------------
# ------------------------------------------------
# 3) Compare one-step training RMSE
# ------------------------------------------------
train_rmse <- accuracy(fit_hw) |>
  select(.model, RMSE)

train_rmse
## # A tibble: 2 × 2
##   .model          RMSE
##   <chr>          <dbl>
## 1 HW_mult         19.5
## 2 HW_mult_damped  19.4
# ------------------------------------------------
# 4) Set best_fit manually (CHOOSE ONE after seeing train_rmse)
#    Uncomment ONE of the two lines below
# ------------------------------------------------

# best_fit <- fit_hw |> select(HW_mult)
best_fit <- fit_hw |> select(HW_mult_damped)

# ------------------------------------------------
# 5) Residual diagnostics for best model
# ------------------------------------------------
gg_tsresiduals(best_fit)
## Warning: `gg_tsresiduals()` was deprecated in feasts 0.4.2.
## ℹ Please use `ggtime::gg_tsresiduals()` instead.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

augment(best_fit) |>
  features(.resid, ljung_box, lag = 24)
## # A tibble: 1 × 3
##   .model         lb_stat lb_pvalue
##   <chr>            <dbl>     <dbl>
## 1 HW_mult_damped    32.2     0.123
# ---------------------------
# 6) Test set RMSE: train through end of 2010, compare to seasonal naïve
# ---------------------------

train <- retail_data |>
  filter(Month <= yearmonth("2010 Dec"))

test <- retail_data |>
  filter(Month > yearmonth("2010 Dec"))

fit_train <- train |>
  model(
    HW_mult        = ETS(Turnover ~ error("M") + trend("A")  + season("M")),
    HW_mult_damped = ETS(Turnover ~ error("M") + trend("Ad") + season("M")),
    SNAIVE         = SNAIVE(Turnover)
  )

fc_test <- fit_train |>
  forecast(h = nrow(test))

test_rmse <- accuracy(fc_test, test) |>
  select(.model, RMSE, MAE, MAPE)

test_rmse
## # A tibble: 3 × 4
##   .model          RMSE   MAE  MAPE
##   <chr>          <dbl> <dbl> <dbl>
## 1 HW_mult         145.  121.  10.5
## 2 HW_mult_damped  240.  191.  16.0
## 3 SNAIVE          286.  232.  19.5
# Which beats seasonal naive?
test_rmse |>
  arrange(RMSE)
## # A tibble: 3 × 4
##   .model          RMSE   MAE  MAPE
##   <chr>          <dbl> <dbl> <dbl>
## 1 HW_mult         145.  121.  10.5
## 2 HW_mult_damped  240.  191.  16.0
## 3 SNAIVE          286.  232.  19.5

8.9

# ------------------------------------------------
# 0) Train/test split (same as before)
# ------------------------------------------------
train <- retail_data |>
  filter(Month <= yearmonth("2010 Dec"))

test <- retail_data |>
  filter(Month > yearmonth("2010 Dec"))

h_test <- nrow(test)

# ------------------------------------------------
# 1) Choose Box-Cox lambda (automatic via Guerrero)
# ------------------------------------------------
lambda <- train |>
  features(Turnover, features = guerrero) |>
  pull(lambda_guerrero)

lambda
## [1] 0.06183721
# ------------------------------------------------
# 2) STL decomposition on Box-Cox(Turnover), then ETS on seasonally adjusted
#    (ETS is fit to the season_adjust component)
# ------------------------------------------------
fit_stl_ets <- train |>
  model(
    STL_ETS = decomposition_model(
      STL(box_cox(Turnover, lambda) ~ season(window = "periodic")),
      ETS(season_adjust)
    )
  )

report(fit_stl_ets)
## Series: Turnover 
## Model: STL decomposition model 
## Transformation: box_cox(Turnover, lambda) 
## Combination: season_adjust + season_year
## 
## ========================================
## 
## Series: season_adjust 
## Model: ETS(A,A,N) 
##   Smoothing parameters:
##     alpha = 0.8077561 
##     beta  = 0.0001000008 
## 
##   Initial states:
##      l[0]        b[0]
##  5.865736 0.007885814
## 
##   sigma^2:  0.0033
## 
##      AIC     AICc      BIC 
## 55.17888 55.35587 74.39660 
## 
## Series: season_year 
## Model: SNAIVE 
## 
## sigma^2: 0
# Forecast on the test period
fc_stl_ets <- fit_stl_ets |>
  forecast(h = h_test)

# ------------------------------------------------
# 3) Fit your previous benchmarks on the SAME training data
# ------------------------------------------------
fit_bench <- train |>
  model(
    HW_mult = ETS(Turnover ~ error("M") + trend("A") + season("M")),
    SNAIVE  = SNAIVE(Turnover)
  )

fc_bench <- fit_bench |>
  forecast(h = h_test)

# ------------------------------------------------
# 4) Compare test-set RMSE across methods
# ------------------------------------------------
acc_compare <- bind_rows(
  accuracy(fc_stl_ets, test),
  accuracy(fc_bench, test)
) |>
  select(.model, RMSE, MAE, MAPE)

acc_compare
## # A tibble: 3 × 4
##   .model   RMSE   MAE  MAPE
##   <chr>   <dbl> <dbl> <dbl>
## 1 STL_ETS  82.6  72.6  7.12
## 2 HW_mult 145.  121.  10.5 
## 3 SNAIVE  286.  232.  19.5
# Optional: just RMSE sorted
acc_compare |>
  arrange(RMSE)
## # A tibble: 3 × 4
##   .model   RMSE   MAE  MAPE
##   <chr>   <dbl> <dbl> <dbl>
## 1 STL_ETS  82.6  72.6  7.12
## 2 HW_mult 145.  121.  10.5 
## 3 SNAIVE  286.  232.  19.5
# ------------------------------------------------
# 5) Visual comparison of forecasts on the test window
# ------------------------------------------------
bind_rows(
  fc_stl_ets |> mutate(.model = "STL_ETS"),
  fc_bench
) |>
  autoplot(retail_data) +
  autolayer(test, Turnover, alpha = 0.6) +
  labs(title = "Test-period forecasts: STL(BoxCox)+ETS vs HW_mult vs SNAIVE",
       y = "Turnover", x = "Month")

#Using an STL decomposition applied to a Box–Cox transformed turnover series, followed by ETS on the seasonally adjusted data, produced a test RMSE of ____. This compares to ____ for Holt-Winters multiplicative (ETS(M,A,M)) and ____ for seasonal naïve. Therefore, STL+ETS (does/does not) improve on the best previous method on the test set.