8.1

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

pigs <- aus_livestock |>
  filter (Animal == 'Pigs' & State == 'Victoria')

pigs |> 
  autoplot(Count)

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.

pigs_fit <- pigs |>
  model(ETS(Count ~ error("A") + trend("N") + season("N")))

report(pigs_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

α = 0.3221247 ℓ0 = 100646.6

pigs_fc <- pigs_fit |>
  forecast(h = 4)

pigs_fc |>
  autoplot(pigs, level = NULL) +
  geom_line(aes(y = .fitted), col="#EEA2AD",
            data = augment(pigs_fit)) +
  guides(color = "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.

pigs_sd <- sd(residuals(pigs_fit)$.resid)
pigs_m <- pigs_fc$.mean[1]
print (paste('Computed 95% prediction interval: ', round(pigs_m -(1.95*pigs_sd),4), round(pigs_m +(1.95*pigs_sd),4)))
## [1] "Computed 95% prediction interval:  76964.4591 113408.6557"
rpi <-pigs_fc|>slice(1) |> hilo(level=95) |>pull()

print(paste('Prediction interval of R', round(rpi,4)))
## [1] "Prediction interval of R [76854.7889, 113518.326]95"

The prediction intervals are slightly different.

8.5

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

ice_exp <- global_economy |>
  filter (Country == 'Iceland')

a.

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

ice_exp |> 
  autoplot(Exports)

There is a larger variation in the first decade, with a steep decline in the late 1960s followed by an increase back to the highest levels in 1970. And then again, another (not quite as steep) decline in the early 1970s.

In the late 1970s through the late 1980s, there is variation between 30 and 40, and then from late 1980s to mid 2000s the variation is smaller - from 30 to 36 or so.

In the mid 2000s, there is a very steep increase to a new peak over 55, after which it starts dropping again.

b.

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

ice_fit <- ice_exp |>
  model(
    ETS(Exports ~ error("A") + trend("N") + season("N"))
  )
ice_fc <- ice_fit |> forecast(h = 15) 

ice_fc |>
  autoplot(ice_exp, level = NULL) +
  geom_line(aes(y = .fitted), col="#87CEFA",
            data = augment(ice_fit)) +
  guides(color = "none")

c. 

Compute the RMSE values for the training data.

accuracy(ice_fit)$RMSE
## [1] 3.474562

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.

ice_exp |>
  model(
    `ANN` = ETS(Exports ~ error("A") +
                       trend("N") + season("N")),
    `AAN` = ETS(Exports ~ error("A")+
                       trend("A") + season("N"))
  ) |>
  forecast(h = 15) |>
  autoplot(ice_exp, level = NULL) +
  guides(color = guide_legend(title = "Forecast"))

e.

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

The Simple - ETS(A,A,N) model predicts an upward trend, while the Holt’s - ETS(A,N,N) model predicts a flat trend. Given the last 4-5 years were a downward trend, I think the ETS(A,N,N) model showing a flat trend is best.

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.

ice_sd <- sd(residuals(ice_fit)$.resid)
ice_m <- ice_fc$.mean[1]
print (paste('Calculated 95% prediction interval for ETS(A,N,N) model : ', round(ice_m -(1.95*ice_sd),4), round(ice_m +(1.95*ice_sd),4)))
## [1] "Calculated 95% prediction interval for ETS(A,N,N) model :  40.1321 53.7958"
#did not save 2nd model, so doing again

ice_fit2 <- ice_exp |>
  model(
    ETS(Exports ~ error("A") + trend("A") + season("N"))
  )
ice_fc2 <- ice_fit2 |> forecast(h = 15) 
ice_sd2 <- sd(residuals(ice_fit2)$.resid)
ice_m2 <- ice_fc2$.mean[1]
print (paste('Calculated 95% prediction interval for ETS(A,A,N) model : ', round(ice_m2 -(1.95*ice_sd2),4), round(ice_m2 +(1.95*ice_sd2),4)))
## [1] "Calculated 95% prediction interval for ETS(A,A,N) model :  40.1452 53.909"
rpi2 <-ice_fc|>slice(1) |> hilo(level=95) |>pull()

print(paste('Prediction interval of R for ETS(A,N,N) model', round(rpi2,4)))
## [1] "Prediction interval of R for ETS(A,N,N) model [40.0334, 53.8945]95"
rpi3 <-ice_fc2|>slice(1) |> hilo(level=95) |>pull()

print(paste('Prediction interval of R for ETS(A,A,N) model', round(rpi3,4)))
## [1] "Prediction interval of R for ETS(A,A,N) model [39.918, 54.1362]95"

Again, the prediction intervals from the models and calculated are close, but not exactly the same.

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

CHN_gdp <- global_economy |>
  filter (Country == 'China')

CHN_gdp |> autoplot(GDP)+
  scale_y_continuous(labels = label_number(suffix = "T", scale = 1e-12))

# will try various combinations but keep seasonality =N
lmbda <- CHN_gdp|>
  features(GDP, features= guerrero) |>
  pull(lambda_guerrero)

CHN_gdp |>
  model(
    `ANN` = ETS(GDP ~ error("A") +
                       trend("N") + season("N")),
    `AAN` = ETS(GDP ~ error("A")+
                       trend("A") + season("N")),
    `AAdN` = ETS(GDP ~ error("A") +
                       trend("Ad") + season("N")),
    `MNN` = ETS(GDP ~ error("M")+
                       trend("N") + season("N")),
    `MAN` = ETS(GDP ~ error("M")+
                       trend("A") + season("N")),
    `MAdN` = ETS(GDP ~ error("M")+
                       trend("Ad") + season("N")),
    `BoxCoxAAN` = ETS(box_cox(GDP,lmbda) ~ error("A")+
                       trend("A") + season("N")),
    `BoxCoxAAdN` = ETS(box_cox(GDP,lmbda) ~ error("A")+
                       trend("Ad") + season("N"))
  ) |>
  forecast(h = 30) |>
  autoplot(CHN_gdp, level = NULL) +
  scale_y_continuous(labels = label_number(suffix = "T", scale = 1e-12))+
  guides(color = guide_legend(title = "Forecast"))

# repeating without either BoxCox so the rest is clearer
# Also, the ANN and MNN are the same, so just leaving ANN 

CHN_gdp |>
  model(
    `ANN` = ETS(GDP ~ error("A") +
                       trend("N") + season("N")),
    `AAN` = ETS(GDP ~ error("A")+
                       trend("A") + season("N")),
    `AAdN` = ETS(GDP ~ error("A") +
                       trend("Ad") + season("N")),
   `MAN` = ETS(GDP ~ error("M")+
                       trend("A") + season("N")),
    `MAdN` = ETS(GDP ~ error("M")+
                       trend("Ad") + season("N"))
  ) |>
  forecast(h = 30) |>
  autoplot(CHN_gdp, level = NULL) +
  scale_y_continuous(labels = label_number(suffix = "T", scale = 1e-12))+
  guides(color = guide_legend(title = "Forecast"))

## 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 <- aus_production |>
  select (Gas)

gas |> autoplot (Gas)

gas |>
  model(
  #first, taking out ANA, as does not seem to fit trend
  #  `ANA` = ETS(Gas ~ error("A") +
  #                     trend("N") + season("A")),
    `AAA` = ETS(Gas ~ error("A")+
                       trend("A") + season("A")),
  #next taking out AAdA
  #  `AAdA` = ETS(Gas ~ error("A") +
  #                    trend("Ad") + season("A")),
   `AAM` = ETS(Gas ~ error("A")+
                       trend("A") + season("M")),
    `AAdM` = ETS(Gas ~ error("A")+
                       trend("Ad") + season("M"))
  ) |>
  forecast(h = 30) |>
  autoplot(gas, level = NULL) 

  guides(color = guide_legend(title = "Forecast"))
## <Guides[1] ggproto object>
## 
## colour : <GuideLegend>

Seasonal multiplicative is effective here to capture the seasonality compounded by the upward trend.

8.8

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

myseries<-aus_retail |>
  filter (`Series ID`=='A3349849A')

myseries |> autoplot(Turnover)

a.

Why is multiplicative seasonality necessary for this series?

Similar to gas, we need to capture the seasonality compounded by the upward trend.

b.

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

ms_fit <- myseries |>
  model(
     'HW'= ETS(Turnover ~ error("M") + trend("A") + season("M")),
     'HW_damped'= ETS(Turnover ~ error("M") + trend("Ad") + season("M")),
  )
ms_fc <- ms_fit |> forecast(h=48)
ms_fc |>
  autoplot(myseries, level = NULL) +
  guides(colour = guide_legend(title = "Forecast"))+
  theme(legend.title = element_text(size=5),
        legend.text = element_text(size=5))

c. 

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

accuracy(ms_fit) |>
  select('.model', 'RMSE')
## # A tibble: 2 × 2
##   .model     RMSE
##   <chr>     <dbl>
## 1 HW         1.78
## 2 HW_damped  1.77

The Root Mean Squared Error (RMSE) for the damped model is smaller, meaning a better model.

d. 

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

ms_fit |>
  select (HW_damped)|>
  gg_tsresiduals()

  • The time plot of the residuals shows that the variation of the residuals is much the same across the historical data and therefore the residual variance can be treated as constant.
  • The acf shows all lags within range.
  • The consistency of variation can also be seen on the histogram of the residuals. The histogram shows that the residuals may not be normal — it appears bimodal. Consequently, forecasts from this method will probably be quite good, but prediction intervals that are computed assuming a normal distribution may be inaccurate.

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?

ms_train <- myseries |>
  filter(year(Month) < 2011)

autoplot(myseries, Turnover) +
  autolayer(ms_train, Turnover, colour = "red")

ms_fit2 <- ms_train |>
  model(
    'HW_damped'= ETS(Turnover ~ error("M") + trend("Ad") + season("M")),
    'HW'= ETS(Turnover ~ error("M") + trend("A") + season("M")),
    'SNaive' = SNAIVE(Turnover))

ms_fc2 <- ms_fit2 |> forecast(new_data = anti_join(myseries, ms_train))
## Joining with `by = join_by(State, Industry, `Series ID`, Month, Turnover)`
ms_fc2 |>
  autoplot(myseries, level = NULL) +
  guides(colour = guide_legend(title = "Forecast"))+
  theme(legend.title = element_text(size=5),
        legend.text = element_text(size=5))

At first I just did damped vs SNAIVE, but damped did not look better than SNAIVE, so I added back in Holt-Winters’ multiplicative (HW) and that does look better.

accuracy(ms_fc2, myseries) |>
  arrange(.model) |>
  select(.model, .type, RMSE )
## # A tibble: 3 × 3
##   .model    .type  RMSE
##   <chr>     <chr> <dbl>
## 1 HW        Test   6.62
## 2 HW_damped Test  11.7 
## 3 SNaive    Test   8.42

The RMSE here for HW_damped is very different than what we saw before splitting to train and test. It does NOT beat the SNaive. But the HW (which is also different than before splitting) DOES beat the SNAIVE.

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?

ms_lmbda <- myseries|>
  features(Turnover, features= guerrero) |>
  pull(lambda_guerrero)


ms_stl <- myseries |>
  model(
    STL(box_cox(Turnover,ms_lmbda) )
  ) |> # letting the trend(window=) and season(window=) both default
  components()

myseries$Turnover_sa <- ms_stl$season_adjust

ms_train2 <- myseries %>%
  filter(year(Month) < 2011)

ms_fit3 <- ms_train2 |>
  model(
    ETS(Turnover_sa ~ error("M") + trend("A") + season("M"))
  )

ms_fc3 <- ms_fit3 |> forecast(new_data = anti_join(myseries, ms_train2))
## Joining with `by = join_by(State, Industry, `Series ID`, Month, Turnover,
## Turnover_sa)`
ms_fc3 |>
  autoplot(myseries, level = NULL) +
  geom_line(aes(y = .fitted), col="#698B22",
            data = augment(ms_fit3)) +
  guides(color = "none")

accuracy(ms_fc3, myseries) |>
  select(RMSE)
## # A tibble: 1 × 1
##    RMSE
##   <dbl>
## 1 0.756

Holt-Winters’ multiplicative (HW) looks better than the ETS with the STL decomposition applied to the Box-Cox transformed series, and the RMSE reflects that.