Excercise 8.1:

Consider the 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 lo , and generate forecasts for the next four months.
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.

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

report(pigs_victoria)

## 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

forec <- pigs_victoria |>
   forecast(h = 4)

forec |>
  autoplot(aus_livestock) +
  labs(y="Count", title="Number of pigs slaughtered in Victoria") +
  guides(colour = "none")

I see that optimal values of alpha is 0.3221247 and for l[0] is 100646.6

95 % interval

## Generated by R
forec %>% hilo() %>% pull(-1)

## <hilo[4]>
## [1] [76854.79, 113518.3]95 [75927.17, 114445.9]95 [75042.22, 115330.9]95
## [4] [74194.54, 116178.6]95

residuals_std <- sd(augment(pigs_victoria)$`.resid`)

# Compute confidence interval
lower_bound <- forec$.mean[1] - (residuals_std * 1.96)
upper_bound <- forec$.mean[1] + (residuals_std * 1.96)

# Display calculated confidence interval
sprintf(
  "Calculated confidence interval: [%f, %f]", 
  lower_bound,
  upper_bound
)

## [1] "Calculated confidence interval: [76871.012478, 113502.102384]"

Values provided by both the method are pretty much same. We see the small difference are due to how the standard deviation and confidence intervals were calculated.

Excercise 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.
b.Use an ETS(A,N,N) model to forecast the series, and plot the forecasts.
c.Compute the RMSE values for the training data.
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.
e.the forecasts from both methods. Which do you think is best?
f.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.

help(global_economy)

## starting httpd help server ... done

data("global_economy")
greece_exports  <- global_economy %>%
  filter(Country == "Greece") %>%
  select(Year, Exports)

head(greece_exports)

## # A tsibble: 6 x 2 [1Y]
##    Year Exports
##   <dbl>   <dbl>
## 1  1960    8.51
## 2  1961    8.65
## 3  1962    9.06
## 4  1963    9.35
## 5  1964    8.54
## 6  1965    8.36

autoplot(greece_exports, Exports) +
  labs(title = "Greece Exports",
       y = "Exports (US$ billions)",
       x = "Year")

The data reveals a consistent upward trend in Greece’s exports from the 1960s to today, signifying gradual growth over the years, especially after 2000. However, there were significant declines in export levels during specific periods, particularly in the late 1970s and early 1990s. The decline in the 1970s is linked to Greece’s dictatorship from 1967 to 1974, which disrupted the economy.

After 2000, exports experienced a notable increase, with sustained growth throughout the 2010s. This surge can be largely attributed to Greece’s entry into the Eurozone and the adoption of the euro in 2001, which enhanced trade and economic integration with other member states.

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

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

# Generate and plot forecasts
fc_ANN <- forecast(fit_ANN, h = "5 years")
autoplot(fc_ANN) +
  labs(title = "Forecasts of Greece Exports (ETS(A,N,N))", y = "Exports (in USD)", x = "Year")

c.Compute the RMSE values for the training data.

rmse_ANN <- accuracy(fit_ANN)$RMSE
cat("The RMSE for the ETS(A,N,N) model is:", rmse_ANN, "\n")

## The RMSE for the ETS(A,N,N) model is: 1.816383

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_AAN <- greece_exports %>%
  model(ETS(Exports ~ error("A") + trend("A") + season("N")))

# Generate and plot forecasts
fc_AAN <- forecast(fit_AAN, h = "5 years")
autoplot(fc_AAN) +
  labs(title = "Forecasts of Greece Exports (ETS(A,A,N))", y = "Exports (in USD)", x = "Year")

rmse_AAN <- accuracy(fit_AAN)$RMSE
cat("The RMSE for the ETS(A,A,N) model is:", rmse_AAN, "\n")

## The RMSE for the ETS(A,A,N) model is: 1.765767

ETS(A,A,N) Pros: This model effectively captures trends, resulting in more accurate forecasts for datasets with significant trends. Cons: Its increased complexity may lead to overfitting, particularly if the trend is not expected to continue.

ETS(A,N,N) Pros: This more straightforward model is easier to interpret and tends to be more resilient to data fluctuations. It can perform sufficiently well when trends are weak or when simplicity is desired for clarity. Cons: Its lack of trend-capturing capabilities may lead to less accurate forecasts in datasets where trends are present.

Given the lower RMSE of the ETS(A,A,N) model (1.765767), it seems to be the superior option for this dataset, as it effectively reflects the upward trend in Greece’s exports. However, if simplicity and interpretability are more important, the ETS(A,N,N) model could still be a reasonable choice, especially when high precision in forecasts is not critical.

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

# Generate forecasts for the next periods
forecast_ANN <- forecast(fit_ANN, h = 4)  # Forecast for the next 4 periods
forecast_AAN <- forecast(fit_AAN, h = 4)
## Calculated by R
forecast_ANN %>% hilo(95) %>% pull('95%') %>% head(1)

## <hilo[1]>
## [1] [29.5966, 36.84272]95

forecast_AAN <- fit_AAN %>% forecast(h = 4)
forecast_ANN <- fit_ANN %>% forecast(h = 4)
y_hat_AAN <- forecast_AAN$.mean[1]
rmse_AAN <- accuracy(fit_AAN)$RMSE

# Get the first forecast value and RMSE for ETS(A,N,N)
y_hat_ANN <- forecast_ANN$.mean[1]
rmse_ANN <- accuracy(fit_ANN)$RMSE

z <- 1.96

# Calculate prediction intervals for ETS(A,N,N)
lower_ANN <- y_hat_ANN - z * rmse_ANN
upper_ANN <- y_hat_ANN + z * rmse_ANN

# Calculate prediction intervals for ETS(A,A,N)
lower_AAN <- y_hat_AAN - z * rmse_AAN
upper_AAN <- y_hat_AAN + z * rmse_AAN

cat("Prediction Interval for ETS(A,N,N): [", lower_ANN, ", ", upper_ANN, "]\n")

## Prediction Interval for ETS(A,N,N): [ 29.65955 ,  36.77977 ]

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.

For the ETS(A,N,N) model, the manual calculations slightly overestimate the lower bound and underestimate the upper bound when compared to the values generated by R. Similarly, the manual calculation results in a higher lower bound and a lower upper bound relative to R.

Both methods yield a fairly close range of estimates, suggesting that our manual calculations are generally accurate, though they may not align perfectly with the intervals from R due to minor differences in RMSE calculations or the approach to interval estimation.

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

Lambda_china <- global_economy %>%
  filter(Country == "China") %>%
  features(GDP, features = guerrero) %>%
  pull(lambda_guerrero)

China_fit <- global_economy %>%
  filter(Country == "China") %>%
  model(`Standard` = ETS(GDP ~ error("A") + trend("N") + season("N")),
        `Holt's method` = ETS(GDP ~ error("A") + trend("A") + season("N")),
        `Damped Holt's method` = ETS(GDP ~ error("A") + trend("Ad", phi = 0.8) + season("N")),
        `Box-Cox` = ETS(box_cox(GDP,Lambda_china) ~ error("A") + trend("Ad") + season("N")),
        `Damped Box-Cox` = ETS(box_cox(GDP,Lambda_china) ~ error("A") + trend("Ad", phi = 0.8) + season("N"))) 

China_fc <- China_fit %>%
  forecast(h = 20)

China_fc %>%
  autoplot(global_economy, level = NULL) +
  labs(title="Chinese GDP") +
  guides(colour = guide_legend(title = "Forecast"))

In our plot, Damped Box-Cox and Damped Holt’s methods yield more conservative GDP growth forecasts than their undamped counterparts. The Box-Cox transformation (red line) leads to higher projections, suggesting faster growth under certain conditions. Holt’s method (blue line) also results in higher growth estimates by capturing the linear trend. The Standard ETS method serves as the baseline model, showing moderate growth without trend dampening or Box-Cox transformation.

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

gas_data <- aus_production %>%
  filter(Gas > 0)

aus_production %>%
  autoplot(Gas)

# Fit different ETS models including multiplicative seasonality and damped trend
fitted_models <- gas_data %>%
  model(
    simple_additive = ETS(Gas ~ error('A') + trend("A") + season("N")),
    simple_mult = ETS(Gas ~ error("M") + trend("A") + season("N")),
    seasonality_add = ETS(Gas ~ error('A') + trend("A") + season("A")),
    seasonality_mult = ETS(Gas ~ error("M") + trend("A") + season("M")),
    mult_season_damp = ETS(Gas ~ error("M") + trend("Ad") + season("M"))
  )
# forecasts for the next 10 time periods
forecasts <- fitted_models %>% forecast(h = 10)

# forecasts with the original data
forecasts %>%
  autoplot(gas_data, level = NULL) +
  labs(y = "Gas Production", title = "Forecasts of Australian Gas Production") +
  facet_wrap(~ .model, scales = "free_y")

accuracy(fitted_models)

## # 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 simple_additive  Traini…  0.442   15.7  11.7   1.05  14.3  2.11  2.07   0.118 
## 2 simple_mult      Traini… -0.0787  16.4  11.9  -1.32  13.5  2.14  2.16   0.0888
## 3 seasonality_add  Traini…  0.00525  4.76  3.35 -4.69  10.9  0.600 0.628  0.0772
## 4 seasonality_mult Traini… -0.115    4.60  3.02  0.199  4.08 0.542 0.606 -0.0131
## 5 mult_season_damp Traini… -0.00439  4.59  3.03  0.326  4.10 0.544 0.606 -0.0217

In the seasonality_mult and mult_season_damp models, multiplicative seasonality is essential, as evidenced by lower RMSE, MAE, and MAPE values compared to additive models (seasonality_add and simple_additive). This indicates that multiplicative seasonality better captures the seasonal patterns in gas production, which show that seasonal fluctuations increase with production.

The mult_season_damp model has a slightly lower RMSE (4.591840) than the seasonality_mult model (4.595113), suggesting that while gas production continues to rise, it may start to stabilize, making a damped trend more suitable for long-term forecasts. Overall, based on performance metrics—especially RMSE and MAPE—the mult_season_damp model is the best choice, as it effectively combines the advantages of multiplicative seasonality with a damped trend for more accurate gas production forecasts.

Excercise 8.8:

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

a.Why is multiplicative seasonality necessary for this series?
b.Apply Holt-Winters’ multiplicative method to the data. Experiment with making the trend damped.
c.Compare the RMSE of the one-step forecasts from the two methods. Which do you prefer?
d.Check that the residuals from the best method look 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?

# Set seed for reproducibility
set.seed(12345)

# Select a random retail series
myseries <- aus_retail |>
  filter(`Series ID` == sample(aus_retail$`Series ID`, 1))

# View the selected series
print(myseries)

## # A tsibble: 441 x 5 [1M]
## # Key:       State, Industry [1]
##    State             Industry                      `Series ID`    Month Turnover
##    <chr>             <chr>                         <chr>          <mth>    <dbl>
##  1 Western Australia Clothing, footwear and perso… A3349825J   1982 Apr     28.8
##  2 Western Australia Clothing, footwear and perso… A3349825J   1982 May     32.1
##  3 Western Australia Clothing, footwear and perso… A3349825J   1982 Jun     28.5
##  4 Western Australia Clothing, footwear and perso… A3349825J   1982 Jul     29  
##  5 Western Australia Clothing, footwear and perso… A3349825J   1982 Aug     25.3
##  6 Western Australia Clothing, footwear and perso… A3349825J   1982 Sep     26.9
##  7 Western Australia Clothing, footwear and perso… A3349825J   1982 Oct     29.1
##  8 Western Australia Clothing, footwear and perso… A3349825J   1982 Nov     32.6
##  9 Western Australia Clothing, footwear and perso… A3349825J   1982 Dec     47.9
## 10 Western Australia Clothing, footwear and perso… A3349825J   1983 Jan     26  
## # ℹ 431 more rows

a.Why is multiplicative seasonality necessary for this series?

set.seed(1234)

myseries <- aus_retail %>%
  filter(`Series ID` == sample(aus_retail$`Series ID`,1))
myseries %>% autoplot()

## Plot variable not specified, automatically selected `.vars = Turnover`

Multiplicative seasonality is essential in this case because seasonal variations increase as the overall level of the series rises. This indicates that the seasonal patterns are proportional to the level of turnover, doing a multiplicative model more suitable for capturing these evolving seasonal effects over time.

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

# Fit the Holt-Winters’ multiplicative models
fit <- myseries %>% model(
  Multi_Season = ETS(Turnover ~ error("M") + trend("A") + season("M")),
  Multi_Season_Damp = ETS(Turnover ~ error("M") + trend("Ad") + season("M"))
)

# Forecast the next 32 months
fc <- fit %>% forecast(h = 32)

#forecasts for both models along with the original data
fc %>%
  autoplot(myseries, level = NULL) +
  labs(y = "AUD") +
  theme_minimal()

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

accuracy(fit) %>%
  select(.model, RMSE)

## # A tibble: 2 × 2
##   .model             RMSE
##   <chr>             <dbl>
## 1 Multi_Season       1.34
## 2 Multi_Season_Damp  1.36

The Multiplicative Seasonal Model has a lower RMSE (3.045007) compared to the Multiplicative Seasonal with Damped Trend Model (3.074567). This indicates that the Multiplicative Seasonal Model is a better fit for the data, as it has smaller forecast errors.

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

best_fit <- myseries %>% model(
  Multi_Season = ETS(Turnover ~ error("M") + trend("A") + season("M"))
)
gg_tsresiduals(best_fit)

The residuals appear to fluctuate around zero with no clear trend or pattern.In the ACF plot, the bars mostly fall within the confidence bands.

The histogram appears roughly centered around zero with no significant skewness, which is another characteristic of normally distributed residuals.

set.seed(1244)
train <- myseries %>%
  filter(year(Month) < 2011)

# Test data: from 2011 onwards
test <- myseries %>%
  filter(year(Month) >= 2011)

fit <- train %>%
  model(
    'Holt Winters Multiplicative Method' = ETS(Turnover ~ error('M') + trend('A') + season('M')),
    'Holt Winters Damped Method' = ETS(Turnover ~ error('M') + trend('Ad') + season('M')),
    'Seasonal Naive' = SNAIVE(Turnover)
  )

comp <- anti_join(myseries, train, by = c('State', 'Industry', 'Series ID', 'Month', 'Turnover'))
fc <- fit %>% forecast(comp)

autoplot(comp, Turnover) +
  autolayer(fc, level = NULL) +
  labs(title = 'Forecast Method Comparison')

# Extract forecast values for each model and convert to a tibble
forecast_values <- fc %>%
  as_tibble() %>%
  select(Month, .model, .mean)

# Convert test data to tibble
test_tibble <- test %>%
  select(Month, Turnover) %>%
  as_tibble()

# Merge forecast values with actual test values by 'Month'
comparison <- test_tibble %>%
  left_join(forecast_values, by = "Month")

# Calculate squared errors between actual and forecast 
comparison <- comparison %>%
  mutate(squared_error = (Turnover - .mean)^2)

# Calculate RMSE manually for each model
rmse_manual <- comparison %>%
  group_by(.model) %>%
  summarise(RMSE = sqrt(mean(squared_error, na.rm = TRUE)))

rmse_manual

## # A tibble: 3 × 2
##   .model                              RMSE
##   <chr>                              <dbl>
## 1 Holt Winters Damped Method          8.60
## 2 Holt Winters Multiplicative Method  3.99
## 3 Seasonal Naive                      9.13

The Seasonal Naive model has the lowest RMSE at 6.12.

The Seasonal Naive model outperformed the Holt Winters models in terms of RMSE.

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

lambda_train <- train %>%
  features(Turnover, features = guerrero) %>%
  pull(lambda_guerrero)

box_cox_myseries <- train %>%
  mutate(
    bc = box_cox(Turnover, lambda_train)
  )

fit2 <- box_cox_myseries %>%
  model(
    'STL Box-Cox' = STL(bc ~ season(window = 'periodic'), robust = TRUE),
    'ETS Box-Cox' = ETS(bc)
  )

fit3 <- box_cox_myseries %>%
  model(
    'Holt Winters Multiplicative Method' = ETS(Turnover ~ error('M') + trend('A') + season('M'))
  )

accuracy(fit2)

## # 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 Tasmania Cafes, … STL B… Trai…  0.00555 0.126 0.0943  0.0299  2.04 0.329 0.351
## 2 Tasmania Cafes, … ETS B… Trai… -0.00673 0.141 0.105  -0.276   2.28 0.367 0.392
## # ℹ 1 more variable: ACF1 <dbl>

The STL Box-Cox method significantly outperformed all previous methods with an RMSE of 0.0562, much higher forecasting accuracy.

ETS Box-Cox method also performed well with an RMSE of 0.0676. The Seasonal Naive method had the next best performance with an RMSE of 6.120023, which is still considerably higher than the Box-Cox methods.

In conclusion, the transition from Holt-Winters methods and Seasonal Naive to STL decomposition with Box-Cox transformation has resulted in a remarkable improvement in forecast accuracy.

Data 624 hw5

Deepa Sharma

2024-10-06