Exercise 8.1

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

# Pull out the pigs slaughtered in Victoria data from the aus_livestock dataset.
pigs <- aus_livestock %>% 
  filter(Animal == 'Pigs' & State == 'Victoria')

# Plot the time series.
pigsSlaughteredInVictoriaPlot <- pigs %>%
  autoplot(Count) +
  labs(title = 'Pigs Slaughtered in Victoria Timeseries')
pigsSlaughteredInVictoriaPlot

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

# Using the ETS() function, Find the optimal values of α and ℓ0.
fit <- pigs %>%
  model(ses = ETS(Count ~ error('A') + trend('N') + season('N')))
optimal_values_report <- fit %>% report()
## 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

A.1 Find the optimal values of α and ℓ0.

According to the above report, the optimal values for α and ℓ0 are as follows:

α = 0.3221247

ℓ0 = 100646.6

A.2 Generate forecasts for the next four months.

# Generate forecasts for the next four months based on the data.
fourMonthForcast <- fit %>%
  forecast(h = 4)
fourMonthForcast
## # A fable: 4 x 6 [1M]
## # Key:     Animal, State, .model [1]
##   Animal State    .model    Month             Count  .mean
##   <fct>  <fct>    <chr>     <mth>            <dist>  <dbl>
## 1 Pigs   Victoria ses    2019 Jan N(95187, 8.7e+07) 95187.
## 2 Pigs   Victoria ses    2019 Feb N(95187, 9.7e+07) 95187.
## 3 Pigs   Victoria ses    2019 Mar N(95187, 1.1e+08) 95187.
## 4 Pigs   Victoria ses    2019 Apr N(95187, 1.1e+08) 95187.
# Plot the Forecast data. So that the forecast data is clearer, we will apply it to data from January 2017 onwards. 
fourMonthForcastPlot <- fit %>%
  forecast(h = 4) %>%
  autoplot(filter(pigs, Month >= yearmonth('2017 Jan'))) +
  labs(title = 'Four Month Forecast Data')
fourMonthForcastPlot

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

B.1 Compute a 95% prediction interval for the first forecast using ^y ± 1.96s where s is the standard deviation of the residuals.

# Get the first forecast.
yHat <- fourMonthForcast %>%
  pull(Count) %>%
  head(1)

# Get the standard deviation of the residuals.
standardDeviation <- augment(fit) %>%
  pull(.resid) %>%
  sd()

# Calculate the lower and upper confidence intervals. 
lowerCi <- yHat - 1.96 * standardDeviation
upperCi <- yHat + 1.96 * standardDeviation
results <- c(lowerCi, upperCi)
names(results) <- c('Lower', 'Upper')
results
## <distribution[2]>
##              Lower              Upper 
##  N(76871, 8.7e+07) N(113502, 8.7e+07)

The 95% prediction interval for the first forecast is from 76871 to 113502.

B.2 Compare your interval with the interval produced by R.

# Use R's hilo() function - https://www.rdocumentation.org/packages/distributional/versions/0.2.2/topics/hilo for the comparison.
hilo(fourMonthForcast$Count, 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 intervals calculated by R are slightly wider than the intervals we calcualted manually.

 

Exercise 8.5

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

# Select Bangladesh for analysis.
bangladeshExports <- global_economy %>%
  filter(Code == 'BGD')
head(bangladeshExports)
## # A tsibble: 6 x 9 [1Y]
## # Key:       Country [1]
##   Country    Code   Year         GDP Growth   CPI Imports Exports Population
##   <fct>      <fct> <dbl>       <dbl>  <dbl> <dbl>   <dbl>   <dbl>      <dbl>
## 1 Bangladesh BGD    1960 4274893913. NA        NA    9.31   10.0    48199747
## 2 Bangladesh BGD    1961 4817580184.  6.06     NA   11.7    10.8    49592802
## 3 Bangladesh BGD    1962 5081413340.  5.45     NA   10.8    10.7    51030137
## 4 Bangladesh BGD    1963 5319458351. -0.456    NA   11.7     9.98   52532417
## 5 Bangladesh BGD    1964 5386054619. 11.0      NA   14.1    10.0    54129100
## 6 Bangladesh BGD    1965 5906636557.  1.61     NA   13.4     9.66   55834038

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

# Plot the series.
bangladeshExports %>%
  autoplot(Exports)  +
  labs(title = 'Bangladesh Annual Exports')

The series displays a general downward trend from 1960 to 1975, and then displays a more or less steady upward trend until around 2012 where exports begin to drop again.

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

fit <- bangladeshExports %>%
  model(ANN = ETS(Exports ~ error('A') + trend('N') + season('N')))

bangladeshExportsForecast <- fit %>%
  forecast(h = 4)

bangladeshExportsForecast %>% autoplot(bangladeshExports) +
  labs(title = 'Bangladesh Annual Exports Forecast')

(c) Compute the RMSE values for the training data.

accuracy(fit)
## # A tibble: 1 x 11
##   Country    .model .type        ME  RMSE   MAE    MPE  MAPE  MASE RMSSE   ACF1
##   <fct>      <chr>  <chr>     <dbl> <dbl> <dbl>  <dbl> <dbl> <dbl> <dbl>  <dbl>
## 1 Bangladesh ANN    Training 0.0869  1.25 0.945 -0.891  12.0 0.983 0.991 0.0148

The RMSE value for the training data is 1.253158.

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

modelComparison <- bangladeshExports %>%
  model(
    ANN = ETS(Exports ~ error('A') + trend('N') + season('N')),
    AAN = ETS(Exports ~ error('A') + trend('A') + season('N'))
  )

accuracy(modelComparison)
## # A tibble: 2 x 11
##   Country    .model .type         ME  RMSE   MAE    MPE  MAPE  MASE RMSSE   ACF1
##   <fct>      <chr>  <chr>      <dbl> <dbl> <dbl>  <dbl> <dbl> <dbl> <dbl>  <dbl>
## 1 Bangladesh ANN    Training 0.0869   1.25 0.945 -0.891  12.0 0.983 0.991 0.0148
## 2 Bangladesh AAN    Training 0.00557  1.25 0.954 -1.90   12.2 0.993 0.989 0.0135

The AAN model results in a slightly lower RMSE which would suggest that it is a more accurate model for this data.

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

modelComparison %>%
  forecast(h = 4) %>%
  autoplot(bangladeshExports, level = NULL) +
  labs(title = 'Bangladesh Annual Exports ANN Vs AAN Forecast Model Comparison')

From the above forecast chart, it looks like the AAN model is better for forcasting this data. The ANN forcast shows a leveling off of the data which does not fit the overal trend of the data. The AAN model shows an upward trend in the data which is sits better with the existing data.

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

standardDeviation <- modelComparison %>%
  select(Country, AAN) %>%
  accuracy() %>%
  transmute(Country, standardDeviation = RMSE)

modelComparison %>%
  select(Country, AAN) %>%
  forecast(h = 1) %>%
  left_join(standardDeviation, by = 'Country') %>%
  mutate(lowerCi = Exports - 1.96 * standardDeviation,
         upperCi = Exports + 1.96 * standardDeviation) %>%
  select(Country, Exports, lowerCi, upperCi)
## # A fable: 1 x 5 [?]
## # Key:     Country [1]
##   Country       Exports    lowerCi    upperCi  Year
##   <fct>          <dist>     <dist>     <dist> <dbl>
## 1 Bangladesh N(15, 1.7) N(13, 1.7) N(18, 1.7)  2018

 

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

# Disable scientific numbers for readability purposes.
options(scipen = 999)

# Extract Chinese GDP data from the global_economy dataset.
chineseGDP <- global_economy %>%
  filter(Country == 'China')

# Create a plot of the data.
chineseGDP %>% autoplot(GDP) +
  labs(title = 'Chinese GDP')

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

# Get the optimized lambda value for BoxCox transformations.
lambda <- chineseGDP %>%
  features(GDP, features = guerrero) %>%
  pull(lambda_guerrero)


# Experiment with various ETS() options.
chineseGDPEtsOptionsComparision <- chineseGDP %>%
  model(
    ETS = ETS(GDP),
    ETSBoxCox = ETS(box_cox(GDP, lambda)),
    ETSDamped = ETS(GDP ~ trend('Ad', phi = 0.9)),
    ETSLog = ETS(log(GDP))
  )

chineseGDPEtsOptionsComparision %>%
  forecast(h = 20) %>%
  autoplot(chineseGDP, level = NULL) +
  labs(title = 'Chinese GDP ETS Forecast Options Comparison')

 

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?

# Plot the data.
aus_production %>%
  autoplot(Gas)

# Create an ETS model.
fit <- aus_production %>%
  model(fit = ETS(Gas))
report(fit)
## Series: Gas 
## Model: ETS(M,A,M) 
##   Smoothing parameters:
##     alpha = 0.6528545 
##     beta  = 0.1441675 
##     gamma = 0.09784922 
## 
##   Initial states:
##      l[0]       b[0]      s[0]    s[-1]    s[-2]     s[-3]
##  5.945592 0.07062881 0.9309236 1.177883 1.074851 0.8163427
## 
##   sigma^2:  0.0032
## 
##      AIC     AICc      BIC 
## 1680.929 1681.794 1711.389
fit %>%
  forecast(h = 4) %>%
  autoplot(aus_production)

Why is multiplicative seasonality necessary here?

multiplicative seasonality is necessary due to the fact that the seasonal variation trends upwards over time.

Experiment with making the trend damped. Does it improve the forecasts?

# Make the trend damped.
fit <- aus_production %>%
  model(fit = ETS(Gas  ~ trend('Ad', phi = 0.9)))

fit %>%
  forecast(h = 4) %>%
  autoplot(aus_production)

Comparing the damped and non damped trend plots above, making the trend damped does not appear to improve the forecast.

 

Exercise 8.8

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

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

myseries %>% autoplot(Turnover)

(a) Why is multiplicative seasonality necessary for this series?

Multiplicative seasonality is necessary for this series as the magnitude of the seasonalilty fluctuation increases over time.

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

fit <- myseries %>%
  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'))
  )

HoltWinters <- fit %>% forecast(h = 10)

HoltWinters %>% autoplot(myseries, level = NULL)

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

accuracy(fit) %>% select('.model', 'RMSE')
## # A tibble: 2 x 2
##   .model                              RMSE
##   <chr>                              <dbl>
## 1 Holt Winters Multiplicative Method  10.9
## 2 Holt Winters Damped Method          11.0

The multiplicative method has a slightly lower RMSE than that of the damped method suggesting that this may be the more accurate choice for the time series.

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

fit %>%
  select('Holt Winters Multiplicative Method') %>%
  gg_tsresiduals()

Both the risiduals plot and historgram above confirm that the residuals for the multiplicative method look like white noise with the execption of a few outliers. The ACF confirms that most of the risiduals are within bounds.

(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)

fit <- myseries_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)
  )

comparison <- anti_join(myseries, myseries_train, by = c('State', 'Industry', 'Series ID', 'Month', 'Turnover'))
forecastResults <- fit %>% forecast(comparison)

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

The above method comparison plot proves that both the Holt Winters Damped and Multiplicative methods beat the Seasonal Niave method, and that the Multiplicative method is the most accurate.

 

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?

# Get the optimized lambda value for BoxCox transformations.
lambda <- myseries_train %>%
  features(Turnover, features = guerrero) %>%
  pull(lambda_guerrero)

training_boxcox <- myseries_train %>%
  mutate(
    bc = box_cox(Turnover, lambda)
  )

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

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

accuracy(fit)
## # A tibble: 2 x 12
##   State Industry .model .type     ME  RMSE   MAE    MPE  MAPE  MASE RMSSE
##   <chr> <chr>    <chr>  <chr>  <dbl> <dbl> <dbl>  <dbl> <dbl> <dbl> <dbl>
## 1 Vict… Hardwar… STL B… Trai… 0.0958  1.97  1.42 0.0281  3.09 0.401 0.441
## 2 Vict… Hardwar… ETS B… Trai… 0.236   2.16  1.62 0.349   3.55 0.456 0.483
## # … with 1 more variable: ACF1 <dbl>
accuracy(multiplicative_best_fit)
## # A tibble: 1 x 12
##   State  Industry .model .type     ME  RMSE   MAE    MPE  MAPE  MASE RMSSE  ACF1
##   <chr>  <chr>    <chr>  <chr>  <dbl> <dbl> <dbl>  <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Victo… Hardwar… Holt … Trai… -0.453  9.86  7.22 -0.547  5.07 0.491 0.524 0.251

Looking at the RMSE values of the STL and ETS Box-Cox methods (0.04560738, and 0.04964458 respectively), we can see that both these methods are more accurate than our previous most accurate 'Holt Winters Multiplicative' method that has an RMSE of 0.6450982.