DATA 624 – ARIMA Model

Exercises 9.1, 9.2, 9.3, 9.5, 9.6, 9.7 and 9.8 from Forecasting: Principles and Practice (3rd ed.), Chapter 9.

---

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?

Yes, they are all white noise. In every plot almost every autocorrelation sits inside the blue dashed lines, which is exactly what white noise looks like (Section 2.9 of the book says we expect about 95% of the spikes to be inside the bounds).

The difference between the figures is just the sample size. With 36 numbers the bounds are wide and the spikes jump around a lot, with more numbers like 1,000 or so the bounds become tight and the spikes are all close to zero. Same type of data, we just have more of it in the later plots.

(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 (blue lines) are \(\pm 1.96/\sqrt{T}\), where \(T\) is the number of observations. More data means a smaller standard error, so the lines just moves closer to zero

# critical values for the three sample sizes in the figure
T_vals <- c(36, 360, 1000)
tibble(T = T_vals, critical_value = round(1.96 / sqrt(T_vals), 3))

## # A tibble: 3 × 2
##       T critical_value
##   <dbl>          <dbl>
## 1    36          0.327
## 2   360          0.103
## 3  1000          0.062

The autocorrelations themselves are different because they are from different samples of a population, so the estimated autocorrelations are random too. The true autocorrelation of white noise is 0 at every lag, but with only 36 points the estimate can wander off pretty far from 0. With 1,000 points the estimates are much closer to the truth.

---

Exercise 9.2

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

# Pull out just Amazon from the gafa_stock data
amazon <- gafa_stock |>
  filter(Symbol == "AMZN")

# Time plot + ACF + PACF in one shot however the textbook used gg_tsdisplay for this
amazon |>
  gg_tsdisplay(Close, plot_type = "partial") +
  labs(title = "Amazon daily closing price")

Each plot gives its own clue that the series is not stationary as we see that the time plot shows a clear upward trend and a level that keeps changing, while a stationary series should be roughly flat with a constant mean and this clearly shows that is not . The ACF provides further evidence as the autocorrelations are all very large, positive, and decay super slowly instead of dropping to zero quickly, and the book points out that a slow decaying ACF with a big \(r_1\) is the classic sign of non stationary data. The PACF has one huge spike at lag 1 (close to 1) and basically nothing after that, meaning today’s price is basically explained by yesterday’s price, so the data behave like a random walk, which is non stationary. For these reasons the series should be differenced, after taking first differences we are modelling the daily changes in price, which should look much more like white noise.

# Quick check: the differenced series looks stationary
amazon |>
  gg_tsdisplay(difference(Close), plot_type = "partial") +
  labs(title = "Amazon daily change in closing price")

---

Exercise 9.3

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

For each series I follow the same recipe from the book:

1. pick a Box-Cox \(\lambda\) with the Guerrero method (guerrero feature),
2. check whether a seasonal difference is needed with unitroot_nsdiffs(),
3. check how many first differences are needed with unitroot_ndiffs() (KPSS test),
4. plot the result to make sure it looks stationary.

(a) Turkish GDP from global_economy

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

# First the raw data
turkey |>
  autoplot(GDP) +
  labs(title = "Turkish GDP (raw)", y = "USD")

# Box-Cox lambda chosen by Guerrero 
lambda_tr <- turkey |>
  features(GDP, features = guerrero) |>
  pull(lambda_guerrero)

# How many first differences does the KPSS procedure suggest? there's annual data, so no seasonal 
nd_tr <- turkey |>
  features(box_cox(GDP, lambda_tr), unitroot_ndiffs) |>
  pull(ndiffs)

tibble(lambda = round(lambda_tr, 3), ndiffs = nd_tr)

## # A tibble: 1 × 2
##   lambda ndiffs
##    <dbl>  <int>
## 1  0.157      1

# Transformed and  differenced series 
turkey |>
  mutate(y = box_cox(GDP, lambda_tr),
         y = if (nd_tr > 0) difference(y, differences = nd_tr) else y) |>
  gg_tsdisplay(y, plot_type = "partial") +
  labs(title = "Turkish GDP: Box-Cox then differenced")

Answer: Guerrero picked \(\lambda \approx\) 0.157 which pulled in the exploding variance of the raw series. So after the transformation the KPSS procedure says we need 1 first differences. The differenced series bounces around a constant mean with no trend that’s obvous, and the ACF drops off quickly, so it looks stationary.

(b) Accommodation takings in the state of Tasmania from aus_accommodation

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

tas |>
  autoplot(Takings) +
  labs(title = "Tasmania accommodation takings (raw)", y = "$ million")

lambda_tas <- tas |>
  features(Takings, features = guerrero) |>
  pull(lambda_guerrero)

# Quarterly data, checks seasonal differencing first at m = 4, then the first differencing
nsd_tas <- tas |>
  features(box_cox(Takings, lambda_tas), unitroot_nsdiffs) |>
  pull(nsdiffs)

nd_tas <- tas |>
  features(difference(box_cox(Takings, lambda_tas), lag = 4), unitroot_ndiffs) |>
  pull(ndiffs)

tibble(lambda = round(lambda_tas, 3), nsdiffs = nsd_tas, ndiffs = nd_tas)

## # A tibble: 1 × 3
##   lambda nsdiffs ndiffs
##    <dbl>   <int>  <int>
## 1  0.002       1      0

# Apply exactly the differences the tests recommend
tas |>
  mutate(y = box_cox(Takings, lambda_tas),
         y = if (nsd_tas > 0) difference(y, lag = 4) else y,
         y = if (nd_tas  > 0) difference(y) else y) |>
  gg_tsdisplay(y, plot_type = "partial") +
  labs(title = "Tasmania takings: Box-Cox + differencing")

Answer: The raw series has strong quarterly seasonality and the seasonal swings get bigger as the level grows, so in this case Box-Cox transformation is needed. Guerrero gives \(\lambda \approx\) 0.002 which is close to lag. The tests then call for 1 seasonal differences and 0 first differences on top of that, after that the series has a stable mean and variance.

(c) Monthly sales from souvenirs

souvenirs |>
  autoplot(Sales) +
  labs(title = "Souvenir shop sales (raw)")

lambda_sv <- souvenirs |>
  features(Sales, features = guerrero) |>
  pull(lambda_guerrero)

# Monthly data so seasonal lag is 12
nsd_sv <- souvenirs |>
  features(box_cox(Sales, lambda_sv), unitroot_nsdiffs) |>
  pull(nsdiffs)

nd_sv <- souvenirs |>
  features(difference(box_cox(Sales, lambda_sv), lag = 12), unitroot_ndiffs) |>
  pull(ndiffs)

tibble(lambda = round(lambda_sv, 3), nsdiffs = nsd_sv, ndiffs = nd_sv)

## # A tibble: 1 × 3
##   lambda nsdiffs ndiffs
##    <dbl>   <int>  <int>
## 1  0.002       1      0

souvenirs |>
  mutate(y = box_cox(Sales, lambda_sv),
         y = if (nsd_sv > 0) difference(y, lag = 12) else y,
         y = if (nd_sv  > 0) difference(y) else y) |>
  gg_tsdisplay(y, plot_type = "partial") +
  labs(title = "Souvenir sales -  Box-Cox and  differencing")

Answer: Souvenir sales have a huge spike in December every year and the spikes grow with the trend, so clearly the variance is y not constant. Guerrero gives \(\lambda \approx\) 0.002, which is basically a log transformation. The tests then suggest 1 seasonal differences and 0 first differences. The plot after transforming and differencing looks stationary.

---

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(12345678)
myseries <- aus_retail |>
  filter(`Series ID` == sample(aus_retail$`Series ID`, 1))

myseries |>
  autoplot(Turnover) +
  labs(title = paste(unique(myseries$State), "-", unique(myseries$Industry)))

# Box-Cox
lambda_rt <- myseries |>
  features(Turnover, features = guerrero) |>
  pull(lambda_guerrero)

# seasonal differencing  and since it's monthly lag is 12
nsd_rt <- myseries |>
  features(box_cox(Turnover, lambda_rt), unitroot_nsdiffs) |>
  pull(nsdiffs)

#  first differencing on the seasonally differenced data
nd_rt <- myseries |>
  features(difference(box_cox(Turnover, lambda_rt), lag = 12), unitroot_ndiffs) |>
  pull(ndiffs)

tibble(lambda = round(lambda_rt, 3), nsdiffs = nsd_rt, ndiffs = nd_rt)

## # A tibble: 1 × 3
##   lambda nsdiffs ndiffs
##    <dbl>   <int>  <int>
## 1  0.083       1      0

# plot the transformed and  differenced series and run a KPSS test on it
retail_diff <- myseries |>
  mutate(y = box_cox(Turnover, lambda_rt),
         y = if (nsd_rt > 0) difference(y, lag = 12) else y,
         y = if (nd_rt  > 0) difference(y) else y)

retail_diff |>
  gg_tsdisplay(y, plot_type = "partial") +
  labs(title = "Retail turnover: Box-Cox + differencing")

# KPSS - thenull hypothesis = stationary, so a p-value of 0.1 which is the max it reports is what we want
retail_diff |>
  features(y, unitroot_kpss)

## # A tibble: 1 × 4
##   State              Industry                              kpss_stat kpss_pvalue
##   <chr>              <chr>                                     <dbl>       <dbl>
## 1 Northern Territory Clothing, footwear and personal acce…     0.407      0.0741

Answer: The retail series is seasonal with a trend, and the seasonal swings grow with the level, so I first transform (Guerrero \(\lambda \approx\) 0.083). The unit root features recommend 1 seasonal differences followed by 0 first differences. The KPSS p-value on the final series is 0.0740653, so we don’t reject stationarity and no further differencing is needed.

---

Exercise 9.6

Simulate and plot some data from simple ARIMA models.

set.seed(624)  # so the plots are reproducible

(a) AR(1) with \(\phi_1 = 0.6\) and \(\sigma^2 = 1\) (code from the book)

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) Time plot. How does the plot change as you change \(\phi_1\)?

sim |>
  autoplot(y) +
  labs(title = "AR(1) with phi = 0.6")

# Wrap the book's provided loop in a function so we can try several phi values
sim_ar1 <- function(phi, n = 100) {
  y <- numeric(n)
  e <- rnorm(n)
  for (i in 2:n) y[i] <- phi * y[i - 1] + e[i]
  tsibble(idx = seq_len(n), y = y, index = idx)
}

# Same random seed for each so the only thing changing is phi
phis <- c(-0.8, 0, 0.6, 0.95)
ar1_all <- bind_rows(lapply(phis, function(p) {
  set.seed(1)
  sim_ar1(p) |> as_tibble() |> mutate(phi = p)   # as_tibble so bind_rows doesn't complain
})) |>
  as_tsibble(index = idx, key = phi)

ar1_all |>
  autoplot(y) +
  facet_wrap(~phi, ncol = 1, labeller = label_both) +
  labs(title = "AR(1) for different values of phi") +
  theme(legend.position = "none")

Answer: \(\phi_1\) controls how much “memory” the series has.

• \(\phi_1 = 0\): it’s just white noise.
• \(\phi_1 = 0.6\): closeby values are correlated, so you get short runs above and under zero.
• \(\phi_1 = 0.95\): very smooth, long swings wandering, it’s getting close to a random walk (\(\phi_1 = 1\)).
• \(\phi_1 = -0.8\): the series flips sign almost every step and constantly swings around zero which like the book said, negative \(\phi_1\) makes the series oscillate around the mean.

(c) MA(1) with \(\theta_1 = 0.6\) and \(\sigma^2 = 1\)

# MA(1): y_t = e_t + theta * e_{t-1}
sim_ma1 <- function(theta, n = 100) {
  y <- numeric(n)
  e <- rnorm(n)
  for (i in 2:n) y[i] <- e[i] + theta * e[i - 1]
  tsibble(idx = seq_len(n), y = y, index = idx)
}

set.seed(624)
ma1 <- sim_ma1(0.6)

(d) Time plot. How does the plot change as you change \(\theta_1\)?

ma1 |>
  autoplot(y) +
  labs(title = "MA(1) with theta = 0.6")

thetas <- c(-0.8, 0, 0.6, 0.95)
ma1_all <- bind_rows(lapply(thetas, function(th) {
  set.seed(1)
  sim_ma1(th) |> as_tibble() |> mutate(theta = th)
})) |>
  as_tsibble(index = idx, key = theta)

ma1_all |>
  autoplot(y) +
  facet_wrap(~theta, ncol = 1, labeller = label_both) +
  labs(title = "MA(1) for different values of theta") +
  theme(legend.position = "none")

Answer: Changing \(\theta_1\) changes the looks much less than changing \(\phi_1\) did. An MA(1) only remembers one error back, so the series always snaps back to zero quickly no matter what \(\theta_1\) is. Positive \(\theta_1\) makes neighbouring points move together a little which create a slightly smoother line, negative \(\theta_1\) makes them tend to move in opposite directions which is sort of a jagged line. \(\theta_1 = 0\) is white noise.

(e) ARMA(1,1) with \(\phi_1 = 0.6\), \(\theta_1 = 0.6\) and \(\sigma^2 = 1\)

# ARMA(1,1): y_t = phi * y_{t-1} + e_t + theta * e_{t-1}
set.seed(624)
n <- 100
y <- numeric(n)
e <- rnorm(n)
for (i in 2:n) {
  y[i] <- 0.6 * y[i - 1] + e[i] + 0.6 * e[i - 1]
}
arma11 <- tsibble(idx = seq_len(n), y = y, index = idx)

(f) AR(2) with \(\phi_1 = -0.8\), \(\phi_2 = 0.3\) and \(\sigma^2 = 1\) (non-stationary)

# AR(2): y_t = phi1 * y_{t-1} + phi2 * y_{t-2} + e_t   (loop starts at 3 because we need two lags)
set.seed(624)
y <- numeric(n)
e <- rnorm(n)
for (i in 3:n) {
  y[i] <- -0.8 * y[i - 1] + 0.3 * y[i - 2] + e[i]
}
ar2 <- tsibble(idx = seq_len(n), y = y, index = idx)

(g) Graph the latter two series and compare them

arma11 |> autoplot(y) + labs(title = "ARMA(1,1): phi = 0.6, theta = 0.6")

ar2    |> autoplot(y) + labs(title = "AR(2): phi1 = -0.8, phi2 = 0.3 (non-stationary)")

Answer: The ARMA(1,1) series is stationary, it wanders a bit but keeps coming back to zero and its variance stays about the same across the entire plot. The AR(2) series does the opposite. It flips back and forth every step because \(\phi_1\) is negative and the swings keep growing until the values are huge. This happens because the parameters break the AR(2) stationarity condition from Section 9.3 (\(\phi_2 - \phi_1 < 1\), but here \(0.3 - (-0.8) = 1.1\)), so the series blows up instead of settling down. fable would never estimate a model like this since it has constraints specifically for this.

---

Exercise 9.7

Consider aus_airpassengers, the total number of passengers (in millions) from Australian air carriers.

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

(a) Use ARIMA() to find an appropriate model. Check the residuals. Plot forecasts for the next 10 periods.

fit_auto <- aus_airpassengers |>
  model(ARIMA(Passengers))

report(fit_auto)

## Series: Passengers 
## Model: ARIMA(0,2,1) 
## 
## Coefficients:
##           ma1
##       -0.8963
## s.e.   0.0594
## 
## sigma^2 estimated as 4.308:  log likelihood=-97.02
## AIC=198.04   AICc=198.32   BIC=201.65

# Residual diagnostics -  time plot, ACF, histogram
fit_auto |> gg_tsresiduals()

# Ljung-Box with dof = number of AR and  MA parameters 
augment(fit_auto) |>
  features(.innov, ljung_box, lag = 10, dof = 1)

## # A tibble: 1 × 3
##   .model            lb_stat lb_pvalue
##   <chr>               <dbl>     <dbl>
## 1 ARIMA(Passengers)    6.70     0.669

fit_auto |>
  forecast(h = 10) |>
  autoplot(aus_airpassengers) +
  labs(title = "ARIMA() forecasts, next 10 years", y = "Millions")

Answer: `ARIMA() picked an ARIMA(0,2,1) model. That means the data were differenced twice. There are no AR terms and there’s also one MA term with \(\theta_1 \approx\) r round(theta1, 3). The residuals look like white noise as the ACF has no spikes outside the bounds, the histogram looks pretty normal, and the Ljung-Box p-value is large, so the model has captured the pattern in the data. Since the model uses two differences with no constant, the book says the long run forecasts will follow a straight line, and that’s exactly what the plot shows, the forecasts keeps that recent upward trend going, with prediction intervals that widen very quickly.

(b) Write the model in terms of the backshift operator

Using Equation 9.2 from the book with \(p = 0\), \(d = 2\), \(q = 1\) and no constant

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

and plugging in the estimate

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

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

# "~ 1" forces the constant (drift) in; pdq() fixes the orders
fit_010 <- aus_airpassengers |>
  model(arima010_drift = ARIMA(Passengers ~ 1 + pdq(0, 1, 0)))

report(fit_010)

## Series: Passengers 
## Model: ARIMA(0,1,0) w/ drift 
## 
## Coefficients:
##       constant
##         1.4191
## s.e.    0.3014
## 
## sigma^2 estimated as 4.271:  log likelihood=-98.16
## AIC=200.31   AICc=200.59   BIC=203.97

fit_010 |>
  forecast(h = 10) |>
  autoplot(aus_airpassengers) +
  labs(title = "ARIMA(0,1,0) with drift", y = "Millions")

Answer: This is the random walk with drift from Table 9.1. The point forecasts are also a straight line, but the slope is the average yearly change over the entire series which is the constant as where the ARIMA(0,2,1) slope is driven by the most recent changes and since growth has been faster lately, the ARIMA(0,2,1) line is a little bit steeper. The prediction intervals are narrower here because \(d = 1\) instead of \(d = 2\) and the book does note that intervals grow faster for higher \(d\).

(d) ARIMA(2,1,2) with drift. Compare to (a) and (c). Remove the constant and see what happens.

fit_212 <- aus_airpassengers |>
  model(arima212_drift = 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
##       -0.5518  -0.7327  0.5895  1.0000    3.2216
## s.e.   0.1684   0.1224  0.0916  0.1008    0.7224
## 
## sigma^2 estimated as 4.031:  log likelihood=-96.23
## AIC=204.46   AICc=206.61   BIC=215.43

# All three models on one plot for comparison... here level = NULL would hide the intervals so the lines are easy to see
aus_airpassengers |>
  model(
    auto_021      = ARIMA(Passengers),
    arima010_drift = ARIMA(Passengers ~ 1 + pdq(0, 1, 0)),
    arima212_drift = ARIMA(Passengers ~ 1 + pdq(2, 1, 2))
  ) |>
  forecast(h = 10) |>
  autoplot(aus_airpassengers, level = NULL) +
  labs(title = "Comparing the three models", y = "Millions")

# Now dropping the constant ("~ 0"). Warnings is  left on here because the warning is the point here.
fit_212_nc <- aus_airpassengers |>
  model(arima212_noconst = ARIMA(Passengers ~ 0 + pdq(2, 1, 2)))

## Warning: 1 error encountered for arima212_noconst
## [1] non-stationary AR part from CSS

fit_212_nc

## # A mable: 1 x 1
##   arima212_noconst
##            <model>
## 1     <NULL model>

Answer: The ARIMA(2,1,2) with drift gives forecasts very close to the other two models seen earlier, just with a slight curve from the AR and MA terms. All three point in the same direction, so the extra parameters don’t add much here. When the constant is removed, fable gives a warning and returns a NULL model instead of a fit, without the constant the model has no way of capturing the upward trend, so the AR part ends up non stationary and fable refuses to return a model that breaks the stationarity conditions, so there’s no model and nothing to forecast.

(e) ARIMA(0,2,1) with a constant. What happens?

fit_021c <- aus_airpassengers |>
  model(arima021_const = ARIMA(Passengers ~ 1 + pdq(0, 2, 1)))

## Warning: Model specification induces a quadratic or higher order polynomial trend.  This
## is generally discouraged, consider removing the constant or reducing the number
## of differences.

report(fit_021c)

## Series: Passengers 
## Model: ARIMA(0,2,1) w/ poly 
## 
## Coefficients:
##           ma1  constant
##       -1.0000    0.0571
## s.e.   0.0585    0.0213
## 
## sigma^2 estimated as 3.855:  log likelihood=-95.1
## AIC=196.21   AICc=196.79   BIC=201.63

fit_021c |>
  forecast(h = 10) |>
  autoplot(aus_airpassengers) +
  labs(title = "ARIMA(0,2,1) with a constant", y = "Millions")

Answer: fable gives a warning that this model creates a quadratic trend. That matches the book when there is a constant and \(d = 2\), the long term forecasts follow a curve instead of a straight line. In the plot the forecast bends upward and as time passes it speeds up. The book said to avoid this because it assumes growth will keep accelerating forever, which is indeed risky. That’s why ARIMA() never adds a constant when \(d = 2\) unless you force it.

---

Exercise 9.8

For the United States GDP series (from global_economy):

us_gdp <- global_economy |>
  filter(Code == "USA")

us_gdp |>
  autoplot(GDP) +
  labs(title = "US GDP (raw)", y = "USD")

(a) If necessary, find a suitable Box-Cox transformation

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 = paste0("US GDP, Box-Cox with lambda = ", round(lambda_us, 3)))

Answer: The raw series curves upward with a sort of exponential growth-ish kind of, so a transformation helps. Guerrero picks \(\lambda \approx\) 0.282. After the transformation the growth is much closer to a straight line, which is easier for a differenced ARIMA model to handle.

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

fit_us_auto <- us_gdp |>
  model(auto = ARIMA(box_cox(GDP, lambda_us)))

report(fit_us_auto)

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

Answer: ARIMA() picks an ARIMA(1,1,0) with drift, one difference to remove the trend, one AR term, and a constant so the forecasts keep growing.

(c) Try some other plausible models by experimenting with the orders

# Have a look at the differenced data to get ideas for p and q
us_gdp |>
  gg_tsdisplay(difference(box_cox(GDP, lambda_us)), plot_type = "partial") +
  labs(title = "Differenced Box-Cox US GDP")

# All candidates use d = 1 so their AICc values are comparable
fits_us <- us_gdp |>
  model(
    auto     = ARIMA(box_cox(GDP, lambda_us)),
    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(fits_us) |>
  arrange(AICc) |>
  select(.model:BIC)

## # A tibble: 5 × 6
##   .model   sigma2 log_lik   AIC  AICc   BIC
##   <chr>     <dbl>   <dbl> <dbl> <dbl> <dbl>
## 1 auto      5479.   -325.  657.  657.  663.
## 2 arima011  5689.   -326.  659.  659.  665.
## 3 arima111  5580.   -325.  659.  659.  667.
## 4 arima210  5580.   -325.  659.  659.  667.
## 5 arima212  5734.   -325.  662.  664.  674.

Answer: The ACF and PACF of the differenced data shows a big spike at lag 1 and not much after, which is why simple models like AR(1) or MA(1) are the candidates that fit naturally. I tried ARIMA(0,1,1), (1,1,1), (2,1,0) and (2,1,2), all with drift. The AICc values are close together and the lowest AICc belongs to auto.

(d) Choose the best model and check the residual diagnostics

# Gets out the best model by AICc
fit_best <- fits_us |> select(Country, all_of(best_us))

fit_best |> gg_tsresiduals()

# dof = number of AR and MA coefficients basically everything except for the constants
K <- tidy(fit_best) |> filter(term != "constant") |> nrow()

augment(fit_best) |>
  features(.innov, ljung_box, lag = 10, dof = K)

## # A tibble: 1 × 4
##   Country       .model lb_stat lb_pvalue
##   <fct>         <chr>    <dbl>     <dbl>
## 1 United States auto      3.81     0.923

Answer: I’m going with auto (lowest AICc). The residual plot shows no pattern, the ACF spikes are all inside the bounds, the histogram is basically normal, and the Ljung-Box p-value is large, so the residuals behave like white noise and the model is good to go.

(e) Produce forecasts. Do they look reasonable?

fit_best |>
  forecast(h = 10) |>
  autoplot(us_gdp) +
  labs(title = paste("US GDP forecasts from", best_us), y = "USD")

Answer: Yes. fable automatically transforms the forecasts back to the original dollar scale, and they continue the steady upward trend of the last few decades at about the same rate. The prediction intervals widen as we go further out, which makes sense for a differenced model (\(d \ge 1\) means the intervals keep growing).

(f) Compare the results with ETS() (no transformation)

fit_ets <- us_gdp |>
  model(ets = ETS(GDP))

report(fit_ets)

## Series: GDP 
## Model: ETS(M,A,N) 
##   Smoothing parameters:
##     alpha = 0.9990876 
##     beta  = 0.5011949 
## 
##   Initial states:
##          l[0]        b[0]
##  448093333334 64917355687
## 
##   sigma^2:  7e-04
## 
##      AIC     AICc      BIC 
## 3190.787 3191.941 3201.089

# Forecasts side by side
bind_rows(
  fit_best |> forecast(h = 10),
  fit_ets  |> forecast(h = 10)
) |>
  autoplot(us_gdp, level = NULL) +
  labs(title = "ARIMA vs ETS forecasts for US GDP", y = "USD")

# AICc can't be directly compared across ARIMA and ETS as they are diffirent model classes 
# so we will pretend the last 10 year dont exist and we will fit bot models up to 2007 then ask each one to forecast
train <- us_gdp |> filter(Year <= 2007)

train |>
  model(
    arima = ARIMA(box_cox(GDP, lambda_us) ~ 1 + pdq(1, 1, 0)),
    ets   = ETS(GDP)
  ) |>
  forecast(h = 10) |>
  fabletools::accuracy(us_gdp) |>
  select(.model, RMSE, MAE, MAPE, MASE)

## # A tibble: 2 × 5
##   .model    RMSE     MAE  MAPE  MASE
##   <chr>    <dbl>   <dbl> <dbl> <dbl>
## 1 arima  1.98e12 1.87e12 11.1   6.32
## 2 ets    1.72e12 1.66e12  9.93  5.61

Answer: Both models forecast GDP going up. They agree for the first few years and then slowly start growing apart and that’s because ETS works directly on the dollar values while ARIMA works on the Box-Cox scale and converts back. The book says AICc can’t be used to compare ETS and ARIMA since they are different types of models, so I tested them another way which was by fitting both using data up to 2007,then forecasted the next 10 years, and then checked which forecasts were closer to the real numbers. The model with the lower RMSE and MASE in the table above is the one I’d pick.