Do exercises 8.1, 8.5, 8.6, 8.7, 8.8, 8.9 in Hyndman. (https://otexts.com/fpp3/)

library(fpp3)
library(tidyverse)

Question 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 a and l_o, and generate forecasts for the next four months.

aus_livestock_Victoria_pig <- aus_livestock %>%
  filter(State=="Victoria" & Animal=="Pigs")

aus_livestock_Victoria_pig %>%
  autoplot(Count)+
  labs(y="Count", title="Pigs Slaughtered in Victoria")  

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

report(fit_Victoria_pig)
## 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

The optimal value for alpha is 0.3221247 and the optimal value for level is 100646.6.

fc_Victoria_pig <- fit_Victoria_pig %>%
  forecast(h = 4)

fc_Victoria_pig %>%
  autoplot(aus_livestock_Victoria_pig) +
  labs(y="Count", title="Pigs Slaughtered in Victoria Forecast (Next 4 Months)") +
  guides(colour = "none")



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.

residuals <- residuals(fit_Victoria_pig)

s <- sd(residuals$.resid)

forecast <-fc_Victoria_pig$.mean[1]

lower_bound <- forecast - 1.96 * s
upper_bound <- forecast + 1.96 * s


cat("95% Prediction Interval:", lower_bound, "to", upper_bound, "\n")
## 95% Prediction Interval: 76871.01 to 113502.1
fc_Victoria_pig %>% 
  hilo(95) %>% 
  pull("95%")
## <hilo[4]>
## [1] [76854.79, 113518.3]95 [75927.17, 114445.9]95 [75042.22, 115330.9]95
## [4] [74194.54, 116178.6]95

The computed 95% prediction interval results closely approximate the interval produced by R which was lower bound 76854.79 and upper bound 113518.3.

Question 8.5

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

Colombia was chosen for analyse.

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

colombia_exports <- global_economy|>
  filter(Country=="Colombia")

colombia_exports %>%
autoplot(Exports)+
  labs(y="Annual Exports", title="Colombia Exports")



It appears that there is no consistent seasonality in the data. There are pronounced decreasing trends from 1960 to 1985, which then stabilize into a consistent pattern of growth and decline from 1985 to the present.

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

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

fc_colombia_export <- fit_colombia_exports %>%
  forecast(h = 10)

fc_colombia_export %>%
  autoplot(colombia_exports) +
  labs(y="Exports", title="Colombia Exports ANN Forecast") +
  guides(colour = "none")

c.Compute the RMSE values for the training data.

accuracy(fit_colombia_exports)  %>%
    select(c(".model", "RMSE",))
## # A tibble: 1 × 2
##   .model  RMSE
##   <chr>  <dbl>
## 1 ANN     1.55

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_colombia_exports_AAN <-colombia_exports %>%
  model(AAN =ETS(Exports ~ error("A") + trend("A") + season("N")))

fc_colombia_export_AAN <- fit_colombia_exports_AAN %>%
  forecast(h = 10)

fc_colombia_export_AAN %>%
  autoplot(colombia_exports) +
  labs(y="Exports", title="Colombia Exports AAN Forecast") +
  guides(colour = "none")

accuracy(fit_colombia_exports_AAN)  %>%
    select(c(".model", "RMSE",))
## # A tibble: 1 × 2
##   .model  RMSE
##   <chr>  <dbl>
## 1 AAN     1.56

The ANN (simple method) model adapts effectively to changes in the data over time by assigning exponentially decreasing weights to past observations. On the other hand, the AAN (Holt’s Method) model allows for various levels of trend damping, enabling it to adjust to different degrees of trend stability within the data. This flexibility makes Holt’s method more adept at capturing shifts in trend behavior.

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

Both models performed adequately considering the unpredictability of the trend in Colombia’s exports, although it appears that the ANN model performed slightly better, as evidenced by its lowest calculated RMSE.

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.

aan_rmse = 1.55567
aan_se = aan_rmse/(sqrt(length(colombia_exports)))

aan_forecast <-fc_colombia_export$.mean[1]

aan_lower_bound <- aan_forecast - 1.96 * aan_se
aan_upper_bound <- aan_forecast + 1.96 * aan_se


cat("95% Prediction Interval:", aan_lower_bound, "to", aan_upper_bound, "\n")
## 95% Prediction Interval: 13.53698 to 15.56972
fc_colombia_export %>% 
  hilo(95) %>% 
  pull("95%") %>%
  head(1)
## <hilo[1]>
## [1] [11.46739, 17.63931]95
ann_rmse = 1.547114 
ann_se = ann_rmse/(sqrt(length(colombia_exports)))

ann_forecast <-fc_colombia_export_AAN$.mean[1]

ann_lower_bound <- ann_forecast - 1.96 * ann_se
ann_upper_bound <- ann_forecast + 1.96 * ann_se


cat("95% Prediction Interval:", ann_lower_bound, "to", ann_upper_bound, "\n")
## 95% Prediction Interval: 13.56471 to 15.58627
fc_colombia_export_AAN %>% 
  hilo(95) %>% 
  pull("95%") %>%
  head(1)
## <hilo[1]>
## [1] [11.41553, 17.73546]95

Both ANN and AAN models calculated 95% prediction intervals within the same range as those produced using R.

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.

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

# Gather the data into a long format
china_economy <- global_economy%>%
  filter(Country == "China")

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

fit_china_economy <- china_economy %>%
  model(`Simple` = ETS(GDP ~ error("A") + trend("N") + season("N")),
        `Holt` = ETS(GDP ~ error("A") + trend("A") + season("N")),
        `Damped Holt` = ETS(GDP ~ error("A") + trend("Ad", phi = 0.8) + season("N")),
        `Box-Cox` = ETS(box_cox(GDP,lambda) ~ error("A") + trend("A") + season("N")),
        `Damped Box-Cox` = ETS(box_cox(GDP,lambda) ~ error("A") + trend("Ad", phi = 0.8) + season("N"))) 

fc_china_economy <- fit_china_economy %>%
  forecast(h = 15)

fc_china_economy %>%
  autoplot(china_economy, level = NULL) +
  labs(y="GDP",title="GDP in China") +
  guides(colour = guide_legend(title = "Forecast"))



Question 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 in this data set because it exhibits a seasonal pattern where the amplitude of the seasonal effect changes over time.

fit_aus_gas <- aus_production %>%
  model(
    "Additive" = ETS(Gas ~ error("A") + trend("A") + season("A")),
    "Multiplicative" = ETS(Gas ~ error("M") + trend("A") + season("M")),
    "Damped Holt" = ETS(Gas ~ error("A") +  trend("Ad", phi = 0.9) + season("M"))
  )

fc_aus_gas  <- fit_aus_gas  %>%
  forecast(h = 15)

fc_aus_gas %>% 
  autoplot(aus_production, level = NULL) +
  labs(title = "Australian Gas Production") +
  guides(colour = guide_legend(title = "Forecasts"))

accuracy(fit_aus_gas)%>%
    select(c(".model", "RMSE",))
## # A tibble: 3 × 2
##   .model          RMSE
##   <chr>          <dbl>
## 1 Additive        4.76
## 2 Multiplicative  4.60
## 3 Damped Holt     4.28

The damped model slightly improves the overall performance, which becomes more evident when calculating the RMSE.

Question 8.8

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

set.seed(12345678)
myseries <- aus_retail %>%
  filter(`Series ID` == sample(aus_retail$`Series ID`,1))


autoplot(myseries,Turnover)



a.Why is multiplicative seasonality necessary for this series? Multiplicative seasonality is necessary in this data set because there is a seasonal pattern.

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

fit_myseries <- myseries %>%
  model("Holt-Winters Multiplicative" = ETS(Turnover ~ error("M") + trend("Ad") + season("M")),
    "Damped Holt-Winters Multiplicative" = ETS(Turnover ~ error("A") +  trend("Ad", phi = 0.9) + season("M")))

fc_myseries <- fit_myseries  %>%
  forecast(h = 15)

fc_myseries %>% 
  autoplot(myseries, level = NULL) +
  guides(colour = guide_legend(title = "Forecasts"))

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

accuracy(fit_myseries)%>%
    select(c(".model", "RMSE",))
## # A tibble: 2 × 2
##   .model                              RMSE
##   <chr>                              <dbl>
## 1 Holt-Winters Multiplicative        0.616
## 2 Damped Holt-Winters Multiplicative 0.603

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

myseries %>%
  model("Damped Holt-Winters Multiplicative" = ETS(Turnover ~ error("A") +  trend("Ad", phi = 0.9) + season("M"))) %>%
  gg_tsresiduals()+
  ggtitle("Damped Holt-Winters Multiplicative")

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?

myseries_train <- myseries %>%
  filter(year(Month) < 2011)
myseries_test <- myseries |>
  filter(year(Month) >= 2011)
fit_myseries_train_test <- myseries_train %>%
  model("Naive" = NAIVE(Turnover),
    "Holt-Winters Multiplicative" = ETS(Turnover ~ error("M") + trend("Ad") + season("M")),
    "Damped Holt-Winters Multiplicative" = ETS(Turnover ~ error("A") +  trend("Ad", phi = 0.9) + season("M")))

fc_myseries_train_test <- fit_myseries_train_test  %>%
  forecast(h = nrow(myseries_test))

fc_myseries_train_test %>% 
  autoplot(myseries, level = NULL) +
  guides(colour = guide_legend(title = "Forecasts"))



The seasonal naïve approach is significantly enhanced by the Holt-Winters Multiplicative method.

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

myseries_train2 <- myseries_train

lambda_myseries_train <- myseries_train2%>%
  features(Turnover, features = guerrero) %>%
  pull(lambda_guerrero)

myseries_stl<- myseries_train2 %>%
  model(STL(box_cox(Turnover,lambda) ~ season(window = "periodic"), robust = TRUE)) %>%
  components() 

myseries_stl%>%
  autoplot()

myseries_test2 <- myseries_test

lambda_myseries_test <- myseries_test2%>%
  features(Turnover, features = guerrero) %>%
  pull(lambda_guerrero)

myseries_stl_test<- myseries_test2 %>%
  model(STL(box_cox(Turnover,lambda) ~ season(window = "periodic"), robust = TRUE)) %>%
  components() 

myseries_stl_test%>%
  autoplot()

myseries_train2$Turnover <- myseries_stl$season_adjust

myseries_test2$Turnover <- myseries_stl_test$season_adjust

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

fc_myseries_season_adjust <- fit_myseries_season_adjust  %>%
  forecast(h = nrow(myseries_test2))

fc_myseries_season_adjust %>% 
  autoplot(myseries_train, level = NULL) +
  guides(colour = guide_legend(title = "Forecasts"))