DATA 624 Homework 5

Exercise 8.1

Part A

data("aus_livestock")

# Filter for pigs in Victoria
pigs_vic <- aus_livestock %>%
  filter(Animal == "Pigs", State == "Victoria")

# estimate parameters
fit <- pigs_vic %>%
  model(ETS(Count ~ error("A") + trend("N") + season("N")))

# extract optimal parameters
params <- tidy(fit) %>%
  select(term, estimate)

# show parameters
params

## # A tibble: 2 × 2
##   term    estimate
##   <chr>      <dbl>
## 1 alpha      0.322
## 2 l[0]  100647.

# Generate forecasts for the next 4 months
fc <- forecast(fit, h = "4 months")

# Plot the forecasts
autoplot(fc, pigs_vic) + 
  labs(
    title = "Pigs Slaughtered in Victoria: SES Forecast",
    y = "Count",
    x = "Month"
  ) +
  theme(
    plot.title = element_text(hjust=0.5, face ='bold')
  )

Part B

# Prediction interval formula: y_hat +/- 1.96 * sd
c <- 1.96

# Get residual sd (standard deviation)
resid_sd <- augment(fit) %>%
  as_tibble() %>%
  summarise(s = sd(.resid, na.rm = TRUE)) %>%
  pull(s)

# first forecast mean
y_hat <- fc[1,] %>%
  pull(.mean)

# calculate the prediction interval
PI_lower <- y_hat - (c * resid_sd)
PI_upper <- y_hat + (c * resid_sd)

manual_PI <- tibble(
  Source = "Manual 95% Prediction Interval",
  "Lower Bound" = PI_lower,
  Mean = y_hat,
  "Upper Bound" = PI_upper
)

# Calculate R's prediction interval
R_PI <- fc %>%
  hilo(95) %>%
  unpack_hilo(`95%`) %>%
  as_tibble() %>%
  slice(1) 

R_PI_tibble <- tibble(
  Source = "R 95% Prediction Interval",
  "Lower Bound" = R_PI$`95%_lower`,
  "Mean" = R_PI$.mean,
  "Upper Bound" = R_PI$`95%_upper`
)

bind_rows(manual_PI, R_PI_tibble)

## # A tibble: 2 × 4
##   Source                         `Lower Bound`   Mean `Upper Bound`
##   <chr>                                  <dbl>  <dbl>         <dbl>
## 1 Manual 95% Prediction Interval        76871. 95187.       113502.
## 2 R 95% Prediction Interval             76855. 95187.       113518.

The manual prediction interval and R’s prediction interval are very close to each other, but R’s prediction interval is ~+/-16 wider than the manual one.

Exercise 8.5

Part A

data("global_economy")

# Filter for Spain Exports
spain_eco <- global_economy %>%
  filter(Country == "Spain") %>%
  select(Year, Exports) 

# Plot the series
autoplot(spain_eco, Exports) +
  labs(
    title = "Spanish Exports of Goods\nand Services (% of GDP)",
    y = "Exports (% of GDP)",
    x = "Year"
  ) +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold")
  )

According to the plot of the times series, there is a clear long-term upward trend. Since the data is annual, there is no seasonality. There are clear rises and falls around major economic periods, such as the 2008 financial crisis.

Part B

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

# Fit model
fit_ses <- spain_eco %>%
  model(ETS(Exports ~ error("A") + trend("N") + season("N")))

# Forecast
fc_ses <- forecast(fit_ses, h = "10 years")

# Plot the forecast
autoplot(fc_ses, spain_eco) +
  labs(
    title = "Spain Exports - SES Method",
    y = "Exports (% of GDP)",
    x = "Year"
  ) +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold")
  )  

Part C

# Get accuracy measures
ses_accuracy <- accuracy(fit_ses)

# Get RMSE
ses_accuracy$RMSE

## [1] 1.298316

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

# Fit models
fit <- spain_eco %>%
  model(
    SES   = ETS(Exports ~ error("A") + trend("N") + season("N")),  # ETS(A,N,N)
    Trend = ETS(Exports ~ error("A") + trend("A") + season("N"))   # ETS(A,A,N)
  )

# Forecast
fc <- forecast(fit, h = "10 years")

# Plot both SES & Trend forecasts on same plot
# Plot the forecast
autoplot(spain_eco, Exports) +
  autolayer(fc, level = 95, aes(colour = .model)) +
  labs(
    title = "Spanish Exports — ETS(A,N,N) vs ETS(A,A,N) Forecasts",
    y = "Exports (% of GDP)",
    x = "Year"
  ) +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold")
  )  

# Get accuracy measures
accuracy(fit)

## # 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 SES    Training  0.444     1.30 0.984  2.18   5.39 0.983 0.991 0.345
## 2 Trend  Training -0.000967  1.22 0.955 -0.474  5.42 0.953 0.934 0.332

The ETS(A,A,N) (Trend) model produce a much lower RMSE than the ETS(A,N,N) (SES) model. The same results are shown across the other accuracy measures as well, such as ME, MAE, MPE, MASE, and RMSSE. The only accuracy measure that has the opposite effect is MAPE, although the difference is minimal. Overall, the trend model provides a more accurate representation of Spain’s export behavior.

Part E

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

The ETS(A,N,N) (SES) model assumes that the data fluctuates around a constant mean with no long-term trend or seasonality. For this series, where there is a clear upward trajectory over time, this model fails to capture the growth pattern.

The ETS(A,A,N) (Additive Trend) model incorporates an additional trend component, which allows the forecasts to be in line with the recent direction of the data. This model better represents Spain’s long-term export growth.

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

# Get y_hat for each model
ses <- fc %>%
  as_tibble() %>%
  filter(.model == "SES") 
first_forecast_ses <- ses[1,]
yhat_ses <- first_forecast_ses$.mean

trend <- fc %>%
  as_tibble() %>%
  filter(.model == "Trend") 
first_forecast_trend <- trend[1,]
yhat_trend <- first_forecast_trend$.mean

# RMSE values from previous table
RMSE_ses <- 1.298316
RMSE_trend <- 1.222947

# 95% manual PIs
z <- 1.96

# Manual PIs:
manual_PI <- tibble(
  Method = "Manual",
  Model = c("SES", "Trend"),
  Lower = c(yhat_ses - (z*RMSE_ses), yhat_trend - (z*RMSE_trend)),
  Mean = c(yhat_ses, yhat_trend),
  Upper = c(yhat_ses + (z*RMSE_ses), yhat_trend + (z*RMSE_trend))
  )

# R PIs:
r_PI <- fc |>
  hilo(95) |>
  unpack_hilo(`95%`) |>
  as_tibble() |>
  group_by(.model) |>
  slice(1) |>
  transmute(
    Method = "R",
    Model = .model,
    Lower = `95%_lower`,
    Mean  = .mean,
    Upper = `95%_upper`
  )

bind_rows(manual_PI, r_PI)%>%
  select(-.model)

## # A tibble: 4 × 5
##   Method Model Lower  Mean Upper
##   <chr>  <chr> <dbl> <dbl> <dbl>
## 1 Manual SES    31.5  34.1  36.6
## 2 Manual Trend  32.2  34.6  37.0
## 3 R      SES    31.5  34.1  36.7
## 4 R      Trend  32.1  34.6  37.0

The manual prediction interval is tighter than R’s prediction interval, but the differences are minimal. Both methods provide a similar prediction interval for the first forecast.

Exercise 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.]

# Select Chinese GDP from global_economy
chinese_gdp <- global_economy %>%
  filter(Country == "China") %>%
  select(Year, GDP)

# plot Chinese GDP
autoplot(chinese_gdp, GDP) +
  labs(
    title = "Chinese GDP ($USD)",
    x = "Year",
    y = "GDP ($USD)"
  ) +
  theme(
    plot.title = element_text(hjust = 0.5, face = 'bold')
  )

# select large h value for forecasting 
h <- "50 years"

# Estimate box-cox lambda
lambda <- chinese_gdp %>%
  features(GDP, features=guerrero)
lambda <- lambda$lambda_guerrero
paste("Lambda: ",lambda)

## [1] "Lambda:  -0.0344628402676986"

# Comparing multiple models
fit_china <- chinese_gdp %>%
  model(
    SES = ETS(GDP ~ error("A") + trend("N")  + season("N")), # ETS(A,N,N)
    Trend = ETS(GDP ~ error("A") + trend("A")  + season("N")), # ETS(A,A,N)
    Damped = ETS(GDP ~ error("A") + trend("Ad") + season("N")), # ETS(A,Ad,N)
    Trend_BC = ETS(box_cox(GDP, lambda) ~ error("A") + trend("A")  + season("N")),
    Damped_BC = ETS(box_cox(GDP, lambda) ~ error("A") + trend("Ad") + season("N"))
  )

# Forecast (selecting large h value)
fc_china <- forecast(fit_china, h = "50 years", bias_adjust = TRUE)

## Warning: There was 1 warning in `mutate()`.
## ℹ In argument: `SES = (function (object, ...) ...`.
## Caused by warning:
## ! The `bias_adjust` argument of `forecast()` is deprecated as of fabletools
##   0.2.0.
## ℹ Please use the `point_forecast` argument instead.
## ℹ The deprecated feature was likely used in the fabletools package.
##   Please report the issue at <https://github.com/tidyverts/fabletools/issues>.

# Plot forecasts
autoplot(fc_china, chinese_gdp) +
  facet_wrap(~ .model, scales = "free_y") +
  labs(
    title = "China GDP - ETS Variants\n(Damped, Damped Box-Cox, SES, Trend, Trend Box-Cox)",
    x = "Year",
    y = "GDP ($USD)"
  ) +
  theme(
    plot.title = element_text(hjust = 0.5, face = 'bold')
  )

The Damped model includes a damping parameter that reduces the slope overtime. The forecast curve initially rises but then starts to level off. The prediction interval is wide.

The SES model is the simplest model with no trend or transformation. It flattens out after the last known observation and assumes no growth over time. It does not reflect the clear upward trend of the series.

The Trend model reflects the series’ upward trend. It produces rapid growth with wide uncertainty. It may over exaggerate future growth as it assumes the trend persists long-term.

It’s hard to see the box cox forecasts, so I will replot them here with a log transformation:

library(scales)
# Estimate box-cox lambda
lambda <- chinese_gdp %>%
  features(GDP, features=guerrero)
lambda <- lambda$lambda_guerrero

# Comparing multiple models
fit_china_bc <- chinese_gdp %>%
  model(
    Trend_BC = ETS(box_cox(GDP, lambda) ~ error("A") + trend("A")  + season("N")),
    Damped_BC = ETS(box_cox(GDP, lambda) ~ error("A") + trend("Ad") + season("N"))
  )

# Forecast (selecting large h value)
fc_china_bc <- forecast(fit_china_bc, h = "50 years", bias_adjust = TRUE)

# Plot forecasts
autoplot(fc_china_bc, chinese_gdp) +
  facet_wrap(~ .model) +
  scale_y_continuous(
    trans = "log10",
    labels = scales::label_number(scale_cut = cut_short_scale())
  ) +
  labs(
    title = "China GDP - ETS Variants\n(Damped Box-Cox, Trend Box-Cox)",
    x = "Year",
    y = "GDP ($USD)"
  ) +
  theme(
    plot.title = element_text(hjust = 0.5, face = 'bold')
  )

The Damped_BC model combines both a damping parameter and a box-cox variance-stabilizing transformation. It produces a wide prediction interval and a stable forecast. It reflects the long-term upward trend, but is conservative by leveling off over time. This seems to be the most realistic model.

The Trend_BC applies the same box-cox transformation without the damping parameter. It reflects the upward trend, but the uncertainty grows rapidly over time.

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?

# select the Gas data from aus_production
gas_data <- aus_production %>%
  select(Quarter, Gas)

# Plot the Data
autoplot(gas_data) +
  labs(
    title = "Australian Gas Production",
    x = "Quarter",
    y = "Gas (Petajoules)"
  ) +
  theme(
    plot.title = element_text(hjust = 0.5, face = 'bold')
  )

There is obvious seasonality in this time series.

# Create the model
fit_gas <- gas_data %>%
  model(
    additive = ETS(Gas ~ error("A") + trend("A") + season("A")),
    multiplicative = ETS(Gas ~ error("A") + trend("A") + season("M")),
    multiplicative_damped = ETS(Gas ~ error("A") + trend("Ad") + season("M"))
  )

# Forecast
fc_gas <- forecast(fit_gas, h = "5 years")

# Plot forecasts
autoplot(fc_gas, gas_data) +
  facet_wrap(~.model) +
  labs(
    title = "Australian Gas Production Forecasts\n(Additive vs. Multiplicative vs.\nMultiplicative-Damped Methods)"
  ) +
  theme(
    plot.title = element_text(hjust = 0.5, face = 'bold')
  )

Multiplicative seasonality is necessary because the seasonal amplitude grows with the level. Additive seasonality assumes a constant seasonal size and understates the overall growth. Multiplicative seasonality allows the seasonality to be proportional to the level.

accuracy(fit_gas) %>%
  select(-.type)

## # A tibble: 3 × 9
##   .model                     ME  RMSE   MAE    MPE  MAPE  MASE RMSSE   ACF1
##   <chr>                   <dbl> <dbl> <dbl>  <dbl> <dbl> <dbl> <dbl>  <dbl>
## 1 additive              0.00525  4.76  3.35 -4.69  10.9  0.600 0.628 0.0772
## 2 multiplicative        0.218    4.19  2.84 -0.920  5.03 0.510 0.553 0.0405
## 3 multiplicative_damped 0.548    4.22  2.81  1.32   4.11 0.505 0.556 0.0265

Additive vs. Multiplicative:

We can see that the multiplicative models outperform the additive model due to the following reasons:

• RMSE drops from the additive to the multiplicative models.
• MAE drops from the additive to the multiplicative models.
• MPE is closer to 0 in the multiplicative models vs. the additive model. Implies that the multiplicative models have lower bias.
• MAPE drops significantly, meaning there are fewer large percentage errors when the level is high.

Multiplicative vs. Multiplicative Damped:

We can see that the damped multiplicative model outperforms the non-damped multiplicative model for the following reasons:

• MAE drops from 2.844 to 2.815
• MASE drops from 0.510 to 0.505
• MAPE drops from 5.033 to 4.110

This implies that the damped model tracks the actual trajectory more closely.

In summary, the damping improves the realism of the model.

Exercise 8.8

# read excel file
retail <- readxl::read_excel("/Users/kristinlussi/Documents/MSDS/DATA 624/retail.xlsx", skip = 1)

# convert data frame to time series object, select just A3349398A time series
myts <- ts(retail[,"A3349398A"], frequency = 12, start = c(1982, 4)) %>% 
  as_tsibble(index=Month)

# Plot the time series
autoplot(myts) +
  labs(
    title = "Australian Retail Turnover",
    x= "Month",
    y = "Turnover"
  ) +
  theme(
    plot.title = element_text(hjust = 0.5, face = 'bold')
  )

Part A

Why is multiplicative seasonality necessary for this series?

The time series has seasonality where the amplitude of seasonal variations increases as the level rises. Multiplicative seasonality is necessary in order to capture these seasonal fluctuations in the forecasts. Additive seasonality would underestimate the variation at higher levels in the series.

Part B

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

# Create the model
fit_retail <- myts %>%
  model(
    HW_multiplicative = ETS(value ~ error("A") + trend("A") + season("M")),
    HW_multiplicative_damped = ETS(value ~ error("A") + trend("Ad") + season("M"))
  )

# Forecast
fc_retail <- forecast(fit_retail, h = 24)

# Plot forecasts
autoplot(fc_retail, myts, level = NULL) +
  facet_wrap(~.model) +
  labs(
    title = "Australian Retail Turnover\n(Multiplicative vs. Multiplicative Damped Methods)",
    x = "Month",
    y = "Turnover"
  ) +
  theme(
    plot.title = element_text(hjust = 0.5, face = 'bold')
  )

Part C

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

accuracy_partc <- accuracy(fit_retail)

tibble(
  Method = accuracy_partc$.model,
  RMSE = accuracy_partc$RMSE
)

## # A tibble: 2 × 2
##   Method                    RMSE
##   <chr>                    <dbl>
## 1 HW_multiplicative         29.4
## 2 HW_multiplicative_damped  29.6

The non-damped model exhibits a lower RMSE than the damped model, indicating a better short-term forecast accuracy. However, the RMSEs are very close (29.43 vs 29.63). The damped model remains valuable for predicting more realistic long-term forecasts by adding long-term stability.

However, for this problem I prefer to use the non-damped model.

Part D

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

# non-damped model
fit_hw_best <- fit_retail %>%
  select(HW_multiplicative)

# Check residuals from the best model
fit_hw_best %>%
  gg_tsresiduals()

According to the residuals plot over time, the residuals fluctuate around zero with no obvious trend or seasonality. There are a few spikes, possibly due to outliers, but there is no systemic pattern with these spikes. The ACF plot shows a few spikes outside the blue dotted lines, but they’re weak. Most of the autocorrelation bars are within the blue dotted lines. We can conclude that there is no significant autocorrelation.

According to the histogram, the residuals are approximately normally distributed.

In conclusion, we can confirm that the residuals behave like white noise. They show no obvious trend, seasonality, or autocorrelation, and are approximately normal. This means that the model is an adequate fit and captures the systemic structure of the data.

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?

# Split into training & testing sets
train <- myts %>%
  filter_index(~"2010 Dec")
test <- anti_join(myts, train, by = "index")

# Multiplicative Model
fit_hw <- train %>%
  model(
    HW_multiplicative = ETS(value ~ error("A") + trend("A") + season("M"))
  )

# Seasonal Naive
fit_snaive <- train %>%
  model(SNaive = SNAIVE(value))

# Forecast on test period
fc_hw <- forecast(fit_hw, new_data = test)
fc_snaive <- forecast(fit_snaive, new_data = test)

# Compute accuracy
hw_accuracy <- accuracy(fc_hw, test)
snaive_accuracy<- accuracy(fc_snaive, test)
bind_rows(hw_accuracy,snaive_accuracy) %>%
  select(-.type)

## # A tibble: 2 × 9
##   .model               ME  RMSE   MAE   MPE  MAPE  MASE RMSSE   ACF1
##   <chr>             <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>  <dbl>
## 1 HW_multiplicative -32.1  50.1  39.4 -1.19  1.48   NaN   NaN 0.0310
## 2 SNaive            160.  180.  160.   6.03  6.03   NaN   NaN 0.569

The Holt-Winters Multiplicative method has a much lower RMSE than the SNAIVE method (50.3010 vs 180.1991). This means that the Holt-Winters Multiplicative method produces much more accurate results on the test set than SNAIVE. So yes, we were able to beat the SNAIVE approach significantly.

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?

# Estimate box-cox lambda
lambda_train <- train %>%
  features(value, features=guerrero)
lambda_train <- lambda_train$lambda_guerrero
paste("Lambda: ",lambda_train)

## [1] "Lambda:  0.115687222101405"

# create model
fit_stl_ets <- train %>%
  model(
    STL_ETS = decomposition_model(
      STL(box_cox(value, lambda_train) ~ trend(window = 7) + season(window = "periodic"),
          robust = TRUE),
      ETS(season_adjust ~ error("A") + trend("A") + season("N"))
    )
  )

# Forecast
fc_stl_ets <- forecast(fit_stl_ets, new_data = test, bias_adjust = TRUE)

# Plot forecast
autoplot(fc_stl_ets, myts) +
  labs(
    title = "Australian Retail Turnover - STL Decomposition + ETS",
    x = "Month",
    y = "Turnover"
  ) +
  theme(
    plot.title = element_text(hjust = 0.5, face = 'bold')
  )

# Test accuracy
stl_ets_accuracy <- accuracy(fc_stl_ets, test)

bind_rows(hw_accuracy,snaive_accuracy,stl_ets_accuracy) %>%
  select(.model, RMSE)

## # A tibble: 3 × 2
##   .model             RMSE
##   <chr>             <dbl>
## 1 HW_multiplicative  50.1
## 2 SNaive            180. 
## 3 STL_ETS           172.

The Holt-Winters Multiplicative method performs the best achieving the lowest RMSE of 50.3010. The STL + ETS method in this exercise performs better than the SNAIVE method, but still worse than the Holt-Winters Multiplicative method.