Objective

Assignment 5 involves answering questions 8.1, 8.5, 8.6, 8.7, 8.8, 8.9 from the textbook Forecasting: principles and practice by Rob J Hyndman and George Athanasopoulos.

Do exercises 8.1, 8.5, 8.6, 8.7, 8.8, 8.9 in Hyndman. Please submit both the link to your Rpubs and the .pdf file with your run code

Exercise 8.1

  1. Consider the the number of pigs slaughtered in Victoria, available in the aus_livestock dataset.

  1. Use the ETS() function to estimate the equivalent model for simple exponential smoothing. Find the optimal values of α and ℓ0, and generate forecasts for the next four months.

  2. Compute a 95% prediction interval for the first forecast using ^y ± 1.96s where s is the standard deviation of the residuals. Compare your interval with the interval produced by R.

Part A

data(aus_livestock)

pigs_victoria <- aus_livestock %>%
  filter(Animal == "Pigs", State == "Victoria")
fit <- pigs_victoria %>%
  model(ETS(Count ~ error("A") + trend("N") + season("N")))

report(fit)
## 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
forecasts <- fit %>%
  forecast(h = "4 months")
forecasts %>%
  autoplot(pigs_victoria) +
  geom_line(aes(y = .fitted), col="#D55E00",
            data = augment(fit)) +
  labs(y="% of GDP", title="Slaughtered Pigs in Victoria") +
  guides(colour = "none")

Part B

When looking at the 95% prediction interval we produced the lower bound is slightly higher, and the upper bound is slightly lower compared the 95% interval produced by R.

# Calculate the residuals from the fitted model
residuals <- fit %>%
  augment() %>%
  pull(.resid)

# Compute the standard deviation of the residuals
residuals_sd <- sd(residuals, na.rm = TRUE)

# Generate the first forecast
first_forecast <- forecasts$.mean[1]

# Compute the 95% prediction interval manually
lower_bound <- first_forecast - 1.96 * residuals_sd
upper_bound <- first_forecast + 1.96 * residuals_sd

print(lower_bound)
## [1] 76871.01
print(upper_bound)
## [1] 113502.1
# Compare with the prediction interval produced by R
forecasts %>%
  hilo(level = 95)
## # A tsibble: 4 x 7 [1M]
## # Key:       Animal, State, .model [1]
##   Animal State   .model    Month             Count  .mean                  `95%`
##   <fct>  <fct>   <chr>     <mth>            <dist>  <dbl>                 <hilo>
## 1 Pigs   Victor… "ETS(… 2019 Jan N(95187, 8.7e+07) 95187. [76854.79, 113518.3]95
## 2 Pigs   Victor… "ETS(… 2019 Feb N(95187, 9.7e+07) 95187. [75927.17, 114445.9]95
## 3 Pigs   Victor… "ETS(… 2019 Mar N(95187, 1.1e+08) 95187. [75042.22, 115330.9]95
## 4 Pigs   Victor… "ETS(… 2019 Apr N(95187, 1.1e+08) 95187. [74194.54, 116178.6]95

Exercise 8.5

Data set global_economy contains the annual Exports from many countries. Select one country to analyse.

  1. Plot the Exports series and discuss the main features of the data.

  2. Use an ETS(A,N,N) model to forecast the series, and plot the forecasts.

  3. Compute the RMSE values for the training data.

  4. Compare the results to those from an ETS(A,A,N) model. (Remember that the trended model is using one more parameter than the simpler model.) Discuss the merits of the two forecasting methods for this data set.

  5. Compare the forecasts from both methods. Which do you think is best?

  6. Calculate a 95% prediction interval for the first forecast for each model, using the RMSE values and assuming normal errors. Compare your intervals with those produced using R.

Part A

We can observe that there is an upward growth trend.

There are some random fluctuations throughout this period with a sharp decline year sometime between 2000 and 2010 - which quickly rebounded in the subsequent years.

data("global_economy")
exports_australia <- global_economy %>%
    filter(Country == "Australia")

autoplot(exports_australia, Exports)

Part B

fit_ets <- exports_australia %>%
  model(ETS(Exports ~ error("A") + trend("N") + season("N")))

report(fit_ets)
## Series: Exports 
## Model: ETS(A,N,N) 
##   Smoothing parameters:
##     alpha = 0.5659948 
## 
##   Initial states:
##      l[0]
##  12.98943
## 
##   sigma^2:  1.3621
## 
##      AIC     AICc      BIC 
## 257.3943 257.8387 263.5756
forecasts_2 <- fit_ets %>%
  forecast(h = "10 years")
forecasts_2 %>%
  autoplot(exports_australia) +
  geom_line(aes(y = .fitted), col="#D55E00",
            data = augment(fit_ets)) +
  labs(y="% of GDP", title="Slaughtered Pigs in Victoria") +
  guides(colour = "none")

Part C

accuracy_metrics <- fit_ets %>%
  accuracy()

print(accuracy_metrics)
## # A tibble: 1 × 11
##   Country   .model        .type    ME  RMSE   MAE   MPE  MAPE  MASE RMSSE   ACF1
##   <fct>     <chr>         <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>  <dbl>
## 1 Australia "ETS(Exports… Trai… 0.232  1.15 0.914  1.09  5.41 0.928 0.928 0.0125
# Print RMSE
accuracy_metrics %>%
  pull(RMSE)
## [1] 1.146794

Part D & E

Looking at metrics like the RMSE (Root Mean Squared Error), and MAE (Mean Absolute Error) the ETS (A,A,N) has scores lower in both metrics (1.116727 and 0.892642) compared to the ETS (A,N,N) model (1.146794 and 0.9135835). This tell us the ETS (A,A,N) model is more accurate.

Overall the ETS (A,A,N) is more accurate in it’s predictions. The ETS (A,N,N) captures the residuals better.

The ETS (A,N,N) model is effective when the data appears stationary and does not display a clear trend and seasonality. Looking at our data visually, we can observe a clear upward trend. This model was not as effective in modeling the data compared to the EST (A,A,N) as seen by the metrics.

The ETS (A,A,N) model is the best model for our data because it properly captures the trend seen in the dataset.

fit_ets_aan <- exports_australia %>%
  model(ETS(Exports ~ error("A") + trend("A") + season("N")))
accuracy_metrics_ets_ann <- fit_ets %>%
  accuracy()

accuracy_metrics_ets_aan <- fit_ets_aan %>%
  accuracy()

print(accuracy_metrics_ets_ann)
## # A tibble: 1 × 11
##   Country   .model        .type    ME  RMSE   MAE   MPE  MAPE  MASE RMSSE   ACF1
##   <fct>     <chr>         <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>  <dbl>
## 1 Australia "ETS(Exports… Trai… 0.232  1.15 0.914  1.09  5.41 0.928 0.928 0.0125
print(accuracy_metrics_ets_aan)
## # A tibble: 1 × 11
##   Country   .model     .type       ME  RMSE   MAE    MPE  MAPE  MASE RMSSE  ACF1
##   <fct>     <chr>      <chr>    <dbl> <dbl> <dbl>  <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Australia "ETS(Expo… Trai… -7.46e-7  1.12 0.893 -0.387  5.39 0.907 0.904 0.109
forecasts_aan <- fit_ets_aan %>%
  forecast(h = "10 years")

forecasts_aan %>%
  autoplot(exports_australia) +
   geom_line(aes(y = .fitted), col="#D55E00",
            data = augment(fit_ets_aan)) +
  labs(y="% of GDP", title="Slaughtered Pigs in Victoria") +
  guides(colour = "none")

Part F

The predictions for both are close. There are slight variations. For example, the EST (A,N,N) model has a slightly lower bound and higher upper bound compared to the intervals produced in R. For the EST (A,A,N) build, it’s the opposite: The lower bound is higher and the upper bound is lower compared to the intervals produced in R.

# Extract the RMSE values from each model
rmse_ets_ann <- 1.146794  # RMSE for ETS(A,N,N)
rmse_ets_aan <- 1.116727  # RMSE for ETS(A,A,N)

# Calculate the first forecasted value for both models
first_forecast_ann <- forecasts_2$.mean[1]  # ETS(A,N,N) first forecast
first_forecast_aan <- forecasts_aan$.mean[1]  # ETS(A,A,N) first forecast

# Calculate the 95% prediction interval for both models
lower_ann <- first_forecast_ann - 1.96 * rmse_ets_ann
upper_ann <- first_forecast_ann + 1.96 * rmse_ets_ann

lower_aan <- first_forecast_aan - 1.96 * rmse_ets_aan
upper_aan <- first_forecast_aan + 1.96 * rmse_ets_aan

cat("Manual 95% Prediction Interval (ETS(A,N,N)): [", lower_ann, ", ", upper_ann, "]\n")
## Manual 95% Prediction Interval (ETS(A,N,N)): [ 18.35944 ,  22.85488 ]
cat("Manual 95% Prediction Interval (ETS(A,A,N)): [", lower_aan, ", ", upper_aan, "]\n")
## Manual 95% Prediction Interval (ETS(A,A,N)): [ 18.64986 ,  23.02743 ]
# For ETS(A,N,N) model
r_prediction_interval_ann <- forecasts_2 %>%
  hilo(level = 95) %>%
  slice(1) %>%
  pull(`95%`)

cat("R's 95% Prediction Interval (ETS(A,N,N)): [", r_prediction_interval_ann$lower, ", ", r_prediction_interval_ann$upper, "]\n")
## R's 95% Prediction Interval (ETS(A,N,N)): [ 18.3197 ,  22.89462 ]
# For ETS(A,A,N) model
r_prediction_interval_aan <- forecasts_aan %>%
  hilo(level = 95) %>%
  slice(1) %>%
  pull(`95%`)

cat("R's 95% Prediction Interval (ETS(A,A,N)): [", r_prediction_interval_aan$lower, ", ", r_prediction_interval_aan$upper, "]\n")
## R's 95% Prediction Interval (ETS(A,A,N)): [ 18.57028 ,  23.107 ]

Question 8.6

Forecast the Chinese GDP from the global_economy data set using an ETS model. Experiment with the various options in the ETS() function to see how much the forecasts change with damped trend, or with a Box-Cox transformation. Try to develop an intuition of what each is doing to the forecasts.

Looking at the data from China from the global_economy dataset, we can see an upward trend that has a wide fluctuations in its peak periods. These increases aren’t of a consistent range. The damped trend best captures the data because it gradually tapers off the upward trend and most accurately represents the characteristics of real world data. An upward trend is going to taper off at some point and the model should take this into effect.

gdp_china <- global_economy %>%
  filter(Country == "China")

gdp_china %>% autoplot(GDP)

lambda <- gdp_china %>% features(GDP, features = guerrero) %>% pull(lambda_guerrero)

p1 <- gdp_china %>% autoplot(GDP) + labs(title = "Original scale", y = element_blank()) 


p2 <- gdp_china %>% autoplot(box_cox(GDP, lambda)) + labs(title = "Box-Cox Transformation", y = element_blank())

p1 + p2

fit <- gdp_china %>%
  model(
    Auto = ETS(GDP),
    AAN = ETS(GDP ~ error("A") + trend("A") + season("N")),
    AAdN = ETS(GDP ~ error("A") + trend("Ad") + season("N")),
    AMN = ETS(GDP ~ error("A") + trend("M") + season("N")),
    Box_Cox = ETS(box_cox(GDP, lambda))
  )

fc <- fit %>% forecast(h = 20)
fc %>%
  autoplot(gdp_china) +
  facet_wrap(vars(.model), scales = "free_y")

accuracy(fit)
## # A tibble: 5 × 11
##   Country .model  .type              ME    RMSE     MAE    MPE  MAPE  MASE RMSSE
##   <fct>   <chr>   <chr>           <dbl>   <dbl>   <dbl>  <dbl> <dbl> <dbl> <dbl>
## 1 China   Auto    Training      4.13e10 2.00e11 9.66e10  2.11   7.72 0.446 0.477
## 2 China   AAN     Training      2.36e10 1.90e11 9.59e10  1.41   7.62 0.442 0.453
## 3 China   AAdN    Training      2.95e10 1.90e11 9.49e10  1.62   7.62 0.438 0.454
## 4 China   AMN     Training     -2.79e10 2.23e11 1.07e11 -0.734  8.20 0.496 0.533
## 5 China   Box_Cox Training     -2.75e10 2.88e11 1.25e11  0.607  7.17 0.578 0.688
## # ℹ 1 more variable: ACF1 <dbl>

Exercise 8.7

Find an ETS model for the Gas data from aus_production and forecast the next few years. Why is multiplicative seasonality necessary here? Experiment with making the trend damped. Does it improve the forecasts?

Multiplicative seasonality is necessary here because the seasonal flunctuations throughout the trend growth gradually increase. If the flunctuations were consistent patterns, then an additive seasonality would be appropriate - but this is not the case.

data(aus_production)

aus_gas <- aus_production%>%
  select(Quarter,Gas)

autoplot(aus_gas,Gas)

Making the trend damped gives us the best predictive model because it give us the most realistic forecast. Dampening the trend, basically slows down the long term growth over time. This is important especially when we don’t expect the growth to continue growing at the same rate.

Looking at the metrics. The Additive Error, Damped Additive Trend, and Multiplicative Seasonality (AAdM) performed the best. It has the lowest Root Mean Squared Error (4.22), Mean Absolute Error (2.81) - which means it minimizes both error magnitude and forecast errors. It’s Mean Absolute Percentage Error is also the lowest at 4.11% and therefore the most accurate.

fit2 <- aus_gas %>%
  model(
    AAA = ETS(Gas ~ error("A") + trend("A") + season("A")),
    AAdA = ETS(Gas ~ error("A") + trend("Ad") + season("A")),
    AAdM = ETS(Gas ~ error("A") + trend("Ad") + season("M")),
    AMA = ETS(Gas ~ error("A") + trend("M") + season("A")),
    AMM = ETS(Gas ~ error("A") + trend("M") + season("M"))
  )

fc2 <- fit2 %>% 
  forecast(h = 20)
fc2 %>%
  autoplot(aus_gas) +
  facet_wrap(vars(.model), scales = "free_y")

accuracy(fit2)
## # A tibble: 5 × 10
##   .model .type          ME  RMSE   MAE     MPE  MAPE  MASE RMSSE   ACF1
##   <chr>  <chr>       <dbl> <dbl> <dbl>   <dbl> <dbl> <dbl> <dbl>  <dbl>
## 1 AAA    Training  0.00525  4.76  3.35 -4.69   10.9  0.600 0.628 0.0772
## 2 AAdA   Training  0.600    4.58  3.16  0.460   7.39 0.566 0.604 0.0660
## 3 AAdM   Training  0.548    4.22  2.81  1.32    4.11 0.505 0.556 0.0265
## 4 AMA    Training -0.444    4.67  3.30 -0.0632  7.85 0.592 0.616 0.0719
## 5 AMM    Training  0.465    4.23  2.84  1.26    4.12 0.509 0.558 0.0321

Exercise 8.8

Recall your retail time series data (from Exercise 7 in Section 2.10).

  1. Why is multiplicative seasonality necessary for this series?

  2. Apply Holt-Winters’ multiplicative method to the data. Experiment with making the trend damped.

  3. Compare the RMSE of the one-step forecasts from the two methods. Which do you prefer?

  4. Check that the residuals from the best method look like white noise.

  5. Now find the test set RMSE, while training the model to the end of 2010. Can you beat the seasonal naïve approach from Exercise 7 in Section 5.11?

Part A

Multiplicative seasonality is necessary for the series because of the fluctuation patterns. There is are seasonal peaks that grow larger throughout time. Capturing this trend through multiplicative seasonality makes sense.

set.seed(12345678)

myseries <- aus_retail |>
  filter(`Series ID` == sample(aus_retail$`Series ID`,1))
autoplot(myseries)
## Plot variable not specified, automatically selected `.vars = Turnover`

Part B

fit3 <- myseries %>%
  model(
    MAM = ETS(Turnover ~ error("M") + trend("A") + season("M")),
    MAdM = ETS(Turnover ~ error("M") + trend("Ad") + season("M"))
  )

fc3 <- fit3 %>% 
  forecast(h = 60)
fc3 %>%
  autoplot(myseries) +
  facet_wrap(vars(.model), scales = "free_y")

accuracy(fit3)
## # A tibble: 2 × 12
##   State      Industry .model .type      ME  RMSE   MAE     MPE  MAPE  MASE RMSSE
##   <chr>      <chr>    <chr>  <chr>   <dbl> <dbl> <dbl>   <dbl> <dbl> <dbl> <dbl>
## 1 Northern … Clothin… MAM    Trai… -0.0128 0.613 0.450 -0.469   5.15 0.513 0.529
## 2 Northern … Clothin… MAdM   Trai…  0.0398 0.616 0.444 -0.0723  5.06 0.507 0.531
## # ℹ 1 more variable: ACF1 <dbl>

Part C

The MAM model has a lower RMSE value at 0.6130 compared to the MAdM model (0.6155). This tells us that the MAM model is slightly better at minimizing large errors in the forecast. However, despite this, the MAdM model has better overall accuracy based on the MAPE metric (5.06). For this reason, I prefer the MAdM model.

Part D

The residuals look consistent, and the fluctuations appear to be random. The histogram of the residuals also appears close to symmetrical. We can safely conclude that the residuals behave like white noise.

fit3 %>%
  select(MAdM) %>%
  gg_tsresiduals()

Part E

Now find the test set RMSE, while training the model to the end of 2010. Can you beat the seasonal naïve approach from Exercise 7 in Section 5.11?

The MAdM model from part D outperforms the seasonal naive approach as it has the lowest RMSE, MAE and MAPE - making it the most accurate approach.

myseries_train <- myseries |>
  filter(year(Month) < 2011)
myseries_test <- myseries %>% anti_join(myseries_train)
## Joining with `by = join_by(State, Industry, `Series ID`, Month, Turnover)`
fits <- myseries_train %>%
  model(
  SNaive = SNAIVE(Turnover),
  MAM = ETS(Turnover ~ error("M") + trend("A") + season("M")),
  MAdM = ETS(Turnover ~ error("M") + trend("Ad") + season("M"))      
)
fc <- fits |>
  forecast(new_data = anti_join(myseries, myseries_train))
## Joining with `by = join_by(State, Industry, `Series ID`, Month, Turnover)`
fc %>% autoplot(myseries)

fc <- fits %>% 
  forecast(new_data = myseries_test)

fc %>% accuracy(myseries)
## # A tibble: 3 × 12
##   .model State Industry .type     ME  RMSE   MAE     MPE  MAPE  MASE RMSSE  ACF1
##   <chr>  <chr> <chr>    <chr>  <dbl> <dbl> <dbl>   <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 MAM    Nort… Clothin… Test  -1.54   1.78 1.60  -12.3   12.6  1.75  1.47  0.495
## 2 MAdM   Nort… Clothin… Test   0.163  1.15 0.878   0.194  6.45 0.960 0.949 0.501
## 3 SNaive Nort… Clothin… Test   0.836  1.55 1.24    5.94   9.06 1.36  1.28  0.601

Exercise 8.9

For the same retail data, try an STL decomposition applied to the Box-Cox transformed series, followed by ETS on the seasonally adjusted data. How does that compare with your best previous forecasts on the test set?

There are larger residuals in the ETS MAdM model than in the Box-Cox model. The Box-Cox transformation has a smoother trend and smaller residuals.

# box cox transformation
lambda <- myseries_train %>% features(Turnover, features = guerrero) %>% pull(lambda_guerrero)

myseries_train %>%
  model(STL(box_cox(Turnover, lambda = lambda))) %>%
  components() %>%
  autoplot()

# ETS model 
mAdm_model <- fits %>%
  select(MAdM)

mAdm_model %>%  
  components() %>%
  autoplot()
## Warning: Removed 12 rows containing missing values or values outside the scale range
## (`geom_line()`).