install.packages("fpp3")
## Installing package into '/cloud/lib/x86_64-pc-linux-gnu-library/4.6'
## (as 'lib' is unspecified)
library(fpp3)
## ── Attaching packages ──────────────────────────────────────────── fpp3 1.0.3 ──
## ✔ tibble      3.3.1     ✔ tsibble     1.2.0
## ✔ dplyr       1.2.1     ✔ tsibbledata 0.4.1
## ✔ tidyr       1.3.2     ✔ ggtime      1.0.0
## ✔ lubridate   1.9.5     ✔ feasts      0.5.0
## ✔ ggplot2     4.0.3     ✔ fable       0.5.0
## ── 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()

Exercise 8.1- Pigs Slaughtered in Victoria

8.1a- Use ETS() to estimate a simple exponential smoothing model for the number of pigs slaughtered in Victoria. Find the optimal values of alpha and initial level, and forecast the next four months.

pigs <- aus_livestock |>
  filter(State == "Victoria", Animal == "Pigs")

fit_pigs <- pigs |>
  model(SES = ETS(Count ~ error("A") + trend("N") + season("N")))

tidy(fit_pigs)
fc_pigs <- fit_pigs |>
  forecast(h = 4)

fc_pigs
fc_pigs |>
  autoplot(pigs)

The estimated alpha is 0.3221, and the initial level is 100,646.6. The forecast for the next four months is about 95,187 pigs per month. Since simple exponential smoothing does not account for trend or seasonality, the forecasts remain constant.

8.1b- Calculate a 95% prediction interval for the first forecast using the forecast ± 1.96 times the residual standard deviation. Compare it with R’s interval.

s <- sd(augment(fit_pigs)$.innov, na.rm = TRUE)

first_forecast <- fc_pigs$.mean[1]

lower <- first_forecast - 1.96 * s
upper <- first_forecast + 1.96 * s

c(lower = lower, upper = upper)
##     lower     upper 
##  76871.01 113502.10
hilo(fc_pigs, 95)

The manually calculated 95% prediction interval is 76,871 to 113,502 pigs. R gives a similar interval of 76,855 to 113,518 pigs. The results are very close, with a small difference due to how the standard deviation is calculated.

Exercise 8.5- Exports (Algeria)

8.5a- Plot the exports series and describe the main features.

algeria <- global_economy |>
  filter(Country == "Algeria")

algeria |>
  autoplot(Exports) +
  labs(title = "Exports: Algeria", y = "% of GDP")

Algeria’s exports fluctuate a lot over time, with several increases and decreases. Exports reach their highest point around 2006–2008 before declining toward the end of the series. There is no clear overall trend.

8.5b- Use ETS(A,N,N) and plot the forecasts

fit_ann <- algeria |>
  model(ANN = ETS(Exports ~ error("A") +
                   trend("N") + season("N")))

fc_ann <- fit_ann |>
  forecast(h = 5)

fc_ann |>
  autoplot(algeria)

The ETS(A,N,N) model produces a flat forecast of about 22.44% of GDP for the next five years. Since the model does not include a trend or seasonality, the forecast stays the same while the prediction intervals become wider over time.

8.5c- Calculate training RMSE

accuracy(fit_ann) |>
  select(.model, RMSE)

The RMSE for the ETS(A,N,N) model is 5.87. This measures how far the fitted values are from the actual export values, on average, using the RMSE calculation.

8.5d- Compare ETS(A,N,N) with ETS(A,A,N)

fit_both <- algeria |>
  model(
    ANN = ETS(Exports ~ error("A") +
                trend("N") + season("N")),
    AAN = ETS(Exports ~ error("A") +
                trend("A") + season("N"))
  )

accuracy(fit_both) |>
  select(.model, RMSE)
glance(fit_both) |>
  select(.model, AICc)

The ETS(A,N,N) model has a lower AICc (447.16) than ETS(A,A,N) (452.33). This suggests the model without a trend is the better choice based on AICc.

8.5e- Compare the forecasts and choose the better model

fc_both <- fit_both |>
  forecast(h = 5)

fc_both |>
  autoplot(algeria)

The ETS(A,N,N) model produces a flat forecast, while ETS(A,A,N) predicts a slight decline. I would choose ETS(A,N,N) because it has the lower AICc and does not require a trend.

8.5f- Calculate and compare 95% prediction intervals

first_forecasts <- fc_both |>
  as_tibble() |>
  group_by(.model) |>
  slice_head(n = 1) |>
  ungroup() |>
  select(.model, .mean)

model_rmse <- accuracy(fit_both) |>
  select(.model, RMSE)

manual_intervals <- first_forecasts |>
  left_join(model_rmse, by = ".model") |>
  mutate(
    lower = .mean - 1.96 * RMSE,
    upper = .mean + 1.96 * RMSE
  ) |>
  select(.model, lower, upper)

manual_intervals
fc_both |>
  hilo(95) |>
  as_tibble() |>
  group_by(.model) |>
  slice_head(n = 1) |>
  ungroup() |>
  select(.model, Year, `95%`)

I calculated the 95% prediction intervals for both models using the forecast plus or minus 1.96 times the RMSE. R gives an interval of 10.75 to 34.14 for ETS(A,N,N) and 10.02 to 33.94 for ETS(A,A,N). Both intervals are similar, with small differences between the manual calculations and R’s results.

Exercise 8.6- China’s GDP

china <- global_economy |>
  filter(Country == "China")

china |>
  autoplot(GDP) +
  labs(title = "China GDP", y = "GDP")

fit_china <- china |>
  model(
    ETS = ETS(GDP),
    Damped = ETS(GDP ~ error("A") +
                   trend("Ad") + season("N")),
    BoxCox = ETS(box_cox(GDP, 0.3) ~
                   error("A") + trend("A") + season("N"))
  )

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

fc_china |>
  autoplot(china, level = NULL) +
  labs(title = "China GDP Forecasts", y = "GDP")

accuracy(fit_china) |>
  select(.model, RMSE)
glance(fit_china) |>
  select(.model, AICc)

China’s GDP increases rapidly over time. The regular ETS model continues this growth in its forecasts, while the damped model slows the trend over time. The Box-Cox transformation produces a different growth pattern. The forecasts show that the choice of model can make a large difference, especially when predicting 20 years into the future.

Exercise 8.7- Australian gas production

gas <- aus_production

gas |>
  autoplot(Gas) +
  labs(title = "Australian Gas Production",
       y = "Gas Production")

fit_gas <- gas |>
  model(
    ETS = ETS(Gas),
    Regular = ETS(Gas ~ error("M") +
                    trend("A") + season("M")),
    Damped = ETS(Gas ~ error("M") +
                   trend("Ad") + season("M"))
  )

fc_gas <- fit_gas |>
  forecast(h = "5 years")

fc_gas |>
  autoplot(gas, level = NULL) +
  labs(title = "Australian Gas Forecasts",
       y = "Gas Production")

accuracy(fit_gas) |>
  select(.model, RMSE)
glance(fit_gas) |>
  select(.model, AICc)
report(fit_gas)
## Warning in report.mdl_df(fit_gas): 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.

Australian gas production increases over time, and the seasonal fluctuations become larger as production increases. This makes multiplicative seasonality appropriate. The regular ETS model has an AICc of 1681.79, compared to 1685.02 for the damped model. Although the damped model has a slightly lower MSE, the regular model has the better AICc. I would choose the regular ETS model because adding a damped trend does not provide a clear improvement.

Exercise 8.8- Australian Retail Data

8.8a- Why is multiplicative seasonality necessary?

set.seed(12345678)

retail <- aus_retail |>
  filter(State == "Northern Territory")

retail <- retail |>
  filter(Industry == sample(unique(Industry), 1))

retail |>
  autoplot(Turnover)

Retail turnover generally increases over time, with larger fluctuations in later years. This suggests multiplicative seasonality may be appropriate because the size of the seasonal changes increases as turnover grows.

8.8b — Apply Holt-Winters’ multiplicative method and a damped trend

fit_retail <- retail |>
  model(
    Regular = ETS(Turnover ~ error("M") +
                    trend("A") + season("M")),
    Damped = ETS(Turnover ~ error("M") +
                   trend("Ad") + season("M"))
  )

fit_retail |>
  forecast(h = 24) |>
  autoplot(retail, level = NULL)

I compared regular and damped Holt-Winters models using multiplicative seasonality. Both models produce similar forecasts, but the damped model reduces the effect of the trend over time.

8.8c — Compare training RMSE

accuracy(fit_retail) |>
  select(.model, RMSE)

The damped model has a slightly lower RMSE (0.882) than the regular model (0.896). This means the damped model fits the training data slightly better.

8.8d — Check the residuals of the better model

best_model <- accuracy(fit_retail) |>
  arrange(RMSE) |>
  pull(.model) |>
  first()

fit_retail |>
  select(all_of(best_model)) |>
  gg_tsresiduals()

fit_retail |>
  select(all_of(best_model)) |>
  augment() |>
  features(.innov, ljung_box, lag = 24, dof = 0)

Most residuals are close to zero, with a few larger errors. The Ljung-Box test gives a p-value of 0.894, which suggests there is no significant remaining autocorrelation. This indicates that the damped model captures the main patterns in the data.

8.8e — Compare test-set RMSE with seasonal naïve

train <- retail |>
  filter(year(Month) <= 2010)

test_fit <- train |>
  model(
    Regular = ETS(Turnover ~ error("M") +
                    trend("A") + season("M")),
    Damped = ETS(Turnover ~ error("M") +
                   trend("Ad") + season("M")),
    SNAIVE = SNAIVE(Turnover)
  )

test_fc <- test_fit |>
  forecast(h = nrow(retail) - nrow(train))

test_fc |>
  accuracy(retail) |>
  select(.model, RMSE)
test_fc |>
  autoplot(retail, level = NULL)

The damped model had the lowest test RMSE at 3.56, followed by the regular model at 3.57 and seasonal naïve at 4.19. This shows that the damped model performed slightly better than the regular model, while both Holt-Winters models performed better than seasonal naïve.

Exercise 8.9- STL decomposition and forecasting

fit_stl <- train |>
  model(
    STL_ETS = decomposition_model(
      STL(box_cox(Turnover, 0.3) ~ season(window = "periodic")),
      ETS(season_adjust ~ season("N"))
    )
  )
fc_stl <- fit_stl |>
  forecast(h = nrow(retail) - nrow(train))

fc_stl |>
  autoplot(retail, level = NULL) +
  labs(title = "STL + ETS Forecast",
       y = "Turnover")

stl_rmse <- fc_stl |>
  accuracy(retail) |>
  select(.model, RMSE)

stl_rmse
test_fc |>
  accuracy(retail) |>
  select(.model, RMSE)
stl_rmse

I used STL decomposition with a Box-Cox transformation and ETS to forecast retail turnover. The STL + ETS model had a test RMSE of 3.70, which was higher than the damped Holt-Winters model at 3.56. This shows that the damped Holt-Winters model performed slightly better on the test data, although both models performed better than seasonal naïve.