library(fpp3)
library(urca)   # for unitroot_kpss / ndiffs

9.1 White noise ACFs (Figure 9.32)

a. The three ACFs indicate white noise: no autocorrelation exceeds the bounds in any panel. The left plot (T = 36) has the largest spikes, the biggest near lag 11 at about 0.3, just inside the ±0.33 bound. As T increases, both the spikes and the bounds shrink toward zero.

b. The bounds are \(\pm 1.96/\sqrt{T}\), so they shrink as \(T\) grows: about \(\pm 0.33\) for \(T=36\), \(\pm 0.10\) for \(T=360\), and \(\pm 0.06\) for \(T=1000\). The sample autocorrelations differ between the plots even though each series is white noise because each is an estimate from a different random sample; sampling variability is larger for small \(T\), so the \(T=36\) spikes are bigger in absolute size.

1.96 / sqrt(c(36, 360, 1000))
## [1] 0.32666667 0.10330107 0.06198064

9.2 Amazon closing prices

amzn <- gafa_stock |>
  filter(Symbol == "AMZN") |>
  mutate(day = row_number()) |>
  update_tsibble(index = day, regular = TRUE)   # trading days are irregular

amzn |> gg_tsdisplay(Close, plot_type = "partial")

amzn |> features(Close, list(unitroot_kpss, unitroot_ndiffs))
## # A tibble: 1 × 4
##   Symbol kpss_stat kpss_pvalue ndiffs
##   <chr>      <dbl>       <dbl>  <int>
## 1 AMZN        14.0        0.01      1
amzn |> gg_tsdisplay(difference(Close), plot_type = "partial")

amzn |> features(difference(Close), ljung_box, lag = 10)
## # A tibble: 1 × 3
##   Symbol lb_stat lb_pvalue
##   <chr>    <dbl>     <dbl>
## 1 AMZN      35.6 0.0000981

So the series should be first-differenced. After differencing the mean is stable, though the variance grows over time and the Ljung–Box test on the differences is still significant (p ≈ 0.0001), so the differences are stationary in mean but not pure white noise.

9.3 Box-Cox and differencing

a. Turkish GDP

turkey <- global_economy |> filter(Country == "Turkey")
lam_tk <- turkey |> features(GDP, guerrero) |> pull(lambda_guerrero)
lam_tk
## [1] 0.1572187
turkey |> features(box_cox(GDP, lam_tk), unitroot_ndiffs)
## # A tibble: 1 × 2
##   Country ndiffs
##   <fct>    <int>
## 1 Turkey       1
turkey |> mutate(d = difference(box_cox(GDP, lam_tk))) |>
    features(d, list(unitroot_kpss, ~ ljung_box(.x, lag = 10)))
## # A tibble: 1 × 5
##   Country kpss_stat kpss_pvalue lb_stat lb_pvalue
##   <fct>       <dbl>       <dbl>   <dbl>     <dbl>
## 1 Turkey     0.0889         0.1    5.84     0.829
turkey |> gg_tsdisplay(difference(box_cox(GDP, lam_tk)), plot_type = "partial")

Guerrero λ ≈ 0.16 (close to a log). One first difference is needed (ndiffs = 1). After differencing, KPSS p ≥ 0.10, so stationarity is not rejected, and Ljung–Box (lag 10) p ≈ 0.83, so the differences also look like white noise.

b. Tasmania accommodation takings

tas <- aus_accommodation |> filter(State == "Tasmania")
tas |> features(Takings, guerrero)          # ~0.002 -> use log
## # A tibble: 1 × 2
##   State    lambda_guerrero
##   <chr>              <dbl>
## 1 Tasmania         0.00182
tas |> features(log(Takings), unitroot_nsdiffs)
## # A tibble: 1 × 2
##   State    nsdiffs
##   <chr>      <int>
## 1 Tasmania       1
tas |> mutate(s = difference(log(Takings), 4)) |>
  features(s, list(unitroot_ndiffs, unitroot_kpss))
## # A tibble: 1 × 4
##   State    ndiffs kpss_stat kpss_pvalue
##   <chr>     <int>     <dbl>       <dbl>
## 1 Tasmania      0     0.199         0.1
tas |> gg_tsdisplay(difference(log(Takings), 4), plot_type = "partial")

λ ≈ 0.002, so a log transformation. The data are quarterly with strong seasonality: one seasonal difference (lag 4) is needed (nsdiffs = 1). After it, ndiffs = 0 and KPSS p ≥ 0.10, so no further differencing is needed. The seasonally differenced series is stationary but not white noise: the ACF has large spikes at lags 1–3 that die out quickly. Transformation: log; differencing: D = 1, d = 0.

c. Souvenir sales

souvenirs |> features(Sales, guerrero)      # ~0.002 -> use log
## # A tibble: 1 × 1
##   lambda_guerrero
##             <dbl>
## 1         0.00212
souvenirs |> features(log(Sales), unitroot_nsdiffs)
## # A tibble: 1 × 1
##   nsdiffs
##     <int>
## 1       1
souvenirs |> mutate(s = difference(log(Sales), 12)) |>
  features(s, list(unitroot_ndiffs, unitroot_kpss))
## # A tibble: 1 × 3
##   ndiffs kpss_stat kpss_pvalue
##    <int>     <dbl>       <dbl>
## 1      0     0.239         0.1
souvenirs |> gg_tsdisplay(difference(log(Sales), 12), plot_type = "partial")
## 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()`).

λ ≈ 0.002, so log. Monthly data: one seasonal difference (lag 12) is needed (nsdiffs = 1). After it, ndiffs = 0 and KPSS p ≥ 0.10, so stationarity is not rejected. The seasonally differenced series still wanders, and the ACF decays fairly slowly over the first 10 lags, so an additional first difference would also be defensible. Transformation: log; differencing: D = 1, d = 0.

9.5 My retail series

set.seed(12345678)
myseries <- aus_retail |>
  filter(`Series ID` == sample(aus_retail$`Series ID`, 1))
distinct(as_tibble(myseries), State, Industry)
## # A tibble: 1 × 2
##   State              Industry                                           
##   <chr>              <chr>                                              
## 1 Northern Territory Clothing, footwear and personal accessory retailing
lam_r <- myseries |> features(Turnover, guerrero) |> pull(lambda_guerrero)
lam_r
## [1] 0.08303631
myseries |> features(box_cox(Turnover, lam_r), unitroot_nsdiffs)
## # A tibble: 1 × 3
##   State              Industry                                            nsdiffs
##   <chr>              <chr>                                                 <int>
## 1 Northern Territory Clothing, footwear and personal accessory retailing       1
myseries |> mutate(s = difference(box_cox(Turnover, lam_r), 12)) |>
  features(s, list(unitroot_ndiffs, unitroot_kpss))
## # A tibble: 1 × 5
##   State              Industry                       ndiffs kpss_stat kpss_pvalue
##   <chr>              <chr>                           <int>     <dbl>       <dbl>
## 1 Northern Territory Clothing, footwear and person…      0     0.407      0.0741
myseries |> gg_tsdisplay(difference(box_cox(Turnover, lam_r), 12), plot_type = "partial")
## 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()`).

myseries |> mutate(s = difference(difference(box_cox(Turnover, lam_r), 12))) |>
  features(s, 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.0165         0.1

With this seed the series is Northern Territory clothing, footwear and personal accessory retailing. Guerrero λ ≈ 0.08 (near log). One seasonal difference is needed (nsdiffs = 1). After it, ndiffs = 0, but KPSS p ≈ 0.074 is borderline, and the plot still wanders with a slowly decaying ACF. An extra first difference brings the KPSS statistic down to 0.017 (p ≥ 0.10). I’d use D = 1, d = 1.

9.6 Simulated ARIMA models

set.seed(624)
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)
}
sim <- sim_ar1(0.6)

b. AR(1) with different \(\phi_1\):

bind_rows(
  as_tibble(sim_ar1(0.6)) |> mutate(phi = "0.6"),
  as_tibble(sim_ar1(0)) |> mutate(phi = "0"),
  as_tibble(sim_ar1(0.95)) |> mutate(phi = "0.95"),
  as_tibble(sim_ar1(-0.6)) |> mutate(phi = "-0.6")
) |>
  as_tibble() |>
  as_tsibble(index = idx, key = phi) |>
  autoplot(y) + facet_wrap(~phi, ncol = 1, scales = "free_y") +
  theme(legend.position = "none")

As \(\phi_1 \to 0\) the series looks like white noise. As \(\phi_1 \to 1\) it becomes smoother and wanders away from the mean for longer stretches (more like a random walk). Negative \(\phi_1\) makes the series flip sign from one period to the next, producing a jagged, oscillating plot.

c. MA(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)
}

d.

bind_rows(
  as_tibble(sim_ma1(0.6)) |> mutate(theta = "0.6"),
  as_tibble(sim_ma1(0)) |> mutate(theta = "0"),
  as_tibble(sim_ma1(0.95)) |> mutate(theta = "0.95"),
  as_tibble(sim_ma1(-0.6)) |> mutate(theta = "-0.6")
) |>
  as_tibble() |>
  as_tsibble(index = idx, key = theta) |>
  autoplot(y) + facet_wrap(~theta, ncol = 1, scales = "free_y") +
  theme(legend.position = "none")

Changes are subtle because an MA(1) only links neighboring observations. Larger positive \(\theta_1\) makes the series slightly smoother and increases its variance (\(1+\theta_1^2\)); negative \(\theta_1\) makes it choppier. It doesn’t show the long excursions an AR(1) with \(\phi_1\) near 1 does.

e. ARMA(1,1) with \(\phi_1 = 0.6\), \(\theta_1 = 0.6\):

y <- numeric(100); e <- rnorm(100)
for (i in 2:100) y[i] <- 0.6 * y[i - 1] + e[i] + 0.6 * e[i - 1]
arma11 <- tsibble(idx = seq_len(100), y = y, index = idx)

f. AR(2) with \(\phi_1 = -0.8\), \(\phi_2 = 0.3\):

y <- numeric(100); e <- rnorm(100)
for (i in 3:100) y[i] <- -0.8 * y[i - 1] + 0.3 * y[i - 2] + e[i]
ar2 <- tsibble(idx = seq_len(100), y = y, index = idx)
Mod(polyroot(c(1, 0.8, -0.3)))   # a root inside the unit circle -> non-stationary
## [1] 0.9274433 3.5941100

g.

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

The ARMA(1,1) is stationary: it wanders around zero with some persistence but stays within a stable band. The AR(2) violates the stationarity condition (\(\phi_2 - \phi_1 = 1.1 > 1\); one characteristic root has modulus ≈ 0.93 < 1), so it oscillates with exploding amplitude; the swings alternate sign because \(\phi_1\) is negative and grow into the hundreds or thousands by the end of 100 observations.

9.7 Australian air passengers

air <- aus_airpassengers |> filter(Year <= 2011)
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
fit_air |> gg_tsresiduals()

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

a. ARIMA() selects ARIMA(0,2,1) with \(\theta_1 = -0.876\) and \(\sigma^2 = 4.67\). The residual ACF shows no significant spikes and the Ljung–Box test gives p = 0.74, so the residuals look like white noise. Forecasts rise roughly linearly from about 66.0 million (2012) to 89.6 million (2021).

b. \((1 - B)^2 y_t = (1 - 0.876B)\,\varepsilon_t\)

c.

fit_c <- air |> model(
  auto  = ARIMA(Passengers ~ pdq(0, 2, 1)),
  drift = ARIMA(Passengers ~ 1 + pdq(0, 1, 0))
)
fit_c |> select(drift) |> report()
## 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_c |> forecast(h = 10) |> autoplot(air, level = NULL)

The random walk with drift (drift = 1.37) forecasts a straight line with the average slope over 1970–2011, reaching about 77.0 million in 2021. That is flatter than the ARIMA(0,2,1), which extrapolates the steeper recent slope.

d.

fit_d <- air |> model(
  a212_drift = ARIMA(Passengers ~ 1 + pdq(2, 1, 2)),
  a212_none  = ARIMA(Passengers ~ 0 + pdq(2, 1, 2))
)
## Warning: 1 error encountered for a212_none
## [1] non-stationary AR part from CSS
fit_d
## # A mable: 1 x 2
##                a212_drift    a212_none
##                   <model>      <model>
## 1 <ARIMA(2,1,2) w/ drift> <NULL model>
fit_d |> select(a212_drift) |> report()
## Series: Passengers 
## Model: ARIMA(2,1,2) w/ drift 
## 
## Coefficients:
##          ar1      ar2      ma1     ma2  constant
##       1.4694  -0.5103  -1.5736  0.6780    0.0650
## s.e.  0.3780   0.3558   0.3081  0.2688    0.0294
## 
## sigma^2 estimated as 4.748:  log likelihood=-87.74
## AIC=187.47   AICc=189.94   BIC=197.75
fit_d |> select(a212_drift) |> forecast(h = 10) |> autoplot(air)

# Removing the constant: the default estimation fails, so try ML
fit_d2 <- air |> model(a212_none = ARIMA(Passengers ~ 0 + pdq(2, 1, 2), method = "ML"))
report(fit_d2)
## Series: Passengers 
## Model: ARIMA(2,1,2) 
## 
## Coefficients:
##          ar1      ar2      ma1     ma2
##       1.8128  -0.8154  -1.9642  0.9999
## s.e.  0.1321   0.1320   0.1295  0.1310
## 
## sigma^2 estimated as 4.203:  log likelihood=-88.19
## AIC=186.37   AICc=188.09   BIC=194.94

ARIMA(2,1,2) with drift forecasts about 65.2 → 88.2 million, very close to part (a) and steeper than part (c). With the constant removed, the default fit returns a NULL model (error: “non-stationary AR part from CSS”). Forcing method = "ML" produces a fit, but \(\hat\phi_1 + \hat\phi_2 \approx 0.997\) (an AR root essentially on the unit circle) and \(\hat\theta_2 \approx 1\) (on the invertibility boundary). Without a drift term, the model can only reproduce the trend by sneaking in a near-second difference, which is a sign the specification is wrong.

e.

fit_e <- air |> model(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_e)
## 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_e |> forecast(h = 10) |> autoplot(air)

R warns that the specification “induces a quadratic or higher order polynomial trend.” A constant with \(d = 2\) adds a quadratic trend, so forecasts curve upward (about 95.4 million by 2021 versus 89.6 in part a). The MA coefficient also hits −1.00, the invertibility boundary. This model is not advisable.

9.8 United States GDP

a.

us <- global_economy |> filter(Code == "USA")
us |> autoplot(GDP)

lam_us <- us |> features(GDP, guerrero) |> pull(lambda_guerrero)
lam_us
## [1] 0.2819443
us |> autoplot(box_cox(GDP, lam_us))

us |> features(box_cox(GDP, lam_us), unitroot_ndiffs)
## # A tibble: 1 × 2
##   Country       ndiffs
##   <fct>          <int>
## 1 United States      1

Guerrero λ ≈ 0.28; the transformed series is close to linear. One difference is needed.

b–c.

fit_us <- us |> model(
  auto = ARIMA(box_cox(GDP, lam_us)),
  full = ARIMA(box_cox(GDP, lam_us), stepwise = FALSE, approximation = FALSE),
  a011 = ARIMA(box_cox(GDP, lam_us) ~ 1 + pdq(0, 1, 1)),
  a111 = ARIMA(box_cox(GDP, lam_us) ~ 1 + pdq(1, 1, 1)),
  a210 = ARIMA(box_cox(GDP, lam_us) ~ 1 + pdq(2, 1, 0))
)
fit_us
## # A mable: 1 x 6
## # Key:     Country [1]
##   Country                          auto                    full
##   <fct>                         <model>                 <model>
## 1 United States <ARIMA(1,1,0) w/ drift> <ARIMA(1,1,0) w/ drift>
## # ℹ 3 more variables: a011 <model>, a111 <model>, a210 <model>
glance(fit_us) |> select(.model, AICc, BIC) |> arrange(AICc)
## # A tibble: 5 × 3
##   .model  AICc   BIC
##   <chr>  <dbl> <dbl>
## 1 auto    657.  663.
## 2 full    657.  663.
## 3 a011    659.  665.
## 4 a111    659.  667.
## 5 a210    659.  667.

Both stepwise and full search choose ARIMA(1,1,0) with drift (\(\hat\phi_1 = 0.459\)). Alternatives with the same \(d\): ARIMA(0,1,1), (1,1,1) and (2,1,0), all with AICc ≈ 659 vs 657.1 for (1,1,0). AICc can only be compared across models with the same differencing, so a \(d = 2\) model is not included in this comparison and all alternatives keep d = 1.

d.

fit_best <- fit_us |> select(Country, auto)
report(fit_best)
## Series: GDP 
## Model: ARIMA(1,1,0) w/ drift 
## Transformation: box_cox(GDP, lam_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
fit_best |> gg_tsresiduals()

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

Residuals are centered on zero, the ACF has no significant spikes, and Ljung–Box p = 0.92, so they are consistent with white noise.

e.

fc_us <- fit_best |> forecast(h = 10)
fc_us |> autoplot(us)

fc_us |> hilo(95) |> filter(Year %in% c(2018, 2027))
## # A tsibble: 2 x 6 [1Y]
## # Key:       Country, .model [1]
##   Country       .model  Year                 GDP   .mean
##   <fct>         <chr>  <dbl>              <dist>   <dbl>
## 1 United States auto    2018   t(N(19996, 5479)) 2.02e13
## 2 United States auto    2027 t(N(21962, 160246)) 2.82e13
## # ℹ 1 more variable: `95%` <hilo>

Forecasts continue the trend from about $20.2 trillion (2018) to $28.2 trillion (2027; 95% interval ≈ $24.7–31.9 trillion). That looks reasonable: it follows the long-run growth path, and the intervals widen with horizon. It does not anticipate recessions such as 2008–09.

f.

fit_ets <- us |> model(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
bind_rows(
  fc_us |> mutate(.model = "ARIMA(1,1,0) drift, Box-Cox"),
  fit_ets |> forecast(h = 10) |> mutate(.model = "ETS(M,A,N)")
) |>
  autoplot(us, alpha = 0.6)

augment(fit_ets) |> features(.innov, ljung_box, lag = 10)
## # A tibble: 1 × 4
##   Country       .model   lb_stat lb_pvalue
##   <fct>         <chr>      <dbl>     <dbl>
## 1 United States ETS(GDP)    3.82     0.955
fit_ets |> forecast(h = 10) |> hilo(95) |> filter(Year == 2027)
## # A tsibble: 1 x 6 [1Y]
## # Key:       Country, .model [1]
##   Country       .model    Year
##   <fct>         <chr>    <dbl>
## 1 United States ETS(GDP)  2027
## # ℹ 3 more variables: GDP <dist>, .mean <dbl>, `95%` <hilo>
# Out-of-sample comparison: train to 2007, test 2008-2017
train_fit <- us |> filter(Year <= 2007) |> model(
  arima = ARIMA(box_cox(GDP, lam_us) ~ 1 + pdq(1, 1, 0)),
  ets   = ETS(GDP)
)
train_fit |> forecast(h = 10) |> accuracy(us) |> select(.model, RMSE, MAPE)
## # A tibble: 2 × 3
##   .model    RMSE  MAPE
##   <chr>    <dbl> <dbl>
## 1 arima  1.98e12 11.1 
## 2 ets    1.72e12  9.93
# 2025 intervals, to compare with actual GDP
bind_rows(
  fc_us |> hilo(95) |> as_tibble() |> filter(Year == 2025) |> select(.model, .mean, `95%`),
  fit_ets |> forecast(h = 10) |> hilo(95) |> as_tibble() |> filter(Year == 2025) |> select(.model, .mean, `95%`)
)
## # A tibble: 2 × 3
##   .model     .mean                          `95%`
##   <chr>      <dbl>                         <hilo>
## 1 auto     2.62e13 [2.334909e+13, 2.928386e+13]95
## 2 ETS(GDP) 2.48e13 [1.556156e+13, 3.402390e+13]95

For comparison, I fit an ETS model to the untransformed data; R chose ETS(M,A,N). Its forecasts are lower than the ARIMA’s (≈ $26.1 vs $28.2 trillion in 2027), with much wider intervals. Both models have white-noise residuals, and AICc can’t be compared because the data differ. ETS did slightly better on a 2008–2017 holdout (RMSE $1.72T vs $1.98T). Since the data end in 2017, the forecasts can now be checked: actual GDP was about $30.8 trillion in 2025 (BEA), outside the ARIMA’s 95% interval but inside ETS’s. ETS is the better choice here.