LFMG Lab 5

Author

Luis Munoz Grass

1 Exponential Smoothing

2 Exercises

2.1 8.1

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

a. 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.

First, we load required libraries.

library(fpp3) # Forecasting package
library(ggplot2) # For visualizations
library(kableExtra) # Nicely formatted tables
library(feasts)

Load and prepare data

pig_data <- aus_livestock %>%
  filter(State == "Victoria", Animal == "Pigs") %>%
  select(Month, Count) %>%
  as_tsibble(index = Month)

autoplot(pig_data, Count) +
  ggtitle("Pigs Slaughtered in Victoria") +
  ylab("Count") +
  xlab("Year")

You can add options to executable code like this

pig_fit <- pig_data %>%
  model(ETS(Count ~ error("A") + trend("N") + season("N")))

# save model parameters
pig_fit_params <- tidy(pig_fit)
pig_fit_params %>%
  kable(caption = "ETS Parameters for Simple Exponential Smoothing") %>%
  kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover"))

ETS Parameters for Simple Exponential Smoothing
.model	term	estimate
ETS(Count ~ error("A") + trend("N") + season("N"))	alpha	3.221247e-01
ETS(Count ~ error("A") + trend("N") + season("N"))	l[0]	1.006466e+05

We see that the alpha is around 0.322, meaning 32.2% of the new observation is used in updating forecasts, while 67.8% is based on past values. The ℓ0 was around 100,647, so that will be the starting point for the smoothing process.

Now that we have optimal parameters we can forecast the next 4 months.

forecast_pig_slag <- pig_fit %>%
  forecast(h = "4 months", level = c(95))

autoplot(pig_data, Count) +
  autolayer(forecast_pig_slag, level = 95, alpha = 0.5) +
  ggtitle("Forecast: Pig Slaughtering in Victoria - 4 months") +
  ylab("Count") +
  xlab("Year")

b. 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.

First we need to calculate the standard deviation of the residuals from the fitted model

residuals_pig <- pig_fit %>% augment() %>% pull(.resid)

s_pig <- sd(residuals_pig, na.rm = TRUE)

s_pig

[1] 9344.666

Now we can extract the first value from the forecast, and compute the manual prediction interval

y_hat_1 <- forecast_pig_slag$.mean[1] 

lower_bound <- y_hat_1 - (1.96 * s_pig)
upper_bound <- y_hat_1 + (1.96 * s_pig)


c("Lower Bound" = lower_bound, "Forecast" = y_hat_1, "Upper Bound" = upper_bound)

Lower Bound    Forecast Upper Bound 
   76871.01    95186.56   113502.10

Let’s compare that now:

intervals <- forecast_pig_slag %>%
  hilo(95) %>% 
  as_tibble()  # Convert to tibble for easier handling

# Separate lower and upper bounds
intervals <- intervals %>%
  separate(`95%`, into = c("Lower", "Upper"), sep = ", ", convert = TRUE)

# Remove brackets and convert to numeric
r_lower <- as.numeric(gsub("\\[", "", intervals$Lower[1]))  # Remove "[" from lower bound
r_upper <- as.numeric(gsub("\\]", "", intervals$Upper[1]))  # Remove "]" from upper bound

# Print values
r_lower

[1] 76854.79

r_upper

[1] 113518.3

comparison <- tibble(
  Method = c("Manual Calculation", "R Forecast"),
  Lower = c(lower_bound, r_lower),
  Upper = c(upper_bound, r_upper)
)

# Display as a table
library(kableExtra)
comparison %>%
  kable(caption = "Comparison of Manual vs. R's Prediction Interval") %>%
  kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover"))

Comparison of Manual vs. R's Prediction Interval
Method	Lower	Upper
Manual Calculation	76871.01	113502.1
R Forecast	76854.79	113518.3

We see that we got almost identical results, meaning that we calculated the standard deviation correctly, and we can assume that the residuals are normally distributed. Differences may be due to rounding differences in standar deviation calculations.

2.2 8.5

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

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

For this exercise I want to focus on exports from Colombia. We’ll start by loading the data, plotting it and getting main descriptive statistics.

colombia_exports <- global_economy %>%
  filter(Country == "Colombia") %>%
  select(Year, Exports) %>%
  as_tsibble(index = Year)

autoplot(colombia_exports, Exports) +
  ggtitle("Annual Exports of Colombia") +
  ylab("Exports (USD Billions)") +
  xlab("Year") +
  theme_minimal()

We can see that the series does not show a clear upward or downward trend over the entire period, and instead exports fluctuate significantly over time. The exports increase and decrease in cycles, with sharp changes in certain periods, and the most significant drops occur around the early 1980s and mid-2010s.

On the other hand, the peaks around 1990 and early 2000s suggest strong export performance in those years, possibly due to high oil prices or trade agreements.

Finally, since the data is annual, there are no visible seasonal patterns.

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

First we will output the estimated smoothing parameter α and the initial level \(ℓ0\)

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

report(colombia_ets)

Series: Exports 
Model: ETS(A,N,N) 
  Smoothing parameters:
    alpha = 0.9707053 

  Initial states:
     l[0]
 15.61444

  sigma^2:  2.479

     AIC     AICc      BIC 
292.1271 292.5716 298.3084

With those values we can get the forecast and plot it

colombia_forecast <- colombia_ets %>%
  forecast(h = "10 years", level = 95)

autoplot(colombia_forecast, colombia_exports) +
  ggtitle("ETS Forecast: Colombia's Exports") +
  ylab("Exports (USD Billions)") +
  xlab("Year") +
  theme_minimal()

c. Compute the RMSE values for the training data.

accuracy_results <- colombia_ets %>%
  accuracy()

rmse_value <- accuracy_results %>%
  filter(.model == "ETS(Exports ~ error(\"A\") + trend(\"N\") + season(\"N\"))") %>%
  pull(RMSE)

rmse_value

[1] 1.547114

d. 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.

First, we fit an ETS (AAN) model

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

report(colombia_ets_2)

Series: Exports 
Model: ETS(A,A,N) 
  Smoothing parameters:
    alpha = 0.9283752 
    beta  = 0.0001000128 

  Initial states:
     l[0]       b[0]
 14.36393 0.01625495

  sigma^2:  2.5994

     AIC     AICc      BIC 
296.7668 297.9207 307.0690

Generate forecast and plot,

colombia_forecast <- colombia_ets %>%
  forecast(h = "10 years", level = 95) %>%
  as_tibble() %>%  # Convert to tibble first
  select(Year, Exports, .mean) %>%  # Ensure Year exists
  mutate(Model = "ETS(A,N,N)")  # Add model label

colombia_forecast_2 <- colombia_ets_2 %>%
  forecast(h = "10 years", level = 95) %>%
  as_tibble() %>%  # Convert to tibble first
  select(Year, Exports, .mean) %>%  # Ensure Year exists
  mutate(Model = "ETS(A,A,N)")  # Add model label

combined_forecasts <- bind_rows(colombia_forecast, colombia_forecast_2) %>%
  as_tsibble(index = Year, key = Model)  # Convert back to tsibble

autoplot(colombia_exports, Exports) +
  geom_line(data = combined_forecasts, aes(x = Year, y = .mean, color = Model), linewidth = 1) +
  ggtitle("ETS Forecast Comparison: With & Without Trend") +
  ylab("Exports (USD Billions)") +
  xlab("Year") +
  theme_minimal()

And get RMS so we can later compare both models:

accuracy_2 <- colombia_ets_2 %>%
  accuracy()

rmse_2 <- accuracy_2 %>%
  pull(RMSE)

comparison_rmse <- tibble(
  Model = c("ETS(A,N,N)", "ETS(A,A,N)"),
  RMSE = c(rmse_value, rmse_2)
)

comparison_rmse %>%
  kable(caption = "RMSE Comparison: ETS(A,N,N) vs ETS(A,A,N)") %>%
  kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover"))

RMSE Comparison: ETS(A,N,N) vs ETS(A,A,N)
Model	RMSE
ETS(A,N,N)	1.547114
ETS(A,A,N)	1.555670

ETS(A,N,N) has a slightly lower RMSE than ETS(A,A,N), meaning it fits the training data marginally better. Adding a trend (ETS(A,A,N)) did not significantly improve accuracy. Since ETS(A,A,N) uses one more parameter, this extra complexity does not justify the small RMSE difference.

Given that export levels fluctuate but do not show a strong, consistent trend, ETS(A,N,N) (Simple Exponential Smoothing) is the preferred model.

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

The forecast for ETS (ANN) is flat (constant) because this model assumes no trend. The prediction remains at a constant level with increasing uncertainty over time. At the same time, the ETS (AAN) forecast shows a very slight upward trend as it may have detected a weak trend in the data, however, since the trend is minimal, it may not significantly improve accuracy over time.

Since the exports data for Colombia does not have a strong, consistent trend, the ETS (ANN) seems like a better choice for this dataset.

f. 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.

Similar to 8.1 we start by extracting first forecasted value for each model.

col_y_hat_1_ets <- colombia_forecast$.mean[1]   # first forecast from ETS(A,N,N)
col_y_hat_1_ets_2 <- colombia_forecast_2$.mean[1]  # first forecast from ETS(A,A,N)

Since we already had RMSE values computed we can plug them to get upper and lower bounds for each model.

col_lower_ets <- col_y_hat_1_ets - (1.96 * rmse_value)
col_upper_ets <- col_y_hat_1_ets + (1.96 * rmse_value)

col_lower_ets_2 <- col_y_hat_1_ets_2 - (1.96 * rmse_2)
col_upper_ets_2 <- col_y_hat_1_ets_2 + (1.96 * rmse_2)

manual_intervals <- tibble(
  Model = c("ETS(A,N,N)", "ETS(A,A,N)"),
  Lower = c(col_lower_ets, col_lower_ets_2),
  Forecast = c(col_y_hat_1_ets, col_y_hat_1_ets_2),
  Upper = c(col_upper_ets, col_upper_ets_2)
)

manual_intervals %>%
  kable(caption = "Manual 95% Prediction Intervals") %>%
  kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover"))

Manual 95% Prediction Intervals
Model	Lower	Forecast	Upper
ETS(A,N,N)	11.52101	14.55335	17.58570
ETS(A,A,N)	11.52638	14.57549	17.62461

Now, we get the values straight from R to later compare.

# Extract 95% prediction intervals for ETS(A,N,N)
colombia_forecast <- colombia_ets %>%
  forecast(h = "10 years", level = 95)

intervals <- colombia_forecast %>%
  hilo(95) %>%
  as_tibble()  # Convert only AFTER extracting intervals

# Separate lower and upper bounds
intervals <- intervals %>%
  separate(`95%`, into = c("Lower", "Upper"), sep = ", ", convert = TRUE)

# Remove brackets and convert to numeric
col_r_lower_ets <- as.numeric(gsub("\\[", "", intervals$Lower[1]))  
col_r_upper_ets <- as.numeric(gsub("\\]", "", intervals$Upper[1]))  

# Extract 95% prediction intervals for ETS(A,N,N)
colombia_forecast_2 <- colombia_ets_2 %>%
  forecast(h = "10 years", level = 95)

intervals_2 <- colombia_forecast_2 %>%
  hilo(95) %>%
  as_tibble()  # Convert only AFTER extracting intervals

# Separate lower and upper bounds
intervals_2 <- intervals_2 %>%
  separate(`95%`, into = c("Lower", "Upper"), sep = ", ", convert = TRUE)

# Remove brackets and convert to numeric
col_r_lower_ets_2 <- as.numeric(gsub("\\[", "", intervals_2$Lower[1]))  
col_r_upper_ets_2 <- as.numeric(gsub("\\]", "", intervals_2$Upper[1]))

Putting R values in a table:

# Create comparison table
col_r_intervals <- tibble(
  Model = c("ETS(A,N,N)", "ETS(A,A,N)"),
  Lower = c(col_r_lower_ets, col_r_lower_ets_2),
  Forecast = c(col_y_hat_1_ets, col_y_hat_1_ets_2),
  Upper = c(col_r_upper_ets, col_r_upper_ets_2)
)

# Display table
col_r_intervals %>%
  kable(caption = "R Forecast 95% Prediction Intervals") %>%
  kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover"))

R Forecast 95% Prediction Intervals
Model	Lower	Forecast	Upper
ETS(A,N,N)	11.46739	14.55335	17.63931
ETS(A,A,N)	11.41553	14.57549	17.73546

So we see in the final comparison:

The forecasted values (.mean) are identical since both methods use the same ETS models.
The lower and upper bounds are slightly different, likely due to small rounding differences and how R calculates variance in forecast distributions.

2.3 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.

[Hint: use a relatively large value of h when forecasting, so you can clearly see the differences between the various options when plotting the forecasts.]

Let’s load and prepare the data

china_gdp <- global_economy %>%
  filter(Country == "China") %>%
  select(Year, GDP) %>%
  as_tsibble(index = Year)

autoplot(china_gdp, GDP) +
  ggtitle("Chinese GDP Over Time") +
  ylab("GDP (USD Billions)") +
  xlab("Year") +
  theme_minimal()

Now let’s start fitting different models:

china_ets_ANN <- china_gdp %>% model(ETS(GDP ~ error("A") + trend("N") + season("N")))
china_ets_AAN <- china_gdp %>% model(ETS(GDP ~ error("A") + trend("A") + season("N")))
china_ets_AAdN <- china_gdp %>% model(ETS(GDP ~ error("A") + trend("Ad") + season("N")))  # Damped trend
china_ets_boxcox <- china_gdp %>% model(ETS(box_cox(GDP, 0.3) ~ error("A") + trend("A") + season("N")))  # Box-Cox transformation

Generate some forecasts for each for the next 30 years:

# Forecasting GDP for 30 years
china_forecast_ANN <- china_ets_ANN %>% forecast(h = "30 years")
china_forecast_AAN <- china_ets_AAN %>% forecast(h = "30 years")
china_forecast_AAdN <- china_ets_AAdN %>% forecast(h = "30 years")
china_forecast_boxcox <- china_ets_boxcox %>% forecast(h = "30 years")

And plotting:

china_forecast_ANN <- china_forecast_ANN %>% mutate(Model = "ETS(A,N,N)")
china_forecast_AAN <- china_forecast_AAN %>% mutate(Model = "ETS(A,A,N)")
china_forecast_AAdN <- china_forecast_AAdN %>% mutate(Model = "ETS(A,Ad,N)")
china_forecast_boxcox <- china_forecast_boxcox %>% mutate(Model = "Box-Cox ETS(A,A,N)")

china_forecast_combined <- bind_rows(
  china_forecast_ANN, china_forecast_AAN, china_forecast_AAdN, china_forecast_boxcox
)

china_forecast_combined <- as_tsibble(china_forecast_combined, index = Year, key = Model)

autoplot(china_gdp, GDP) +
  geom_line(data = china_forecast_combined, aes(x = Year, y = .mean, color = Model), linewidth = 1)+
  ggtitle("Comparison of ETS Model Forecasts for Chinese GDP") +
  ylab("GDP (USD Billions)") +
  xlab("Year") +
  theme_minimal()

It seems models with “A” (Additive Trend) predict continued GDP growth, as ANN is the only one showing a flat line (no growth), which isn’t realistic.

Also, damped trends (Ad) seems to assume growth slows over time, as opposed to the Box-Cox which stabilizes variance but also can introduce strong curvature in forecast.

Personally, I would say that for this dataset ETS(A,Ad,N) (Damped trend) is the most realistic for forecasting long-term GDP.

2.4 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?

Let’s load and prepare this dataset:

gas_data <- aus_production %>%
  select(Quarter, Gas) %>%
  as_tsibble(index = Quarter)


autoplot(gas_data, Gas) +
  ggtitle("Australian Quarterly Gas Production") +
  ylab("Gas Production (Petajoules)") +
  xlab("Year") +
  theme_minimal()

So far we can see from the plot that the data has a strong seasonal pattern, and the magnitude of seasonal fluctuations increases over time, suggesting multiplicative seasonality.

This can guide us in choosing a model, since seasonality is not only clearly present but also the variation increases over time, we need to try ETS models with multiplicative seasonality.

If we used additive seasonality it would assume seasonal variations remain constant over time and would not capture the increasing magnitude of fluctuations seen in later years.

So the next step is fitting the models and forecasting.

gas_ets_AANM <- gas_data %>% model(ETS(Gas ~ error("A") + trend("A") + season("M"))) # Multiplicative
gas_ets_AAdM <- gas_data %>% model(ETS(Gas ~ error("A") + trend("Ad") + season("M")))  # Multiplicative damped 

gas_forecast_AANM <- gas_ets_AANM %>% forecast(h = "8 years")
gas_forecast_AAdM <- gas_ets_AAdM %>% forecast(h = "8 years")

And then we plot

gas_forecast_AANM <- gas_forecast_AANM %>% mutate(Model = "ETS(A,A,M)")
gas_forecast_AAdM <- gas_forecast_AAdM %>% mutate(Model = "ETS(A,Ad,M)")

gas_forecast_combined <- bind_rows(gas_forecast_AANM, gas_forecast_AAdM)


gas_forecast_combined <- as_tsibble(gas_forecast_combined, index = Quarter, key = Model)

autoplot(gas_data, Gas) +
  geom_line(data = gas_forecast_combined, aes(x = Quarter, y = .mean, color = Model), linewidth = 1) +
  ggtitle("Comparison of ETS Models for Australian Gas Production") +
  ylab("Gas Production (Petajoules)") +
  xlab("Year") +
  theme_minimal()

From the plot we can see that when using a damped trend (blue), growth slows over time and future gas production does not rise exponentially. This makes sense because physical and economic limits can restrict unlimited gas production growth, making it more realistic for long-term forecasting.

2.5 8.8

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

Once again this retail time series making a comeback.

So let’s load and prepare the data used back then (random seed and everything).

# seed
set.seed(86863)

# Use same random time series from aus_retail
myseries <- aus_retail %>%
  filter(`Series ID` == sample(aus_retail$`Series ID`, 1))

glimpse(myseries)

Rows: 441
Columns: 5
Key: State, Industry [1]
$ State       <chr> "New South Wales", "New South Wales", "New South Wales", "…
$ Industry    <chr> "Takeaway food services", "Takeaway food services", "Takea…
$ `Series ID` <chr> "A3349792X", "A3349792X", "A3349792X", "A3349792X", "A3349…
$ Month       <mth> 1982 Apr, 1982 May, 1982 Jun, 1982 Jul, 1982 Aug, 1982 Sep…
$ Turnover    <dbl> 85.4, 84.8, 80.7, 82.4, 80.7, 82.1, 87.3, 87.4, 97.2, 93.0…

autoplot(myseries, Turnover) +
  labs(title = "Monthly Retail Turnover",
       x = "Year", y = "Turnover") +
  theme_minimal()

a) Why is multiplicative seasonality necessary for this series?

For one, seasonal variations increase over time. If we used additive seasonality, we would assume the seasonal variations remain constant over time, which clearly isn’t the case here.

But also retail sales exhibit proportional growth, and multiplicative seasonality allows seasonal patterns to scale appropriately with the increasing trend.

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

Fitting the models:

retail_ets_AAM <- myseries %>%
  model(ETS(Turnover ~ error("A") + trend("A") + season("M")))

retail_ets_AAdM <- myseries %>%
  model(ETS(Turnover ~ error("A") + trend("Ad") + season("M")))  # Damped trend

Forecasting:

retail_forecast_AAM <- retail_ets_AAM %>% forecast(h = "5 years")
retail_forecast_AAdM <- retail_ets_AAdM %>% forecast(h = "5 years")

And plotting:

retail_forecast_AAM <- retail_forecast_AAM %>% mutate(Model = "ETS(A,A,M)")
retail_forecast_AAdM <- retail_forecast_AAdM %>% mutate(Model = "ETS(A,Ad,M)")

retail_forecast_combined <- bind_rows(retail_forecast_AAM, retail_forecast_AAdM)

retail_forecast_combined <- as_tsibble(retail_forecast_combined, index = Month, key = Model)

autoplot(myseries, Turnover) +
  geom_line(data = retail_forecast_combined, aes(x = Month, y = .mean, color = Model), linewidth = 1) +
  ggtitle("Comparison of Holt-Winters Multiplicative Models for Retail Turnover") +
  ylab("Turnover") +
  xlab("Year") +
  theme_minimal()

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

accuracy_AAM <- accuracy(retail_ets_AAM)
accuracy_AAdM <- accuracy(retail_ets_AAdM)

rmse_AAM <- accuracy_AAM %>% filter(.model == "ETS(Turnover ~ error(\"A\") + trend(\"A\") + season(\"M\"))") %>% pull(RMSE)
rmse_AAdM <- accuracy_AAdM %>% filter(.model == "ETS(Turnover ~ error(\"A\") + trend(\"Ad\") + season(\"M\"))") %>% pull(RMSE)

rmse_comparison <- tibble(
  Model = c("ETS(A,A,M)", "ETS(A,Ad,M)"),
  RMSE = c(rmse_AAM, rmse_AAdM)
)

rmse_comparison %>%
  kable(caption = "RMSE Comparison: Multiplicative vs. Damped") %>%
  kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover"))

RMSE Comparison: Multiplicative vs. Damped
Model	RMSE
ETS(A,A,M)	9.719140
ETS(A,Ad,M)	9.509444

For this data-set the damped trend model produces more accurate forecasts for one-step-ahead predictions as it accounts for the possibility of slowing growth in retail trends.

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

In this case it means checking the residuals of the ETS (A,Ad,M) method.

residuals_AAdM <- augment(retail_ets_AAdM) %>%
  select(Month, .resid)

Plotting residuals:

# Residual time series plot
autoplot(residuals_AAdM, .resid) +
  ggtitle("Residuals from ETS(A,Ad,M) Model") +
  ylab("Residuals") +
  xlab("Year") +
  theme_minimal()

# Histogram of residuals
ggplot(residuals_AAdM, aes(x = .resid)) +
  geom_histogram(bins = 30, fill = "lightblue", color = "black") +
  ggtitle("Histogram of Residuals") +
  xlab("Residuals") +
  theme_minimal()

# ACF plot to check for autocorrelation
residuals_AAdM %>%
  ACF(.resid) %>%
  autoplot() +
  ggtitle("Autocorrelation of Residuals") +
  theme_minimal()

# Perform Ljung-Box test
ljung_box_test <- residuals_AAdM %>%
  features(.resid, ljung_box, lag = 10)

ljung_box_test

# A tibble: 1 × 2
  lb_stat lb_pvalue
    <dbl>     <dbl>
1    17.2    0.0709

Residuals over time fluctuate around zero, with no clear trend suggesting the model is capturing most of the structure in the data.

The histogram shows residuals as approximately normally distributed.

ACF plot shows most bars falling within the blue significance limits, meaning no strong autocorrelation.

Finally, the p-value is bigger than 0.05, which fails to reject null hypothesis, confirming that the residuals behave like white noise.

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?

To find the RMSE we have to train the model up to 2010 first.

# Training set
train_data <- myseries %>%
  filter(Month <= yearmonth("2010 Dec"))

# Test set
test_data <- myseries %>%
  filter(Month > yearmonth("2010 Dec"))

# Train the model 
retail_ets_train <- train_data %>%
  model(ETS(Turnover ~ error("A") + trend("Ad") + season("M")))

Then forecast for ETS (A,Ad,M) model.

ets_forecast <- retail_ets_train %>%
  forecast(h = nrow(test_data))  # Forecast same length as test set

We will also compute the SNAIVE, assuming that the forecast is equal to the value from the same season in the previous year.

# seasonal naive model
seasonal_naive_model <- train_data %>%
  model(SNAIVE(Turnover ~ lag("year")))

# seasonal naive forecasts
snaive_forecast <- seasonal_naive_model %>%
  forecast(h = nrow(test_data))

Finally we compute the RMSE for both models

# accuracy for ETS model
accuracy_ets <- accuracy(ets_forecast, test_data)

# accuracy for SNAIVE model
accuracy_snaive <- accuracy(snaive_forecast, test_data)

# Extract RMSE values
rmse_ets <- accuracy_ets %>% pull(RMSE)
rmse_snaive <- accuracy_snaive %>% pull(RMSE)

rmse_comparison <- tibble(
  Model = c("ETS(A,Ad,M)", "Seasonal Naïve"),
  RMSE = c(rmse_ets, rmse_snaive)
)

rmse_comparison %>%
  kable(caption = "Test Set RMSE: ETS(A,Ad,M) vs. Seasonal Naïve") %>%
  kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover"))

Test Set RMSE: ETS(A,Ad,M) vs. Seasonal Naïve
Model	RMSE
ETS(A,Ad,M)	78.32102
Seasonal Naïve	96.80193

The ETS(A,Ad,M) has a lower RMSE (78.32 vs. 96.80) than the SNAIVE model. This means the ETS model forecasts retail turnover more accurately than the SNAIVE approach, probably because it captures both trend and seasonality, improving predictive power.

2.6 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?

Let’s start by applying STL decomposition with Box-Cox transformation.

# Find optimal Box-Cox lambda for variance stabilization
lambda_retail <- myseries %>%
  features(Turnover, features = guerrero) %>%
  pull(lambda_guerrero)

# STL decomposition with Box-Cox transformation
Decomp <- myseries %>%
  model(STL(box_cox(Turnover, lambda_retail) ~ season(window = "periodic"), robust = TRUE)) %>%
  components()

# Plot STL decomposition
autoplot(Decomp) +
  ggtitle("STL with Box-Cox Transformation")

Now extracting the seasonally adjusted data.

myseries$Turnover_sa <- Decomp$season_adjust

Splitting the training and test data.

train_data <- myseries %>% filter(year(Month) < 2011)
test_data <- myseries %>% filter(year(Month) >= 2011)

At the same time we fit the ETS (A,Ad,M) model on seasonally adjusted data

fit_stl_ets <- train_data %>%
  model(ETS(Turnover_sa ~ error("A") + trend("Ad") + season("M")))

# Check residuals to confirm model validity
fit_stl_ets %>% gg_tsresiduals() +
  ggtitle("Residuals of ETS(A,Ad,M) Model on Seasonally Adjusted Data")

Then we forecast test set

myseries_test <- myseries %>%
  filter(year(Month) >= 2011) %>%
  select(State, Industry, Month, Turnover_sa)

fc_stl_ets <- fit_stl_ets %>%
  forecast(new_data = myseries_test)

Now, train the previous ETS(A,Ad,M) model on raw data and forecast

# Train the previous ETS(A,Ad,M) model on raw data
retail_ets_train <- train_data %>%
  model(ETS(Turnover ~ error("A") + trend("Ad") + season("M")))

# Forecast the test set with the previous ETS(A,Ad,M) model
ets_forecast <- retail_ets_train %>%
  forecast(h = nrow(test_data))

Finally we compute RMSE and compare to the previous model

# Compute RMSE for previous ETS(A,Ad,M) model on test set
rmse_ets_test <- ets_forecast %>%
  accuracy(test_data) %>%
  filter(.type == "Test") %>%
  pull(RMSE)

# Compute RMSE for STL + ETS(A,Ad,M) on test set
rmse_stl_ets_test <- fc_stl_ets %>%
  accuracy(test_data) %>%
  filter(.type == "Test") %>%
  pull(RMSE)

# Compute RMSE for ETS(A,Ad,M) on training set
rmse_stl_ets <- fit_stl_ets %>%
  accuracy() %>%
  filter(.type == "Training") %>%
  pull(RMSE)

# Compare RMSE
rmse_comparison <- tibble(
  Model = c("ETS(A,Ad,M) on Raw Data", "STL + ETS(A,Ad,M)"),
  RMSE_Training = c(rmse_ets, rmse_stl_ets),
  RMSE_Test = c(rmse_ets_test, rmse_stl_ets_test)
)

# Display RMSE comparison
rmse_comparison %>%
  kable(caption = "RMSE Comparison: ETS(A,Ad,M) vs. STL + ETS(A,Ad,M)") %>%
  kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover"))

RMSE Comparison: ETS(A,Ad,M) vs. STL + ETS(A,Ad,M)
Model	RMSE_Training	RMSE_Test
ETS(A,Ad,M) on Raw Data	78.3210168	79.2024235
STL + ETS(A,Ad,M)	0.0435141	0.1558501

The results indicate a huge improvement in RMSE when using STL + ETS(A,Ad,M) compared to the raw ETS(A,Ad,M) model, meaning the new model generalizes better to unseen data.