This assignment covers Exercises 9.1, 9.2, 9.3, 9.5, 9.6, 9.7, and 9.8 from Forecasting: Principles and Practice (3rd edition).

Exercise 9.1: ACF Plots of White Noise

Exercise 9.1(a): Differences Among the Plots

All three plots are consistent with white noise. At nonzero lags, the sample autocorrelations fluctuate around zero without a systematic pattern. As the sample size increases from 36 to 360 to 1,000, the sample autocorrelations generally become smaller because they are estimated more precisely.

An occasional spike outside the 95% bounds does not necessarily contradict white noise. Approximately 5% of sample autocorrelations may fall outside these bounds by chance, even when the underlying series is white noise.

Exercise 9.1(b): Critical Values and Sample Autocorrelations

The approximate 95% bounds for white-noise autocorrelations are:

\[ \pm \frac{1.96}{\sqrt{T}}, \]

where \(T\) is the number of observations. The bounds are approximately:

  • For 36 observations: ±0.327.
  • For 360 observations: ±0.103.
  • For 1,000 observations: ±0.062.

Larger samples produce narrower bounds because there is less sampling uncertainty.

Although the theoretical autocorrelations of white noise are zero at every nonzero lag, sample autocorrelations are not exactly zero. Each random sample produces different estimates because of sampling variation.

Exercise 9.2: Amazon Stock Prices

amazon_stock <- gafa_stock |>
  filter(Symbol == "AMZN") |>
  arrange(Date) |>
  mutate(TradingDay = row_number()) |>
  update_tsibble(index = TradingDay, regular = TRUE)

amazon_stock |>
  gg_tsdisplay(Close, plot_type = "partial") +
  labs(
    title = "Amazon Daily Closing Stock Prices",
    x = "Trading day",
    y = "Closing price (US dollars)"
  )

The time plot shows a strong overall upward movement in Amazon’s closing price, followed by a sharp decline near the end. Prices do not fluctuate around a stable mean, indicating that the series is non-stationary.

The ACF has large positive auto correlations that remain close to one and decay slowly across the displayed lags. This persistence is characteristic of a non-stationary price series.

The PACF shows a dominant spike at lag 1, close to one, with much smaller values at subsequent lags. This is consistent with unit-root-like behavior, although the PACF alone cannot establish non-stationary.

Together, these plots suggest taking a first difference to remove the changing level. The differences series represents the dollar change between consecutive trading days. Its stationarity should then be checked before fitting an ARIMA model.

Exercise 9.3: Transformations and Differencing

Exercise 9.3(a): Turkish GDP

turkey_gdp <- global_economy |>
  filter(Code == "TUR")

turkey_gdp |>
  autoplot(GDP) +
  labs(
    title = "Turkey: Annual GDP",
    x = "Year",
    y = "GDP (current US dollars)"
  )

turkey_gdp |>
  gg_tsdisplay(log(GDP), plot_type = "partial") +
  labs(
    title = "Turkey: Log GDP",
    y = "Log GDP"
  )

### Exercise 9.3(a): First Difference of Log GDP

turkey_gdp |>
  gg_tsdisplay(
    difference(log(GDP)),
    plot_type = "partial"
  ) +
  labs(
    title = "Turkey: First Difference of Log GDP",
    y = "Change in log GDP"
  )
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_point()`).

turkey_gdp |>
  features(log(GDP), unitroot_ndiffs)
## # A tibble: 1 × 2
##   Country ndiffs
##   <fct>    <int>
## 1 Turkey       1

A logarithmic transformation, equivalent to a Box–Cox transformation with λ = 0, reduces the increasing scale of fluctuations in Turkish GDP. However, log GDP still has an upward trend and a slowly decaying ACF, indicating non-stationarity.

After one ordinary difference, the logged series fluctuates around a relatively stable level without a sustained trend. Its ACF no longer shows the strong persistence seen before differencing. The unitroot_ndiffs() function also recommends one difference.

Therefore, a suitable choice is a Box–Cox transformation with λ = 0 followed by first differencing:

\[ w_t = \log(GDP_t) - \log(GDP_{t-1}). \] No seasonal difference is needed because these are annual observations. The transformed and differences series appears approximately stationary, although some large changes remain.

Exercise 9.3(b): Accommodation Takings in Tasmania

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

tasmania_takings |>
  autoplot(Takings) +
  labs(
    title = "Tasmania: Quarterly Accommodation Takings",
    x = "Year",
    y = "Takings (millions of Australian dollars)"
  )

tasmania_lambda <- tasmania_takings |>
  features(Takings, features = guerrero) |>
  pull(lambda_guerrero)

tasmania_lambda
## [1] 0.001819643
tasmania_takings |>
  autoplot(box_cox(Takings, tasmania_lambda)) +
  labs(
    title = "Tasmania: Box–Cox Transformed Accommodation Takings",
    x = "Year",
    y = "Transformed takings"
  )

Exercise 9.3(b): Seasonal Differencing

tasmania_diff <- tasmania_takings |>
  mutate(
    Transformed = box_cox(Takings, tasmania_lambda),
    SeasonalDifference = difference(Transformed, lag = 4)
  )

tasmania_diff |>
  gg_tsdisplay(SeasonalDifference, plot_type = "partial") +
  labs(
    title = "Tasmania: Seasonally Differenced Transformed Takings",
    y = "Seasonal difference"
  )
## Warning: Removed 4 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 4 rows containing missing values or values outside the scale range
## (`geom_point()`).

tasmania_diff |>
  features(SeasonalDifference, unitroot_ndiffs)
## # A tibble: 1 × 2
##   State    ndiffs
##   <chr>     <int>
## 1 Tasmania      0

The original accommodation takings show an upward trend and seasonal fluctuations that increase with the series level. Guerrero’s method selects a Box–Cox parameter of approximately λ = 0.00182, which is very close to a logarithmic transformation. This makes the seasonal fluctuations more consistent in size.

A seasonal difference at lag 4 removes the repeating quarterly pattern and much of the changing level. The resulting ACF declines relatively quickly rather than remaining close to one. Some short-term autocorrelation remains, but stationary does not require the series to be white noise.

After seasonal difference, unitroot_ndiffs() returns 0, suggesting that no additional ordinary difference is needed. Therefore, a reasonable choice is the estimated Box–Cox transformation followed by one seasonal difference: d = 0 and D = 1, with seasonal period 4.

Writing the Box–Cox transformed takings as z_t, the resulting series is:

\[ w_t = z_t - z_{t-4}. \] ### Exercise 9.3(c): Monthly Souvenir Sales

souvenirs |>
  autoplot(Sales) +
  labs(
    title = "Monthly Souvenir Sales",
    x = "Year",
    y = "Sales (Australian dollars)"
  )

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

souvenir_lambda
## [1] 0.002118221
souvenirs |>
  autoplot(box_cox(Sales, souvenir_lambda)) +
  labs(
    title = "Box–Cox Transformed Souvenir Sales",
    x = "Year",
    y = "Transformed sales"
  )

### Exercise 9.3(c): Seasonal Differencing

souvenir_diff <- souvenirs |>
  mutate(
    Transformed = box_cox(Sales, souvenir_lambda),
    SeasonalDifference = difference(Transformed, lag = 12)
  )

souvenir_diff |>
  gg_tsdisplay(SeasonalDifference, plot_type = "partial") +
  labs(
    title = "Souvenir Sales: Seasonally Differenced Box–Cox Series",
    y = "Seasonal difference"
  )
## Warning: Removed 12 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 12 rows containing missing values or values outside the scale range
## (`geom_point()`).

souvenir_diff |>
  features(SeasonalDifference, unitroot_ndiffs)
## # A tibble: 1 × 1
##   ndiffs
##    <int>
## 1      0

Souvenir sales show an upward trend and annual seasonal peaks whose size increases with the level of sales. Guerrero’s method selects a Box–Cox parameter of approximately λ = 0.00212, which is close to a logarithmic transformation. The transformation makes the seasonal variation more consistent in size.

After seasonal differencing at lag 12, the strong repeating annual pattern and sustained upward trend are substantially reduced. Some changes in level and short-term autocorrelation remain, so stationarity is only approximate.

The unitroot_ndiffs() function returns 0 for the seasonally differenced series. Based on this result and the plots, I select one seasonal difference and no additional ordinary difference: d = 0 and D = 1, with seasonal period 12.

Writing the Box–Cox transformed sales as z_t, the resulting series is:

\[ w_t = z_t - z_{t-12}. \]

Exercise 9.5: Differencing the Retail Series

retail_data <- aus_retail |>
  filter(
    State == "New South Wales",
    Industry == "Takeaway food services"
  )

retail_lambda <- retail_data |>
  features(Turnover, features = guerrero) |>
  pull(lambda_guerrero)

retail_lambda
## [1] 0.002144737
retail_data |>
  autoplot(Turnover) +
  labs(
    title = "NSW Takeaway Food Services: Monthly Turnover",
    x = "Year",
    y = "Turnover (millions of Australian dollars)"
  )

retail_data |>
  autoplot(box_cox(Turnover, retail_lambda)) +
  labs(
    title = "NSW Takeaway Food Services: Box–Cox Transformation",
    x = "Year",
    y = "Transformed turnover"
  )

### Exercise 9.5: Check Seasonal Differencing

retail_diff <- retail_data |>
  mutate(
    Transformed = box_cox(Turnover, retail_lambda),
    SeasonalDifference = difference(Transformed, lag = 12)
  )

retail_diff |>
  gg_tsdisplay(SeasonalDifference, plot_type = "partial") +
  labs(
    title = "NSW Takeaway Turnover: Seasonal Difference",
    y = "Seasonal difference"
  )
## Warning: Removed 12 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 12 rows containing missing values or values outside the scale range
## (`geom_point()`).

retail_diff |>
  features(SeasonalDifference, unitroot_ndiffs)
## # A tibble: 1 × 3
##   State           Industry               ndiffs
##   <chr>           <chr>                   <int>
## 1 New South Wales Takeaway food services      0

Guerrero’s method selects a Box–Cox parameter of approximately λ = 0.00214, close to a logarithmic transformation. This reduces the increase in fluctuation size as turnover grows.

After a seasonal difference at lag 12, the series fluctuates around a relatively stable level without the original sustained upward trend. The ACF shows strong short-term dependence but declines and crosses zero rather than remaining close to one.

The unitroot_ndiffs() function returns 0 for the seasonally differenced series, recommending no additional ordinary differencing. Therefore, I select d = 0 and D = 1, with seasonal period 12, after the Box–Cox transformation.

Writing transformed turnover as z_t, the resulting series is:

\[ w_t = z_t - z_{t-12}. \]

The series appears approximately stationary but is not white noise. Its remaining autocorrelation can be represented by AR and MA terms rather than automatically applying further differences.

Exercise 9.6: Simulating ARIMA Models

Exercise 9.6(a): Simulate an AR(1) Model

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
)

The simulated model is:

\[ y_t = 0.6y_{t-1} + \varepsilon_t, \qquad \varepsilon_t \sim N(0,1), \]

with initial value \(y_1 = 0\).

Exercise 9.6(b): Compare AR(1) Time Plots

phi_values <- c(-0.6, 0, 0.6, 0.9)
ar1_series <- vector("list", length(phi_values))

for (j in seq_along(phi_values)) {
  y <- numeric(100)

  for (i in 2:100) {
    y[i] <- phi_values[j] * y[i - 1] + e[i]
  }

  ar1_series[[j]] <- tibble(
    idx = seq_len(100),
    y = y,
    Phi = paste("phi =", phi_values[j])
  )
}

ar1_comparison <- bind_rows(ar1_series) |>
  as_tsibble(key = Phi, index = idx)

ar1_comparison |>
  autoplot(y) +
  facet_wrap(vars(Phi), ncol = 2) +
  labs(
    title = "AR(1) Simulations with Different Coefficients",
    x = "Time",
    y = "Simulated value"
  ) +
  guides(colour = "none")

With φ = −0.6, the series frequently alternates above and below zero because consecutive observations are negatively correlated.

With φ = 0, each observation after the initial value equals its random error, giving white noise with no dependence on previous observations.

With φ = 0.6, nearby observations tend to move together, producing longer runs above or below zero.

With φ = 0.9, shocks persist longer, creating more sustained movements and larger excursions. Although this can resemble a trend over short periods, the coefficient satisfies |φ| < 1, so the process is mean-reverting rather than a random walk.

All simulations start at zero and use the same random errors, allowing the effect of changing φ to be compared directly.

Exercise 9.6(c): Simulate an MA(1) Model

set.seed(625)

ma_errors <- rnorm(101)
ma_values <- numeric(100)

for (i in 1:100) {
  ma_values[i] <- ma_errors[i + 1] + 0.6 * ma_errors[i]
}

ma1_sim <- tsibble(
  idx = seq_len(100),
  y = ma_values,
  index = idx
)

The simulated model is:

\[ y_t = \varepsilon_t + 0.6\varepsilon_{t-1}, \qquad \varepsilon_t \sim N(0,1). \]

The extra random error supplies the lagged error needed for the first observation.

Exercise 9.6(d): Compare MA(1) Time Plots

theta_values <- c(-0.6, 0, 0.6, 0.9)
ma1_series <- vector("list", length(theta_values))

for (j in seq_along(theta_values)) {
  ma_values <- numeric(100)

  for (i in 1:100) {
    ma_values[i] <- ma_errors[i + 1] +
      theta_values[j] * ma_errors[i]
  }

  ma1_series[[j]] <- tibble(
    idx = seq_len(100),
    y = ma_values,
    Theta = paste("theta =", theta_values[j])
  )
}

ma1_comparison <- bind_rows(ma1_series) |>
  as_tsibble(key = Theta, index = idx)

ma1_comparison |>
  autoplot(y) +
  facet_wrap(vars(Theta), ncol = 2) +
  labs(
    title = "MA(1) Simulations with Different Coefficients",
    x = "Time",
    y = "Simulated value"
  ) +
  guides(colour = "none")

With θ = −0.6, consecutive observations tend to move in opposite directions, producing a more jagged, alternating pattern.

With θ = 0, the series is white noise because each observation equals its current random error.

With θ = 0.6 or 0.9, consecutive observations are positively correlated, producing more similar neighbouring values. The θ = 0.9 series also has somewhat larger fluctuations.

Unlike an AR(1) model, an MA(1) model does not propagate shocks indefinitely. Each error affects only the current and following observation, so the theoretical autocorrelation is zero beyond lag 1. Increasing θ therefore does not create the long-lasting persistence seen when an AR(1) coefficient approaches one.

The simulations use the same random errors, so differences between the plots reflect changes in θ.

Exercise 9.6(e): Simulate an ARMA(1,1) Model

set.seed(626)

arma_errors <- rnorm(100)
arma_values <- numeric(100)

for (i in 2:100) {
  arma_values[i] <- 0.6 * arma_values[i - 1] +
    arma_errors[i] + 0.6 * arma_errors[i - 1]
}

arma11_sim <- tsibble(
  idx = seq_len(100),
  y = arma_values,
  index = idx
)

The simulated model is:

\[ y_t = 0.6y_{t-1} + \varepsilon_t + 0.6\varepsilon_{t-1}, \qquad \varepsilon_t \sim N(0,1), \]

with initial value \(y_1 = 0\).

Exercise 9.6(f): Simulate a Non-Stationary AR(2) Model

set.seed(627)

ar2_errors <- rnorm(100)
ar2_values <- numeric(100)

for (i in 3:100) {
  ar2_values[i] <- -0.8 * ar2_values[i - 1] +
    0.3 * ar2_values[i - 2] + ar2_errors[i]
}

ar2_sim <- tsibble(
  idx = seq_len(100),
  y = ar2_values,
  index = idx
)

The simulated model is:

\[ y_t = -0.8y_{t-1} + 0.3y_{t-2} + \varepsilon_t, \qquad \varepsilon_t \sim N(0,1), \]

with initial values \(y_1 = y_2 = 0\). These coefficients produce a non-stationary process.

Exercise 9.6(g): Compare ARMA(1,1) and AR(2)

simulation_comparison <- bind_rows(
  arma11_sim |>
    as_tibble() |>
    mutate(Model = "ARMA(1,1)"),
  ar2_sim |>
    as_tibble() |>
    mutate(Model = "Non-stationary AR(2)")
) |>
  as_tsibble(key = Model, index = idx)

simulation_comparison |>
  autoplot(y) +
  facet_wrap(vars(Model), ncol = 1, scales = "free_y") +
  labs(
    title = "Comparison of ARMA(1,1) and AR(2) Simulations",
    x = "Time",
    y = "Simulated value"
  ) +
  guides(colour = "none")

The ARMA(1,1) simulation fluctuates around zero within a relatively stable range. Its positive AR coefficient creates persistence, but shocks gradually decay because |φ₁| = 0.6 < 1.

In contrast, the AR(2) simulation alternates between positive and negative values with rapidly increasing magnitude. The coefficients φ₁ = −0.8 and φ₂ = 0.3 produce an explosive, non-stationary process rather than fluctuations with stable variance.

The panels use different vertical scales. The ARMA(1,1) values remain within a few units of zero, whereas the AR(2) values reach thousands in magnitude. This highlights the difference between stable, mean-reverting behaviour and explosive behavior.

Exercise 9.7: Australian Air Passengers

Exercise 9.7(a): Select an ARIMA Model

air_fit <- aus_airpassengers |>
  model(
    Auto = ARIMA(Passengers)
  )

report(air_fit)
## 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

Exercise 9.7(a): Residual Diagnostics

air_fit |>
  gg_tsresiduals()

air_fit |>
  augment() |>
  features(.innov, ljung_box, lag = 10, dof = 1) |>
  knitr::kable(digits = 4)
.model lb_stat lb_pvalue
Auto 6.6964 0.6687

The selected model is ARIMA(0,2,1), with an estimated MA coefficient of −0.8963 and no constant. It uses second-order difference and one moving-average term.

The residuals fluctuate around zero, and the displayed ACF values lie within the approximate 95% bounds. The Ljung–Box test gives Q* = 6.6964 and p = 0.6687 at lag 10, accounting for the estimated MA coefficient. We therefore fail to reject the null hypothesis of no residual autocorrelation through these lags.

The residuals are reasonably consistent with being uncorrelated, but several large residuals and the histogram’s tails suggest that normality and constant variance may be imperfect. Consequently, prediction intervals should be interpreted cautiously.

Exercise 9.7(a): Ten-Year Forecasts

air_fc <- air_fit |>
  forecast(h = 10)

air_fc |>
  autoplot(aus_airpassengers) +
  labs(
    title = "Australian Air Passengers: Ten-Year ARIMA Forecasts",
    x = "Year",
    y = "Passengers (millions)"
  )

The ARIMA(0,2,1) model forecasts continued growth in Australian air passenger numbers over the next ten years. The point forecasts follow a straight upward path, extending the model’s estimated recent growth rate. The 80% and 95% prediction intervals widen with the forecast horizon, indicating increasing uncertainty further into the future.

Exercise 9.7(b): Backshift Notation

The fitted ARIMA(0,2,1) model can be written as

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

where \(y_t\) represents annual passenger numbers and \(B\) is the backshift operator, so \(By_t=y_{t-1}\).

Expanding gives

\[ y_t-2y_{t-1}+y_{t-2} = \varepsilon_t-0.8963\varepsilon_{t-1}. \]

The model has no constant, and the estimated innovation variance is \(\hat{\sigma}^2=4.308\).

Exercise 9.7(c): ARIMA(0,1,0) with Drift

air_drift_fit <- aus_airpassengers |>
  model(
    Auto = ARIMA(Passengers),
    Drift = ARIMA(Passengers ~ 1 + pdq(0, 1, 0))
  )

air_drift_fit |>
  select(Drift) |>
  report()
## 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
air_drift_fc <- air_drift_fit |>
  forecast(h = 10)

air_drift_fc |>
  autoplot(aus_airpassengers, level = NULL) +
  labs(
    title = "Australian Air Passengers: Automatic ARIMA versus Drift",
    x = "Year",
    y = "Passengers (millions)"
  )

The ARIMA(0,1,0) model with drift estimates an average annual increase of 1.4191 million passengers. Its forecasts form a straight line, increasing by this amount each year.

Both models forecast continued growth, but the automatic ARIMA(0,2,1) model produces a steeper forecast path. The drift model uses the average historical annual change, whereas the automatic model allows the growth rate to evolve over time and extends its estimated recent rate.

Exercise 9.7(d): ARIMA(2,1,2) With and Without Drift

air_models <- aus_airpassengers |>
  model(
    Auto = ARIMA(Passengers),
    Drift = ARIMA(Passengers ~ 1 + pdq(0, 1, 0)),
    ARIMA212_Drift = ARIMA(Passengers ~ 1 + pdq(2, 1, 2)),
    ARIMA212_NoDrift = ARIMA(
  Passengers ~ 0 + pdq(2, 1, 2),
  method = "ML")
  )

air_models |>
  select(ARIMA212_Drift) |>
  report()
## 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
air_models |>
  select(ARIMA212_NoDrift) |>
  report()
## Series: Passengers 
## Model: ARIMA(2,1,2) 
## 
## Coefficients:
##          ar1      ar2      ma1     ma2
##       1.8081  -0.8111  -1.9642  1.0000
## s.e.  0.1121   0.1126   0.1091  0.1103
## 
## sigma^2 estimated as 3.9:  log likelihood=-97.19
## AIC=204.38   AICc=205.88   BIC=213.52
air_models |>
  forecast(h = 10) |>
  autoplot(aus_airpassengers, level = NULL) +
  labs(
    title = "Australian Air Passengers: Comparing ARIMA Forecasts",
    x = "Year",
    y = "Passengers (millions)"
  )

The ARIMA(2,1,2) model with drift produces forecasts close to those from the ARIMA(0,1,0) drift model and below the automatic ARIMA(0,2,1) forecasts.

Removing the constant changes the estimated AR and MA coefficients substantially. Over the ten-year horizon, the no-drift ARIMA(2,1,2) model produces the highest forecasts. Therefore, removing drift does not necessarily produce flat short-term forecasts: the fitted ARMA dynamics can sustain growth for some time.

The no-drift model has an AICc of 205.88, slightly below 206.61 for the model with drift. However, the with-drift model has an MA coefficient of 1.0000, indicating a fit at the invertibility boundary and warranting caution.

Exercise 9.7(e): ARIMA(0,2,1) with a Constant

air_constant_fit <- aus_airpassengers |>
  model(
    Constant = 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(air_constant_fit)
## 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
air_constant_fit |>
  forecast(h = 10) |>
  autoplot(aus_airpassengers) +
  labs(
    title = "Australian Air Passengers: ARIMA(0,2,1) with a Constant",
    x = "Year",
    y = "Passengers (millions)"
  )

Including a constant in the ARIMA(0,2,1) model produces a quadratic forecast trend. The estimated constant is 0.0571, so the forecasts imply an increasing annual growth rate rather than a constant annual increase.

The forecast path curves upward, while the prediction intervals widen with the horizon. R warns that quadratic or higher-order polynomial trends are generally discouraged because extrapolating this acceleration can produce unrealistic long-term forecasts.

The estimated MA coefficient is -1.0000, placing the fit at the invertibility boundary and providing another reason for caution.

Exercise 9.8(a): US GDP and Box–Cox Transformation

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

us_gdp |>
  autoplot(GDP) +
  labs(
    title = "United States: Annual GDP",
    x = "Year",
    y = "GDP (current US dollars)"
  )

us_lambda <- us_gdp |>
  features(GDP, features = guerrero) |>
  pull(lambda_guerrero)

us_lambda
## [1] 0.2819443
us_gdp |>
  autoplot(box_cox(GDP, us_lambda)) +
  labs(
    title = "United States: Box–Cox Transformed GDP",
    x = "Year",
    y = "Transformed GDP"
  )

US GDP shows a strong upward trend, with larger absolute changes at higher GDP levels. Guerrero’s method selects a Box–Cox parameter of approximately 0.2819.

The transformation compresses the larger observations and makes the trend more nearly linear. It is therefore a reasonable choice for reducing level-dependent variation before fitting an ARIMA model. However, the transformed series still trends upward, so transformation alone does not make it stationary; differencing must also be considered.

Exercise 9.8(b): Automatic ARIMA Model

us_auto_fit <- us_gdp |>
  model(
    Auto = ARIMA(box_cox(GDP, us_lambda))
  )

report(us_auto_fit)
## Series: GDP 
## Model: ARIMA(1,1,0) w/ drift 
## Transformation: box_cox(GDP, us_lambda) 
## 
## 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

Automatic selection chooses an ARIMA(1,1,0) model with drift for Box–Cox transformed US GDP. The model uses first-order differencing and one autoregressive term, with an estimated AR coefficient of 0.4586. Its AICc is 657.10.

Exercise 9.8(c): Compare Alternative ARIMA Models

us_models <- us_gdp |>
  model(
    Auto = ARIMA(box_cox(GDP, us_lambda)),
    ARIMA010 = ARIMA(
      box_cox(GDP, us_lambda) ~ 1 + pdq(0, 1, 0)
    ),
    ARIMA011 = ARIMA(
      box_cox(GDP, us_lambda) ~ 1 + pdq(0, 1, 1)
    ),
    ARIMA110 = ARIMA(
      box_cox(GDP, us_lambda) ~ 1 + pdq(1, 1, 0)
    ),
    ARIMA111 = ARIMA(
      box_cox(GDP, us_lambda) ~ 1 + pdq(1, 1, 1)
    ),
    ARIMA210 = ARIMA(
      box_cox(GDP, us_lambda) ~ 1 + pdq(2, 1, 0)
    )
  )

us_models |>
  glance() |>
  select(.model, AIC, AICc, BIC) |>
  arrange(AICc) |>
  knitr::kable(digits = 4)
.model AIC AICc BIC
Auto 656.6476 657.1005 662.7768
ARIMA110 656.6476 657.1005 662.7768
ARIMA011 658.7317 659.1845 664.8608
ARIMA111 658.6443 659.4135 666.8165
ARIMA210 658.6449 659.4141 666.8171
ARIMA010 667.5345 667.7568 671.6206

The ARIMA(1,1,0) model with drift has the lowest AICc, 657.1005, among the models compared. This agrees with automatic selection. The ARIMA(0,1,1) model is relatively competitive, but adding further AR or MA terms does not improve AICc. I therefore retain ARIMA(1,1,0) with drift, subject to checking its residuals.

Exercise 9.8(d): Residual Diagnostics

us_best_fit <- us_models |>
  select(Country, ARIMA110)

us_best_fit |>
  gg_tsresiduals()

us_best_fit |>
  augment() |>
  features(.innov, ljung_box, lag = 10, dof = 1) |>
  knitr::kable(digits = 4)
Country .model lb_stat lb_pvalue
United States ARIMA110 3.8137 0.9232

The residual ACF values lie within the approximate 95% bounds. The Ljung–Box test gives a statistic of 3.8137 and a p-value of 0.9232, so we fail to reject the null hypothesis of no residual autocorrelation through lag 10.

The residuals fluctuate around zero, although several large negative observations and a left-tailed histogram suggest departures from normality. Overall, the diagnostics support retaining ARIMA(1,1,0) with drift, while interpreting its prediction intervals cautiously.

Exercise 9.8(e): US GDP Forecasts

us_best_fit <- us_models |>
  select(Country, ARIMA110)

us_arima_fc <- us_best_fit |>
  forecast(h = 10)

us_arima_fc |>
  autoplot(us_gdp) +
  labs(
    title = "United States GDP: Ten-Year ARIMA Forecasts",
    x = "Year",
    y = "GDP (current US dollars)"
  )

The forecasts show continued growth in US GDP over the next ten years, broadly extending the historical upward pattern. On the original GDP scale, the forecast path curves upward, and the 80% and 95% prediction intervals widen with the forecast horizon.

These forecasts appear plausible as an extrapolation of historical nominal GDP growth. However, they assume that the historical dynamics continue and cannot anticipate recessions or structural changes. The intervals also depend on model assumptions, so they should not be interpreted as capturing every source of economic uncertainty.

Exercise 9.8(f): Compare ARIMA and ETS Forecasts

us_final_models <- us_gdp |>
  model(
    ARIMA = ARIMA(
      box_cox(GDP, us_lambda) ~ 1 + pdq(1, 1, 0)
    ),
    ETS = ETS(GDP)
  )

us_final_models |>
  select(Country, ETS) |>
  report()
## 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
us_final_fc <- us_final_models |>
  forecast(h = 10)

us_final_fc |>
  autoplot(us_gdp, level = NULL) +
  labs(
    title = "United States GDP: ARIMA versus ETS Forecasts",
    x = "Year",
    y = "GDP (current US dollars)"
  )

Automatic ETS selection chooses ETS(M,A,N): multiplicative errors, an additive trend, and no seasonality, fitted to un-transformed GDP.

Both models forecast continued growth. The ETS point forecasts follow a straight upward path, while the Box–Cox ARIMA forecasts curve upward and become higher over the ten-year horizon.

The difference reflects their treatment of trend and transformation. The plot alone does not establish which model forecasts more accurately; that would require out-of-sample evaluation.