#chunks
knitr::opts_chunk$set(eval=TRUE, message=FALSE, warning=FALSE, fig.height=5, fig.align='center')

#libraries
library(tidyverse)
library(fpp3)
library(latex2exp)
library(seasonal)
library(GGally)
library(gridExtra)
library(reshape2)
library(Hmisc)
library(corrplot)
library(e1071)
library(caret)
library(VIM)
library(forecast)
#random seed
set.seed(42)

Exercise 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 ℓ_0, and generate forecasts for the next four months.

The dataset contains monthly records from July 1972 to December 2018 (558 observations). The variables are Month, Animal (pigs), State (Victoria), and Count (the number of pigs slaughtered).

A simple exponential smoothing model was applied using the ETS() function with the specification ETS(A,N,N): Additive errors (A), no trend (N), no seasonality (N). The model estimation produced the smoothing parameter α = 0.3221247 and initial level ℓ_0 = 100646.6. This α value means that the model updates the forecast with new data, but it still relies heavily on previous data to make predictions (since α is closer to 0 than to 1).

Next, we used the model to generate forecasts for the next four months (Jan to Apr 2019). The forecasts are centered on 95,187 slaughtered pigs, with increasing variability as time passes. We used a model with an additive error component but no trend or seasonal components, so the forecasts will remain constant over time, centered on the same value 95,187. The plot below depicts the observed data as well as the fitted and forecasted values for the next four months.

#Filter data
vict_pigs_df <- aus_livestock %>% 
  filter(State == "Victoria", Animal == "Pigs")

#Check data
head(vict_pigs_df)
## # A tsibble: 6 x 4 [1M]
## # Key:       Animal, State [1]
##      Month Animal State     Count
##      <mth> <fct>  <fct>     <dbl>
## 1 1972 Jul Pigs   Victoria  94200
## 2 1972 Aug Pigs   Victoria 105700
## 3 1972 Sep Pigs   Victoria  96500
## 4 1972 Oct Pigs   Victoria 117100
## 5 1972 Nov Pigs   Victoria 104600
## 6 1972 Dec Pigs   Victoria 100500
str(vict_pigs_df)
## tbl_ts [558 × 4] (S3: tbl_ts/tbl_df/tbl/data.frame)
##  $ Month : mth [1:558] 1972 Jul, 1972 Aug, 1972 Sep, 1972 Oct, 1972 Nov, 1972 Dec...
##  $ Animal: Factor w/ 7 levels "Bulls, bullocks and steers",..: 6 6 6 6 6 6 6 6 6 6 ...
##  $ State : Factor w/ 8 levels "Australian Capital Territory",..: 7 7 7 7 7 7 7 7 7 7 ...
##  $ Count : num [1:558] 94200 105700 96500 117100 104600 ...
##  - attr(*, "key")= tibble [1 × 3] (S3: tbl_df/tbl/data.frame)
##   ..$ Animal: Factor w/ 7 levels "Bulls, bullocks and steers",..: 6
##   ..$ State : Factor w/ 8 levels "Australian Capital Territory",..: 7
##   ..$ .rows : list<int> [1:1] 
##   .. ..$ : int [1:558] 1 2 3 4 5 6 7 8 9 10 ...
##   .. ..@ ptype: int(0) 
##   ..- attr(*, ".drop")= logi TRUE
##  - attr(*, "index")= chr "Month"
##   ..- attr(*, "ordered")= logi TRUE
##  - attr(*, "index2")= chr "Month"
##  - attr(*, "interval")= interval [1:1] 1M
##   ..@ .regular: logi TRUE
#Simple exponential smoothing
fit <- vict_pigs_df %>%
  model(ETS(Count ~ error("A") + trend("N") + season("N")))

#Optimal values: alpha=0.3221247, l0=100646.6
report(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
#Forecast 4 months
fc <- fit %>%
  forecast(h = "4 months")

#Show results of forecasting
print(fc)
## # 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 "ETS(Count ~ error(\"A\") +… 2019 Jan N(95187, 8.7e+07) 95187.
## 2 Pigs   Victoria "ETS(Count ~ error(\"A\") +… 2019 Feb N(95187, 9.7e+07) 95187.
## 3 Pigs   Victoria "ETS(Count ~ error(\"A\") +… 2019 Mar N(95187, 1.1e+08) 95187.
## 4 Pigs   Victoria "ETS(Count ~ error(\"A\") +… 2019 Apr N(95187, 1.1e+08) 95187.
#Plot data + forecast and fitted values
fc |>
  autoplot(vict_pigs_df) +
  geom_line(aes(y = .fitted), col="#D55E00",
            data = augment(fit)) +
  labs(y="Count", title="Number of pigs slaughtered in Victoria state together with the forecast, \nJul 1976 - Apr 2019") +
  guides(colour = "none")

b.

Compute a 95% prediction interval for the first forecast using ŷ±1.96s where s is the standard deviation of the residuals. Compare your interval with the interval produced by R.

The manually calculated 95% prediction interval was [76,871.01, 113,502.1]. R’s hilo() function produced an interval of [76,854.79, 113,518.3]. This close alignment indicates that the manual method accurately captures forecast variability, confirming the interval produced by R’s built-in functions. The intervals produced by R are slightly wider, indicating that the forecast package handles uncertainty more nuancedly.

#Extract standard deviation of the residuals from the model
s <- augment(fit) %>%
  pull(.resid) %>%
  sd()

#Extract first fc value
first_fc <- fc %>%
  filter(Month == yearmonth("2019 Jan")) %>%
  pull(.mean)

#Calculate the 95% prediction interval for the first forecast using ŷ±1.96s
lower_bound <- first_fc - 1.96 * s
upper_bound <- first_fc + 1.96 * s

#The interval produced by R
r_interval <- fc %>%
  filter(Month == yearmonth("2019 Jan")) %>%
  mutate(interval = hilo(Count, level = 95)) %>%
  mutate(lower = interval[[1]]$lower, upper = interval[[1]]$upper)

#Print the results
cat("95% Prediction Interval (Manual Calculation): [", lower_bound, ", ", upper_bound, "]\n")
## 95% Prediction Interval (Manual Calculation): [ 76871.01 ,  113502.1 ]
cat("95% Prediction Interval (Produced by R): [", r_interval$lower, ", ", r_interval$upper, "]\n")
## 95% Prediction Interval (Produced by R): [ 76854.79 ,  113518.3 ]

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

Argentina’s exports rose steadily between 1960 and 2000, with some fluctuations.Exports increased dramatically between 2000 and 2002. This increase was primarily caused by the global rise in commodity prices, particularly for Argentina’s key exports of soybeans, wheat, and beef. Following this surge, exports began to decline noticeably in 2003. This decrease can be explained by economical crisis at the beginning of 2000s. Also, following an increase in global commodity prices, the market experienced a correction. Argentina, which is heavily reliant on agricultural exports, was impacted by these price fluctuations, which added to the decline in export revenues. Because the data is aggregated once a year, there is no explicit seasonality.

#Filter data
argentina_df <- global_economy %>% 
  filter(Country == "Argentina") 

#Check data
head(argentina_df)
## # 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 Argentina ARG    1960          NA  NA        NA    7.60    7.60   20619075
## 2 Argentina ARG    1961          NA   5.43     NA    5.99    5.99   20953077
## 3 Argentina ARG    1962 24450604878. -0.852    NA    9.38    4.69   21287682
## 4 Argentina ARG    1963 18272123664. -5.31     NA    7.89    7.89   21621840
## 5 Argentina ARG    1964 25605249382. 10.1      NA    5.56    5.56   21953929
## 6 Argentina ARG    1965 28344705967. 10.6      NA    4.15    6.23   22283390
str(argentina_df)
## tbl_ts [58 × 9] (S3: tbl_ts/tbl_df/tbl/data.frame)
##  $ Country   : Factor w/ 263 levels "Afghanistan",..: 9 9 9 9 9 9 9 9 9 9 ...
##  $ Code      : Factor w/ 263 levels "ABW","AFG","AGO",..: 8 8 8 8 8 8 8 8 8 8 ...
##  $ Year      : num [1:58] 1960 1961 1962 1963 1964 ...
##  $ GDP       : num [1:58] NA NA 2.45e+10 1.83e+10 2.56e+10 ...
##  $ Growth    : num [1:58] NA 5.428 -0.852 -5.308 10.13 ...
##  $ CPI       : num [1:58] NA NA NA NA NA NA NA NA NA NA ...
##  $ Imports   : num [1:58] 7.6 5.99 9.38 7.89 5.56 ...
##  $ Exports   : num [1:58] 7.6 5.99 4.69 7.89 5.56 ...
##  $ Population: num [1:58] 20619075 20953077 21287682 21621840 21953929 ...
##  - attr(*, "key")= tibble [1 × 2] (S3: tbl_df/tbl/data.frame)
##   ..$ Country: Factor w/ 263 levels "Afghanistan",..: 9
##   ..$ .rows  : list<int> [1:1] 
##   .. ..$ : int [1:58] 1 2 3 4 5 6 7 8 9 10 ...
##   .. ..@ ptype: int(0) 
##   ..- attr(*, ".drop")= logi TRUE
##  - attr(*, "index")= chr "Year"
##   ..- attr(*, "ordered")= logi TRUE
##  - attr(*, "index2")= chr "Year"
##  - attr(*, "interval")= interval [1:1] 1Y
##   ..@ .regular: logi TRUE
#Plot the data
argentina_df %>%
  autoplot(Exports) +
  labs(y="Exports, % of GDP", title="Exports of goods and services from Argetina,  1960 - 2017")

b.

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

The plot depicts the original series of Argentina’s exports, as well as the forecasts generated by the ETS(A, N) model for the next 4 years.

The optimal smoothing parameter (alpha) is around 0.8855. The initial level (l[0]) is around 7.4074. The model produced a sigma^2 of 7.9819, indicating variance in the residuals.

We used the model to forecast the next 4 years. The plot below shows the historical data, fitted values (in red), and forecasted values with prediction intervals. The plot indicates that the ETS(A, N, N) model predicts constant export percentages, with no upward or downward trends, does not account for potential fluctuations caused by future economic changes, as it is supposed to be for the simple exponential smoothing.

#Simple exponential smoothing
fit <- argentina_df %>%
  model(ETS(Exports ~ error("A") + trend("N") + season("N")))

#Optimal values: alpha=0.8855249, l0=7.407428
report(fit)
## Series: Exports 
## Model: ETS(A,N,N) 
##   Smoothing parameters:
##     alpha = 0.8855249 
## 
##   Initial states:
##      l[0]
##  7.407428
## 
##   sigma^2:  7.9819
## 
##      AIC     AICc      BIC 
## 359.9466 360.3911 366.1279
#Forecast 4 years
fc <- fit %>%
  forecast(h = 4)

#Plot data + forecast and fitted values
fc |>
  autoplot(argentina_df) +
  geom_line(aes(y = .fitted), col="#D55E00",
            data = augment(fit)) +
  labs(y="Exports, % of GDP", title="Exports of goods and services from Argetina \ntogether with the ETS(A,N,N) forecast, 1960 - 2019") +
  guides(colour = "none")

c.

Compute the RMSE values for the training data.

The training data has a root mean squared error value of approximately 2.776. It represents the average deviation between the model’s fitted values and the actual export.

#Accuracy of the model, rmse=2.776087
accuracy(fit)
## # A tibble: 1 × 11
##   Country   .model      .type     ME  RMSE   MAE   MPE  MAPE  MASE RMSSE    ACF1
##   <fct>     <chr>       <chr>  <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>   <dbl>
## 1 Argentina "ETS(Expor… Trai… 0.0762  2.78  1.62 -1.73  15.7 0.983 0.986 0.00902

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.

ETS(A,N,N) and ETS(A,A,N) have similar RMSE values of around 2.78. It means that adding a trend component does not significantly improve forecast accuracy. The ETS(A,N,N) has a slightly lower MAE (1.62), it appears that the series lacks a strong, consistent trend component. The MAE for ETS(A,A,N) is slightly higher (1.64), implying that the trend component may not provide much value. The difference in ACF1 values (0.00902 for ETS(A,N,N) vs. 0.0271 for ETS(A,A,N)) implies that the residuals are slightly more correlated in the ETS(A,A,N) model, which could indicate overfitting. The ETS(A,N,N) model, being simpler, appears to handle the data well without adding unnecessary complexity. Despite the addition of a parameter, the ETS(A,A,N) model improves forecasting accuracy for this dataset only marginally. Therefore, in this case, the simpler model may be more appropriate.

#ETS(A,A,N) exponential smoothing
fit_aan <- argentina_df %>%
  model(ETS(Exports ~ error("A") + trend("A") + season("N")))

#Optimal values: alpha=0.1050697, beta  = 0.1050697, l0=25.02638, b0=0.7827546
report(fit_aan)
## Series: Exports 
## Model: ETS(A,A,N) 
##   Smoothing parameters:
##     alpha = 0.8629392 
##     beta  = 1e-04 
## 
##   Initial states:
##      l[0]       b[0]
##  6.794833 0.07159122
## 
##   sigma^2:  8.2795
## 
##      AIC     AICc      BIC 
## 363.9608 365.1147 374.2630
#Forecast 4 years
fc_aan <- fit_aan %>% forecast(h = 4)

#Plot data + forecast and fitted values
fc_aan |>
  autoplot(argentina_df) +
  geom_line(aes(y = .fitted), col="#D55E00",
            data = augment(fit_aan)) +
  labs(y="Exports, % of GDP", title="Exports of goods and services from Argetina \ntogether with the ETS(A,A,N) forecast, 1960 - 2019") +
  guides(colour = "none")

#Accuracy of the model, rmse=2.776427
accuracy(fit_aan)
## # A tibble: 1 × 11
##   Country   .model      .type      ME  RMSE   MAE   MPE  MAPE  MASE RMSSE   ACF1
##   <fct>     <chr>       <chr>   <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>  <dbl>
## 1 Argentina "ETS(Expor… Trai… 0.00795  2.78  1.64 -2.51  15.9 0.994 0.986 0.0271

e.

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

The ETS(A,N,N) model produces a forecast that remains relatively flat, indicating that the model does not anticipate a clear trend moving forward. This behavior may be appropriate if the future of the time series is expected to oscillate around a stable level. The ETS(A,A,N) model generates a forecast with a slight trend component, which may not be consistent with the nature of the data. Given the similar RMSE values, lower residual autocorrelation, and the advantage of simplicity, the ETS(A,N,N) model appears to be a better fit for this dataset. The exports series does not show a consistent, strong trend that justifies the use of a more complex model such as ETS(A, A, N).

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

ETS(A,N,N) Model:

  • Manual calculation: [5.88, 16.76]
  • Produced by R: [5.79, 16.86]

ETS(A,A,N) Model:

  • Manual calculation: [5.99, 16.87]
  • Produced by R: [5.79, 17.07]

The manually calculated intervals using RMSE values differ slightly from the intervals generated by R. This difference arises because the manual calculation assumes constant variance (RMSE) and normal errors, whereas the R takes into account the internal error structure and model-specific uncertainty, including level and trend components (if present). The intervals produced by R are slightly wider, indicating that the forecast package handles uncertainty more nuancedly.

#Extract first fc value
first_fc_ann <- fc %>%
  filter(Year == 2018) %>%
  pull(.mean)

first_fc_aan <- fc_aan %>%
  filter(Year == 2018) %>%
  pull(.mean)

z_value <- qnorm(0.975) 

#Calculate the 95% prediction interval for the first forecast using rmse
lower_ann <- first_fc_ann - z_value * accuracy(fit)$RMSE
upper_ann <- first_fc_ann + z_value * accuracy(fit)$RMSE

lower_aan <- first_fc_aan - z_value * accuracy(fit_aan)$RMSE
upper_aan <- first_fc_aan + z_value * accuracy(fit_aan)$RMSE

#The interval produced by R
r_interval_ann <- fc %>%
  filter(Year == 2018) %>%
  mutate(interval = hilo(Exports, level = 95)) %>%
  mutate(lower = interval[[1]]$lower, upper = interval[[1]]$upper)

r_interval_aan <- fc_aan %>%
  filter(Year == 2018) %>%
  mutate(interval = hilo(Exports, level = 95)) %>%
  mutate(lower = interval[[1]]$lower, upper = interval[[1]]$upper)

#Print the results
cat("95% Prediction Interval for ETS(A,N,N) (Manual Calculation): [", lower_ann, ", ", upper_ann, "]\n")
## 95% Prediction Interval for ETS(A,N,N) (Manual Calculation): [ 5.882074 ,  16.76414 ]
cat("95% Prediction Interval for ETS(A,N,N) (Produced by R): [", r_interval_ann$lower, ", ", r_interval_ann$upper, "]\n")
## 95% Prediction Interval for ETS(A,N,N) (Produced by R): [ 5.785765 ,  16.86044 ]
cat("95% Prediction Interval for ETS(A,A,N) (Manual Calculation): [", lower_aan, ", ", upper_aan, "]\n")
## 95% Prediction Interval for ETS(A,A,N) (Manual Calculation): [ 5.989515 ,  16.87291 ]
cat("95% Prediction Interval for ETS(A,N,N) (Produced by R): [", r_interval_aan$lower, ", ", r_interval_aan$upper, "]\n")
## 95% Prediction Interval for ETS(A,N,N) (Produced by R): [ 5.791571 ,  17.07085 ]

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

The plot for China’s GDP has a distinct upward trend, particularly after the 1990, which suggests an exponential growth pattern. This growth is distinguished by an increase in variability over time, which a regular additive trend model may find difficult to capture. We should investigate a variety of ETS models in order to account for the variability and dynamics.

4 models were used to generate a 15-year forecast for China’s GDP:

  • Additive trend (A,A,N): It assumes linear growth with constant errors. This model serves as a baseline, capturing consistent growth patterns but potentially underestimating exponential growth. RMSE = 1.90e+11; MAE = 9.59e+10, AIC = 3258.

  • Additive damped trend (A,Ad,N): Damping moderates future growth, assuming a slowdown in China’s GDP growth rate over time. RMSE = 1.90e+11, MAE = 9.49e + 10, AIC = 3260. Close AIC values for additive trend models indicate that introducing a damped trend only slightly modifies future growth while not significantly improving the model’s fit.

  • Box-Cox transform with additive damped trend: Box-Cox transformation stabilizes variance and an additive damped trend. It helps to deal with changing variability in GDP more effectively while assuming a potential slowing of growth. RMSE = 1.96e+11, MAE = 1.02e+11, AIC = -135. The error metrics are slightly higher than those for simple additive models, indicating that growth has been overestimated. The significantly lower AIC indicates a better overall fit to the data, demonstrating the model’s ability to adapt to changing variance while maintaining growth moderation.

  • Multiplicative trend and error (M,M,N): It takes into account proportional changes in the data to better capture the exponential nature of GDP growth. However, the model may be overly optimistic about future growth. RMSE = 2.67e+11, MAE = 1.16e+11, AIC = 3096. The higher error metrics indicate that this model may overfit the data and predict overly optimistic growth. The lower than additive models AIC means the model’s aggressive projection may pose a risk of overestimation in the long run.

  • Box-Cox transform with multiplicative trend and error (Box-Cox M,M,N): It combines the benefits of variance stabilization with multiplicative components to handle changing variance and growth dynamics. RMSE = 2.97e+11, MAE = 1.29e+11, AIC = -136. The models’ highest error metrics indicate that, while they capture exponential trends, they may be overly aggressive. This model has the lowest AIC, indicating the best fit to the data, but its higher error metrics suggest that it may predict too strong a growth trajectory.

Additive models have similar accuracy metrics, with the damped trend model performing slightly better in terms of MAE. Their higher AIC values compared to the Box-Cox models suggest that they do not capture the variability in GDP growth as well. However, their lower RMSE and MAE suggest that they provide more conservative and potentially realistic forecasts. The Box-Cox damped model accounts for changing GDP variance using the Box-Cox transformation and incorporates damping to moderate future growth. While the RMSE and MAE are slightly higher than the simple additive models, the significantly lower AIC (-135) indicates a better fit to the data. The model has some balance between exponential growth and a possible future slowdown.

Multiplicative models are more aggressive in capturing exponential growth. Despite having the lowest AIC values, their higher RMSE and MAE indicate a tendency to overestimate future GDP. The Box-Cox multiplicative model, in particular, has the highest RMSE and MAE values, indicating potential risks in forecasting overly optimistic growth.

The choice of model depends on the assumptions about China’s future GDP growth. The Box-Cox damped model takes a balanced approach, capturing the changing variance in historical GDP data while also including a damping effect to moderate future growth. It outperforms simpler additive models in terms of model fit (lower AIC), but it may slightly overestimate growth due to higher error metrics. If we expect China’s GDP to continue its exponential growth but at a slower rate, the Box-Cox damped model appears to be the best fit.

#Filter data
china_df <- global_economy %>% filter(Country == "China")

#Plot data
autoplot(china_df, GDP) +
  labs(y="GDP, $USD", title="China, Gross domestic product, 1960 - 2017") 

#Box-Cox lambda
lambda <- china_df %>%
  features(GDP, features = guerrero) %>%
  pull(lambda_guerrero)

#ETS models
ets_model <- china_df %>%
  model(ets_aan = ETS(GDP ~ error("A") + trend("A") + season("N")), #additive trend
        ets_aan_damped = ETS(GDP ~ error("A") + trend("Ad") + season("N")), #damped trend
        ets_boxcox_damped = ETS(box_cox(GDP, lambda) ~ error("A") + trend("Ad") + season("N")), #boxcox damped
        ets_mmn = ETS(GDP ~ error("M") + trend("M") + season("N")), # Multiplicative trend and error
        ets_boxcox_mmn = ETS(box_cox(GDP, lambda) ~ error("M") + trend("M") + season("N"))) # Box-Cox with multiplicative trend and error
        

#accuracy
accuracy(ets_model)
## # A tibble: 5 × 11
##   Country .model        .type       ME    RMSE     MAE     MPE  MAPE  MASE RMSSE
##   <fct>   <chr>         <chr>    <dbl>   <dbl>   <dbl>   <dbl> <dbl> <dbl> <dbl>
## 1 China   ets_aan       Trai…  2.36e10 1.90e11 9.59e10  1.41    7.62 0.442 0.453
## 2 China   ets_aan_damp… Trai…  2.95e10 1.90e11 9.49e10  1.62    7.62 0.438 0.454
## 3 China   ets_boxcox_d… Trai…  6.35e 9 1.96e11 1.02e11  1.77    7.26 0.472 0.468
## 4 China   ets_mmn       Trai… -1.08e10 2.67e11 1.16e11 -0.0321  7.29 0.534 0.637
## 5 China   ets_boxcox_m… Trai… -3.25e10 2.97e11 1.29e11  0.925   7.32 0.594 0.709
## # ℹ 1 more variable: ACF1 <dbl>
report(ets_model)
## # A tibble: 5 × 10
##   Country .model     sigma2 log_lik   AIC  AICc   BIC      MSE     AMSE      MAE
##   <fct>   <chr>       <dbl>   <dbl> <dbl> <dbl> <dbl>    <dbl>    <dbl>    <dbl>
## 1 China   ets_aan  3.88e+22 -1624.  3258. 3259. 3268. 3.61e+22 1.31e+23 9.59e+10
## 2 China   ets_aan… 3.96e+22 -1624.  3260. 3262. 3273. 3.62e+22 1.30e+23 9.49e+10
## 3 China   ets_box… 1.50e- 3    73.4 -135. -133. -123. 1.37e- 3 4.56e- 3 2.95e- 2
## 4 China   ets_mmn  9.31e- 3 -1543.  3096. 3097. 3106. 7.13e+22 3.04e+23 7.20e- 2
## 5 China   ets_box… 4.84e- 6    73.2 -136. -135. -126. 1.34e- 3 3.68e- 3 1.69e- 3
#Forecast for the 15 years
fc <- ets_model %>% forecast(h = 15)

#Choose data to plot forecast clearer
china_df_1990 <- china_df %>% filter(Year >= 1990)

#Plot forecasts
fc |>
  autoplot(china_df_1990, level=NULL) +
  labs(y="GDP, $USD", title="China, Gross domestic product with forecast for 15 years") +
  scale_color_manual(values = c("red", "darkgreen", "purple", "blue", "lightblue"), labels = c("A,A,N","A,Ad,N", "Box-Cox A,Ad,N", "M,M,N", "Box-Cox M,M,N"))

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?

The plot of Australian gas production shows a clear upward trend as well as seasonal fluctuations. Multiplicative seasonality is required here because seasonal fluctuations in gas production appear to increase as overall production rises (the magnitude of the seasonal pattern in the data depends on the magnitude of the data).

3 models were used to generate a 4-year forecast for Australian gas production: Additive trend and additive seasonality (A,A,A); additive trend and multiplicative seasonality (M,A,M); damped additive trend with multiplicative seasonality (M,Ad,M). We used multiplicative errors for M,A,M and M,Ad,M models because they adapt to the changing scale of the data, which shows increasing variability with the trend. Using an additive seasonality model would fail to capture mentioned proportional changes. The multiplicative models account for this proportionality, and are a better fit to the data. Damped trend helps to moderate future growth if the data indicates that growth will eventually slow. The models with multiplicative seasonality outperform the model with additive seasonality. For example, the A,A,A model has an AIC of 1872, an MAE of 3.35, and an RMSE of 4.76, while the M,A,M model has an AIC of 1681, an MAE of 3.02, and an RMSE of 4.60.

The M,Ad,M model’s accuracy is slightly better than for M,A,M model because of the introduction of a damped trend (RMSE of 4.59 vs 4.6, MAE of 3.03 vs 3.02). However, the AIC is slightly higher (1684 vs 1681) that for the model without damped trend. The damped trend model has a slightly better accuracy, but the difference is negligible, and the increased AIC implies that it introduces additional complexity without significantly improving the forecast. As a result, the M,A,M model is the most suitable for capturing the seasonal patterns and underlying trend in this dataset.

#Filter data
gas_df <- aus_production %>% select(Quarter, Gas)

#Plot data
autoplot(gas_df, Gas) +
  labs(y="Gas production, petajoules", title="Australian Gas Production")

#ETS models
ets_model <- gas_df %>%
  model(ets_aaa = ETS(Gas ~ error("A") + trend("A") + season("A")), #additive seasonality
        ets_aam = ETS(Gas ~ error("M") + trend("A") + season("M")), #multiplicative seasonality and errors
        ets_aam_damped = ETS(Gas ~ error("M") + trend("Ad") + season("M"))) #added damped trend

#Forecast for the 16 quarters (4 years)
fc <- ets_model %>% forecast(h = 16)

#Accuracy
accuracy(ets_model)
## # A tibble: 3 × 10
##   .model         .type          ME  RMSE   MAE    MPE  MAPE  MASE RMSSE    ACF1
##   <chr>          <chr>       <dbl> <dbl> <dbl>  <dbl> <dbl> <dbl> <dbl>   <dbl>
## 1 ets_aaa        Training  0.00525  4.76  3.35 -4.69  10.9  0.600 0.628  0.0772
## 2 ets_aam        Training -0.115    4.60  3.02  0.199  4.08 0.542 0.606 -0.0131
## 3 ets_aam_damped Training -0.00439  4.59  3.03  0.326  4.10 0.544 0.606 -0.0217
report(ets_model)
## # A tibble: 3 × 9
##   .model           sigma2 log_lik   AIC  AICc   BIC   MSE  AMSE    MAE
##   <chr>             <dbl>   <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>  <dbl>
## 1 ets_aaa        23.6       -927. 1872. 1873. 1903.  22.7  29.7 3.35  
## 2 ets_aam         0.00324   -831. 1681. 1682. 1711.  21.1  32.2 0.0413
## 3 ets_aam_damped  0.00329   -832. 1684. 1685. 1718.  21.1  32.0 0.0417
#Plot forecasts
gas_df_1980 <- gas_df %>% filter(Quarter >= yearquarter("1980 Q1"))
fc |>
  autoplot(gas_df_1980, level=NULL) +
  labs(y="Gas production, petajoules", title="Australian Gas Production with forecast for 4 years") +
  facet_grid(.model~.) + 
  scale_color_manual(values = c("red", "purple", "blue"), labels = c("A,A,A", "M,A,M", "M,Ad,M"))

Exercise 8.8

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

a.

Why is multiplicative seasonality necessary for this series?

The turnover data plot shows both an upward trend and seasonality, with the seasonal variations increasing in magnitude as the series’ overall level rises. This indicates that the seasonal effects are proportional to the series level, implying that multiplicative seasonality is a good fit for this data. In contrast, an additive seasonality model assumes constant seasonal effects at all series levels, which does not account for the observed pattern.

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

#Plot data
autoplot(myseries, Turnover) +
  labs(y="Turnover,  $Million AUD", title="Tasmanian retail trade turnover from cafes, restaurants \nand takeaway food services, Apr 1982 - Dec 2018")

b.

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

To capture seasonal variations and trends, we used models below to generate a 2-year forecast for Tasmanian retail turnover. :

  • Multiplicative seasonality with additive trend (M,A,M): It assumes linear growth with seasonal variations proportional to series level. RMSE = 1.34, MAE = 0.943, AIC = 2875. Multiplicative errors are used because the variability in turnover appears to increase proportionally to the series level.

  • Multiplicative seasonality with damped additive trend (M,Ad,M): It also includes damping to slow the trend over time. RMSE = 1.36, MAE = 0.946, AIC = 2887.

The metrics show that both models are very close in accuracy, with the non-damped model slightly outperforming the damped model in terms of RMSE and AIC. The close values show that including a damped trend does not significantly improve the model’s fit. The forecast plot demonstrates that both models capture the trend and seasonality accurately. While the damped trend model modifies the future growth rate, the differences in accuracy between the two models are minor. As a result, the choice between these models may be determined by whether moderated (dampened) future growth is deemed more realistic or not. However, because the non-damped model outperforms the damped version slightly, it may be preferable for forecasting Tasmanian retail turnover in this scenario.

#Holt-Winters' multiplicative method and damped trend
ets_model <- myseries %>%
  model(ets_mam = ETS(Turnover ~ error("M") + trend("A") + season("M")),
        ets_mam_damped = ETS(Turnover ~ error("M") + trend("Ad") + season("M")))

#Forecast for the 2 years (24 months)
fc <- ets_model %>% forecast(h = 24)

#Accuracy
accuracy(ets_model)
## # 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, res… ets_m… Trai… 3.10e-4  1.34 0.943 -0.249  3.77 0.393 0.435
## 2 Tasmania Cafes, res… ets_m… Trai… 1.18e-1  1.36 0.946  0.280  3.74 0.394 0.440
## # ℹ 1 more variable: ACF1 <dbl>
report(ets_model)
## # A tibble: 2 × 11
##   State    Industry  .model  sigma2 log_lik   AIC  AICc   BIC   MSE  AMSE    MAE
##   <chr>    <chr>     <chr>    <dbl>   <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>  <dbl>
## 1 Tasmania Cafes, r… ets_m… 0.00277  -1420. 2875. 2876. 2944.  1.80  2.68 0.0378
## 2 Tasmania Cafes, r… ets_m… 0.00288  -1426. 2887. 2889. 2961.  1.84  2.80 0.0380
#Plot forecasts
myseries_2000 <- myseries %>% filter(Month >= yearmonth("2010 Jan"))
fc |>
  autoplot(myseries_2000, level=NULL) +
  labs(y="Turnover,  $Million AUD", title="Tasmanian retail trade turnover from cafes, restaurants a\nnd takeaway food services with forecast for 6 years") +
  scale_color_manual(values = c("red", "purple"), labels = c("M,A,M", "M,Ad,M"))

c.

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

The RMSE values for the single-step forecasts: Holt-Winters multiplicative is 4.63; Damped trend model: 4.35.

The damped model had a slightly lower RMSE. For a one-step (short-term) forecast, the damped model predicted the next data point in the series more accurately.

#Forecast for the 1 month
fc_onestep <- ets_model %>% forecast(h = 1)

# Calculate RMSE for one-step forecasts
rmse_mam <- sqrt(mean((tail(myseries, 1)$Turnover - fc_onestep$.mean[1])^2))
rmse_damped <- sqrt(mean((tail(myseries, 1)$Turnover - fc_onestep$.mean[2])^2))

cat("Holt-Winters One-Step RMSE:", rmse_mam, "\n")
## Holt-Winters One-Step RMSE: 4.630392
cat("Damped Trend Model One-Step RMSE:", rmse_damped, "\n")
## Damped Trend Model One-Step RMSE: 4.354338

d.

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

Although the damped model has slightly higher one-step forecast accuracy, the performance differences between the two models are small. The non-damped model slightly outperforms in terms of overall RMSE and AIC, so it is reasonable to choose it for this assignment.

Residuals over time plot shows residuals fluctuating around zero without any discernible pattern or trend. ACF plot has the majority of the bars fall within the blue dashed lines, indicating that there is no significant autocorrelation in the residuals. One spike exceeds the confidence intervals, but it is within acceptable limits. The histogram of residuals shows a roughly symmetric distribution centered on zero. This implies that the residuals are normally distributed. The Ljung-Box test showed a p-value of 0.410, which is higher than the typical significance level (0.05). This implies that we cannot reject the null hypothesis that the residuals are white noise. There is no significant autocorrelation in the residual data.

Based on the residual plots and the Ljung-Box test results, it appears that the residuals resemble white noise, the chosen non-damped is appropriate for the data.

#Check residuals
ets_model %>% select(ets_mam) %>%
  gg_tsresiduals() +
  ggtitle("Residuals' plots, Australian Turnover")

#Ljung-Box test 
ets_model %>% select(ets_mam) %>%
  augment(ets_model) %>%  features(.innov, ljung_box, lag = 10)
## # A tibble: 1 × 3
##   .model  lb_stat lb_pvalue
##   <chr>     <dbl>     <dbl>
## 1 ets_mam    10.4     0.410

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?

The SNaive model produced a test set RMSE of 9.13, an MAE of 7.58, and a MAPE of 14.4%. The ETS(M, A, M) model had a test set RMSE of 3.99, MAE of 3.23, and MAPE of 6.31%. In contrast, the damped model had a higher RMSE (8.60), MAE (7.18), and MAPE (13.6%). The ETS(M, A, M) model outperforms the seasonal naive and damped models based on these metrics. The forecast plot shows that the ETS(M, A, M) model (in red) closely tracks the upward trend of retail turnover, whereas the damped ETS(M, Ad, M) model (in purple) generates more conservative forecasts. In contrast, the SNaive model (in blue) reproduces the seasonal pattern observed during the training period but fails to capture the upward trend effectively, as evidenced by its higher RMSE, MAE, and MAPE values. Furthermore, residual analysis revealed that the residuals of the ETS(M, A, M) model were low in autocorrelation, indicating that the model accurately captured the data’s underlying patterns.

#naïve approach from Exercise 7 in Section 5.11
myseries_train <- myseries |>
  filter(year(Month) < 2011)

#MAM models
ets_models <- myseries_train %>%
  model(ets_mam = ETS(Turnover ~ error("M") + trend("A") + season("M")),
        ets_mam_damped = ETS(Turnover ~ error("M") + trend("Ad") + season("M")))

#SNAIVE model
snaive_model <- myseries_train %>%
  model(SNAIVE(Turnover))

fc_ets <- ets_models |>
  forecast(new_data = anti_join(myseries, myseries_train))

fc_snaive <- snaive_model |>
  forecast(new_data = anti_join(myseries, myseries_train))

#Accuracy
ets_models |> accuracy()
## # 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, re… ets_m… Trai… -0.00758  1.18 0.823 -0.305  4.04 0.371 0.405
## 2 Tasmania Cafes, re… ets_m… Trai…  0.111    1.18 0.815  0.316  3.98 0.368 0.406
## # ℹ 1 more variable: ACF1 <dbl>
snaive_model |> accuracy()
## # A tibble: 1 × 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 Tasmania Cafes, … SNAIV… Trai…  1.33  2.90  2.22  6.31  10.7     1     1 0.800
report(ets_models)
## # A tibble: 2 × 11
##   State    Industry  .model  sigma2 log_lik   AIC  AICc   BIC   MSE  AMSE    MAE
##   <chr>    <chr>     <chr>    <dbl>   <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>  <dbl>
## 1 Tasmania Cafes, r… ets_m… 0.00314  -1016. 2065. 2067. 2130.  1.38  1.99 0.0406
## 2 Tasmania Cafes, r… ets_m… 0.00328  -1020. 2077. 2079. 2146.  1.39  2.01 0.0405
report(snaive_model)
## Series: Turnover 
## Model: SNAIVE 
## 
## sigma^2: 6.6608
fc_ets |> accuracy(myseries)
## # A tibble: 2 × 12
##   .model    State Industry .type    ME  RMSE   MAE   MPE  MAPE  MASE RMSSE  ACF1
##   <chr>     <chr> <chr>    <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 ets_mam   Tasm… Cafes, … Test   2.38  3.99  3.23  4.20  6.31  1.46  1.37 0.807
## 2 ets_mam_… Tasm… Cafes, … Test   6.73  8.60  7.18 12.5  13.6   3.24  2.96 0.910
fc_snaive |> accuracy(myseries)
## # A tibble: 1 × 12
##   .model    State Industry .type    ME  RMSE   MAE   MPE  MAPE  MASE RMSSE  ACF1
##   <chr>     <chr> <chr>    <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 SNAIVE(T… Tasm… Cafes, … Test   7.12  9.13  7.58  13.2  14.4  3.42  3.15 0.863
#Plot forecasts
combined_fc <- bind_rows(fc_ets, fc_snaive)

combined_fc %>%
  autoplot(myseries, level = NULL) +
  labs(y="Turnover,  $Million AUD", title="Forecasts from ETS and SNAIVE Models") + 
  scale_color_manual(values = c("red", "purple", "blue"), labels = c("M,A,M", "M,Ad,M", "SNAIVE"))

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?

STL Decomposition + ETS(A, A, N): the model has a training accuracy of RMSE=0.0699, MAE=0.0515, AIC=190.10, and a test set accuracy of RMSE=6.6863, MAE=5.5723, MAPE=11.22%.

Previous best model ETS(M, A, M): the model has a training accuracy of RMSE=1.18, MAE=0.823, AIC=2065.0, and a test set accuracy of RMSE=3.99, MAE=3.23, MAPE=6.31%.

The STL decomposition + ETS(A, A, N) approach did not outperform the previous best model (ETS(M, A, M)) on the test set. The new model has a higher RMSE on the test set (6.6863) than the previous model (3.99). The new model’s MAE and MAPE (5.5723 and 11.22%, respectively) are also lower than those of the previous model (3.23 and 6.31%). By decomposing the series and removing seasonality before fitting the ETS model, we may have missed important seasonal patterns that the original ETS(M, A, M) model could detect. This may explain why the new approach performed poorly on the test set. While using an STL decomposition followed by an ETS model is an acceptable method for dealing with seasonality and non-constant variance, it did not outperform the simpler ETS model with multiplicative seasonality in this case. The retail data is likely to contain seasonal variations that are critical for accurate forecasting, which the ETS(M, A, M) model can handle more effectively.

#Box-Cox transformation
lambda <- BoxCox.lambda(myseries_train$Turnover)
myseries_train_transformed <- myseries_train %>%
  mutate(Turnover_transformed = BoxCox(Turnover, lambda))

#STL decomposition
stl_fit <- myseries_train_transformed %>%
  model(stl_bx = STL(Turnover_transformed ~ season(window = "periodic")))

#Seasonally adjusted data
seasonally_df <- components(stl_fit)

#ETS model, seasonally adjusted data
ets_stl <- seasonally_df %>%
  model(ets_adj = ETS(season_adjust ~ error("A") + trend("A") + season("N")))

# Forecast for the test period
fc_stl <- ets_stl %>% 
  select(ets_adj) %>% 
  forecast(h = nrow(anti_join(myseries, myseries_train)))

# Inverse Box-Cox transformation of the forecasts
fc_stl_transformed <- fc_stl %>%
  mutate(.mean = InvBoxCox(.mean, lambda))

#Accuracy
ets_stl %>% select(ets_adj) %>%  accuracy()
## # A tibble: 1 × 10
##   .model  .type           ME   RMSE    MAE      MPE  MAPE  MASE RMSSE    ACF1
##   <chr>   <chr>        <dbl>  <dbl>  <dbl>    <dbl> <dbl> <dbl> <dbl>   <dbl>
## 1 ets_adj Training 0.0000290 0.0699 0.0515 -0.00444  1.60 0.347 0.375 0.00403
report(ets_stl)
## Series: season_adjust 
## Model: ETS(A,A,N) 
##   Smoothing parameters:
##     alpha = 0.7458769 
##     beta  = 0.0001000022 
## 
##   Initial states:
##      l[0]        b[0]
##  1.852001 0.007525362
## 
##   sigma^2:  0.0049
## 
##      AIC     AICc      BIC 
## 190.1021 190.2791 209.3198
myseries_test <- anti_join(myseries, myseries_train)
fc_stl_transformed$.mean |> accuracy(myseries_test$Turnover)
##                 ME     RMSE      MAE       MPE     MAPE
## Test set -5.134939 6.686284 5.572324 -10.37392 11.22215