Exercise 9.1

set.seed(624)
wn <- function(n) tsibble(idx = seq_len(n), y = rnorm(n), index = idx)

wn(36)   |> ACF(y) |> autoplot() + labs(title = "White noise, n = 36")

wn(360)  |> ACF(y) |> autoplot() + labs(title = "White noise, n = 360")

wn(1000) |> ACF(y) |> autoplot() + labs(title = "White noise, n = 1000")

(a) The three plots differ in how wide the dashed bounds are and in how large the spikes are. The spikes get smaller and the bounds get narrower as the sample size grows. All three indicate white noise, because nearly all of the autocorrelations fall inside the bounds. About 5% of spikes are expected to cross the bounds by chance, so the few small crossings in the \(n = 1000\) plot do not change that.

(b) The critical values are \(\pm 1.96/\sqrt{T}\), so they shrink as the series length \(T\) grows: about \(\pm 0.33\) for \(T = 36\), \(\pm 0.10\) for \(T = 360\), and \(\pm 0.06\) for \(T = 1000\). The autocorrelations differ because each series is a different random sample. The true autocorrelation is zero, but each sample estimate has random error, and that error is larger for short series.

Exercise 9.2

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

amzn |>
  gg_tsdisplay(Close, plot_type = "partial") +
  labs(title = "Amazon daily closing price")

All three plots show the series is non-stationary:

  • Time plot: a strong upward trend, so the mean changes over time.
  • ACF: the autocorrelations are large, positive and decay very slowly, which is the signature of a trended, non-stationary series.
  • PACF: a single very large spike at lag 1 (close to 1) and little after it, meaning each value is almost exactly the previous value plus noise, like a random walk.

After one difference:

amzn |>
  gg_tsdisplay(difference(Close), plot_type = "partial") +
  labs(title = "Amazon closing price, first difference")

amzn |> features(Close, unitroot_ndiffs)
## # A tibble: 1 × 2
##   Symbol ndiffs
##   <chr>   <int>
## 1 AMZN        1

The first difference removes the trend, and the KPSS-based check confirms one difference is enough. The differenced series is not perfect white noise: a few ACF spikes are just outside the bounds, and the variance grows near the end as the price level rises.

Exercise 9.3

lam <- function(data, var) {
  data |> features({{ var }}, features = guerrero) |> pull(lambda_guerrero)
}

(a) Turkish GDP

turkey <- global_economy |> filter(Country == "Turkey")
l_turkey <- lam(turkey, GDP)
l_turkey
## [1] 0.1572187
turkey |> features(box_cox(GDP, l_turkey), unitroot_ndiffs)
## # A tibble: 1 × 2
##   Country ndiffs
##   <fct>    <int>
## 1 Turkey       1
turkey |>
  gg_tsdisplay(difference(box_cox(GDP, l_turkey)), plot_type = "partial") +
  labs(title = "Turkish GDP, Box-Cox + first difference")

A Box-Cox transformation with \(\lambda \approx 0.16\) stabilizes the variance, and one first difference makes the series stationary. The data are annual, so no seasonal differencing is needed.

(b) Tasmanian accommodation takings

tas <- aus_accommodation |> filter(State == "Tasmania")
l_tas <- lam(tas, Takings)
l_tas
## [1] 0.001819643
tas |> features(box_cox(Takings, l_tas), unitroot_nsdiffs)
## # A tibble: 1 × 2
##   State    nsdiffs
##   <chr>      <int>
## 1 Tasmania       1
tas |>
  mutate(d = difference(box_cox(Takings, l_tas), 4)) |>
  features(d, unitroot_ndiffs)
## # A tibble: 1 × 2
##   State    ndiffs
##   <chr>     <int>
## 1 Tasmania      0
tas |>
  gg_tsdisplay(difference(box_cox(Takings, l_tas), 4), plot_type = "partial") +
  labs(title = "Tasmanian takings, log + seasonal difference")

\(\lambda \approx 0\), so a log transformation is suitable. The quarterly data need one seasonal difference (lag 4), and after that no further first difference is needed.

(c) Souvenir sales

l_souv <- lam(souvenirs, Sales)
l_souv
## [1] 0.002118221
souvenirs |> features(box_cox(Sales, l_souv), unitroot_nsdiffs)
## # A tibble: 1 × 1
##   nsdiffs
##     <int>
## 1       1
souvenirs |>
  mutate(d = difference(box_cox(Sales, l_souv), 12)) |>
  features(d, unitroot_ndiffs)
## # A tibble: 1 × 1
##   ndiffs
##    <int>
## 1      0
souvenirs |>
  gg_tsdisplay(difference(box_cox(Sales, l_souv), 12), plot_type = "partial") +
  labs(title = "Souvenir sales, log + seasonal difference")

Again \(\lambda \approx 0\), so a log transformation, followed by one seasonal difference (lag 12) for the monthly data. No further first difference is needed.

Exercise 9.5

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

l_retail <- lam(myseries, Turnover)
l_retail
## [1] 0.08303631
myseries |> features(box_cox(Turnover, l_retail), unitroot_nsdiffs)
## # A tibble: 1 × 3
##   State              Industry                                            nsdiffs
##   <chr>              <chr>                                                 <int>
## 1 Northern Territory Clothing, footwear and personal accessory retailing       1
myseries |>
  mutate(d = difference(box_cox(Turnover, l_retail), 12)) |>
  features(d, unitroot_ndiffs)
## # A tibble: 1 × 3
##   State              Industry                                            ndiffs
##   <chr>              <chr>                                                <int>
## 1 Northern Territory Clothing, footwear and personal accessory retailing      0
myseries |>
  gg_tsdisplay(difference(box_cox(Turnover, l_retail), 12),
               plot_type = "partial", lag_max = 36) +
  labs(title = "Retail turnover, Box-Cox + seasonal difference")

After a Box-Cox transformation (\(\lambda \approx 0.08\)), the series needs one seasonal difference (lag 12). The KPSS-based check then reports zero further differences, so the seasonally differenced series is treated as stationary. Its ACF still decays fairly slowly over the first year of lags, so adding a first difference as well would also be a defensible choice.

Exercise 9.6

sim_to_tsibble <- function(y) tsibble(idx = seq_along(y), y = y, index = idx)

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]
  sim_to_tsibble(y)
}

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]
  sim_to_tsibble(y)
}

(a) AR(1) with \(\phi_1 = 0.6\), \(\sigma^2 = 1\):

set.seed(624)
ar1 <- sim_ar1(0.6)

(b) Time plot, and how it changes with \(\phi_1\):

set.seed(624)
bind_rows(lapply(c(-0.9, 0, 0.6, 0.95), \(p) as_tibble(sim_ar1(p)) |> mutate(phi = p))) |>
  ggplot(aes(idx, y)) +
  geom_line() +
  facet_wrap(~ phi, ncol = 1, labeller = label_both, scales = "free_y") +
  labs(title = "AR(1) series for different phi_1")

With \(\phi_1 = 0\) the series is white noise. As \(\phi_1\) moves toward 1, each value depends more on the last one, so the series wanders in longer, smoother runs away from zero. A negative \(\phi_1\) makes the series flip sign from one observation to the next, so it looks jagged and oscillates rapidly around zero.

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

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

(d) Time plot, and how it changes with \(\theta_1\):

set.seed(624)
bind_rows(lapply(c(-0.9, 0, 0.6, 0.9), \(t) as_tibble(sim_ma1(t)) |> mutate(theta = t))) |>
  ggplot(aes(idx, y)) +
  geom_line() +
  facet_wrap(~ theta, ncol = 1, labeller = label_both, scales = "free_y") +
  labs(title = "MA(1) series for different theta_1")

Changing \(\theta_1\) has a much smaller effect than changing \(\phi_1\), since an MA(1) only carries one period of memory. A positive \(\theta_1\) makes neighbouring values a little more alike and slightly smoother, and a negative \(\theta_1\) makes the series a little more jagged. A larger \(|\theta_1|\) also increases the variance, \((1 + \theta_1^2)\sigma^2\).

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

set.seed(624)
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 <- sim_to_tsibble(y)

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

set.seed(624)
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 <- sim_to_tsibble(y)

(g) Compare:

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

ar2 |> autoplot(y) + labs(title = "AR(2): phi_1 = -0.8, phi_2 = 0.3")

The ARMA(1,1) series is stationary. It wanders around zero with smooth, positively correlated runs, but its level and spread stay roughly constant.

The AR(2) series is non-stationary. Its parameters break the stationarity condition \(\phi_2 - \phi_1 < 1\) (here \(0.3 + 0.8 = 1.1\)), so the oscillations grow without bound. The negative \(\phi_1\) makes it flip sign every period, and the swings grow exponentially over time.

Exercise 9.7

(a) Automatic ARIMA:

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
fit_auto |> gg_tsresiduals()

augment(fit_auto) |> features(.innov, ljung_box, lag = 10)
## # A tibble: 1 × 3
##   .model            lb_stat lb_pvalue
##   <chr>               <dbl>     <dbl>
## 1 ARIMA(Passengers)    6.70     0.754
fit_auto |> forecast(h = 10) |>
  autoplot(aus_airpassengers) +
  labs(title = "Air passengers: ARIMA(0,2,1) forecast", y = "Millions")

ARIMA() selects ARIMA(0,2,1). The residuals look like white noise: the ACF has no significant spikes and the Ljung-Box p-value is large (about 0.75).

(b) In backshift notation:

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

(c) ARIMA(0,1,0) with drift:

fit_cmp <- aus_airpassengers |> model(
  `ARIMA(0,2,1)` = ARIMA(Passengers ~ pdq(0, 2, 1)),
  `ARIMA(0,1,0) drift` = ARIMA(Passengers ~ 1 + pdq(0, 1, 0))
)

fit_cmp |> forecast(h = 10) |>
  autoplot(aus_airpassengers, level = NULL) +
  labs(title = "ARIMA(0,2,1) vs random walk with drift", y = "Millions")

The random walk with drift extends a straight line along the average growth over the whole series. The ARIMA(0,2,1) forecasts slope up more steeply, because the second difference lets the trend follow the faster growth in recent years.

(d) ARIMA(2,1,2) with drift, then without the constant:

fit_212 <- aus_airpassengers |> model(
  `ARIMA(0,2,1)` = ARIMA(Passengers ~ pdq(0, 2, 1)),
  `ARIMA(0,1,0) drift` = ARIMA(Passengers ~ 1 + pdq(0, 1, 0)),
  `ARIMA(2,1,2) drift` = ARIMA(Passengers ~ 1 + pdq(2, 1, 2))
)

fit_212 |> forecast(h = 10) |>
  autoplot(aus_airpassengers, level = NULL) +
  labs(title = "Adding ARIMA(2,1,2) with drift", y = "Millions")

aus_airpassengers |> model(ARIMA(Passengers ~ 0 + pdq(2, 1, 2)))
## # A mable: 1 x 1
##   `ARIMA(Passengers ~ 0 + pdq(2, 1, 2))`
##                                  <model>
## 1                           <NULL model>

ARIMA(2,1,2) with drift gives forecasts very close to the random walk with drift. It has the same linear long-run trend from the constant, with a little short-term movement from the AR and MA terms.

Removing the constant fails: the result is a <NULL model>, because the estimated AR part is non-stationary. Without a constant, a \(d = 1\) model has no way to produce a trend, so the optimizer pushes the AR coefficients to the non-stationary boundary to try to fit the strong upward trend, and estimation breaks down.

(e) ARIMA(0,2,1) with a constant:

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

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

fable warns that the model “induces a quadratic or higher order polynomial trend”. With \(d = 2\), a constant adds a quadratic trend, so the forecasts curve upward faster and faster. That is rarely realistic and is why fitting a constant with \(d = 2\) is discouraged.

Exercise 9.8

(a) Box-Cox transformation:

us <- global_economy |> filter(Country == "United States")
l_us <- lam(us, GDP)
l_us
## [1] 0.2819443
us |> autoplot(box_cox(GDP, l_us)) +
  labs(title = "US GDP, Box-Cox transformed")

US GDP grows exponentially, so a transformation helps. Guerrero selects \(\lambda \approx 0.28\), which makes the growth close to linear.

(b) Automatic ARIMA on the transformed data:

fit_us <- us |> model(ARIMA(box_cox(GDP, l_us)))
report(fit_us)
## Series: GDP 
## Model: ARIMA(1,1,0) w/ drift 
## Transformation: box_cox(GDP, l_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

ARIMA() selects ARIMA(1,1,0) with drift.

(c) Other plausible models:

fits_us <- us |> model(
  `ARIMA(1,1,0) drift` = ARIMA(box_cox(GDP, l_us) ~ 1 + pdq(1, 1, 0)),
  `ARIMA(0,1,1) drift` = ARIMA(box_cox(GDP, l_us) ~ 1 + pdq(0, 1, 1)),
  `ARIMA(2,1,0) drift` = ARIMA(box_cox(GDP, l_us) ~ 1 + pdq(2, 1, 0)),
  `ARIMA(1,1,1) drift` = ARIMA(box_cox(GDP, l_us) ~ 1 + pdq(1, 1, 1)),
  `ARIMA(0,2,1)` = ARIMA(box_cox(GDP, l_us) ~ 0 + pdq(0, 2, 1))
)

glance(fits_us) |> select(.model, AICc) |> arrange(AICc)
## # A tibble: 5 × 2
##   .model              AICc
##   <chr>              <dbl>
## 1 ARIMA(0,2,1)        650.
## 2 ARIMA(1,1,0) drift  657.
## 3 ARIMA(0,1,1) drift  659.
## 4 ARIMA(1,1,1) drift  659.
## 5 ARIMA(2,1,0) drift  659.

(d) Among the \(d = 1\) models, ARIMA(1,1,0) with drift has the lowest AICc, so it is the best choice. ARIMA(0,2,1) shows a lower AICc, but AICc cannot compare models with different orders of differencing, because the likelihood is computed on different data. That number does not make it the better model.

best_us <- fits_us |> select(Country, `ARIMA(1,1,0) drift`)
best_us |> gg_tsresiduals()

augment(best_us) |> features(.innov, ljung_box, lag = 10)
## # A tibble: 1 × 4
##   Country       .model             lb_stat lb_pvalue
##   <fct>         <chr>                <dbl>     <dbl>
## 1 United States ARIMA(1,1,0) drift    3.81     0.955

The residuals look like white noise: no significant ACF spikes and a Ljung-Box p-value of about 0.96. The histogram is roughly normal apart from one large negative residual in 2009, from the financial crisis.

(e) Forecasts:

best_us |> forecast(h = 10) |>
  autoplot(us) +
  labs(title = "US GDP: ARIMA(1,1,0) with drift forecast", y = "US$")

The forecasts look reasonable. They continue the steady growth of recent decades, and the prediction intervals widen with the horizon as expected.

(f) Compare with ETS on the untransformed data:

us_both <- us |> model(
  ARIMA = ARIMA(box_cox(GDP, l_us) ~ 1 + pdq(1, 1, 0)),
  ETS = ETS(GDP)
)
us_both
## # 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(us_both) |> select(.model, RMSE, MAPE, MASE)
## # A tibble: 2 × 4
##   .model          RMSE  MAPE  MASE
##   <chr>          <dbl> <dbl> <dbl>
## 1 ARIMA  147232166284.  1.52 0.257
## 2 ETS    166906260796.  1.91 0.307
us_both |> forecast(h = 10) |>
  autoplot(us, level = NULL) +
  labs(title = "US GDP: ARIMA vs ETS", y = "US$")

ETS() selects ETS(M,A,N), multiplicative errors with an additive trend. The ARIMA model fits the training data better (lower RMSE, MAPE and MASE). Its forecasts also grow a little faster, because the Box-Cox transformation turns the linear trend into curved, compounding growth on the original scale. ETS extends a straight-line trend. Both forecasts are plausible, but the ARIMA model has the better fit and white-noise residuals, so I prefer it.

Session info

sessionInfo()
## R version 4.6.1 (2026-06-24)
## Platform: aarch64-apple-darwin25.4.0
## Running under: macOS Tahoe 26.7.1
## 
## Matrix products: default
## BLAS:   /opt/homebrew/Cellar/openblas/0.3.34/lib/libopenblasp-r0.3.34.dylib 
## LAPACK: /opt/homebrew/Cellar/r/4.6.1/lib/R/lib/libRlapack.dylib;  LAPACK version 3.12.1
## 
## locale:
## [1] C.UTF-8/C.UTF-8/C.UTF-8/C/C.UTF-8/C.UTF-8
## 
## time zone: America/New_York
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] fable_0.5.0       feasts_0.5.0      fabletools_0.8.0  ggtime_1.0.0     
##  [5] tsibbledata_0.4.1 tsibble_1.2.0     ggplot2_4.0.3     lubridate_1.9.5  
##  [9] tidyr_1.3.2       dplyr_1.2.1       tibble_3.3.1      fpp3_1.0.3       
## 
## loaded via a namespace (and not attached):
##  [1] ggdist_3.3.3         sass_0.4.10          rappdirs_0.3.4      
##  [4] utf8_1.2.6           generics_0.1.4       anytime_0.3.13      
##  [7] lattice_0.22-9       digest_0.6.39        magrittr_2.0.5      
## [10] evaluate_1.0.5       grid_4.6.1           timechange_0.4.0    
## [13] RColorBrewer_1.1-3   mixtime_0.3.0        fastmap_1.2.0       
## [16] jsonlite_2.0.0       tinytex_0.60         purrr_1.2.2         
## [19] scales_1.4.0         jquerylib_0.1.4      numDeriv_2016.8-1.1 
## [22] cli_3.6.6            rlang_1.3.0          crayon_1.5.3        
## [25] vecvec_1.3.0         cachem_1.1.0         withr_3.0.3         
## [28] yaml_2.3.12          tools_4.6.1          tzdb_0.5.0          
## [31] vctrs_0.7.3          R6_2.6.1             lifecycle_1.0.5     
## [34] pkgconfig_2.0.3      urca_1.3-4           bslib_0.12.0        
## [37] progressr_1.0.0      pillar_1.11.1        gtable_0.3.6        
## [40] glue_1.8.1           Rcpp_1.1.2           xfun_0.60           
## [43] tidyselect_1.2.1     rstudioapi_0.19.0    knitr_1.52          
## [46] farver_2.1.2         htmltools_0.5.9      nlme_3.1-169        
## [49] rmarkdown_2.32       labeling_0.4.3       compiler_4.6.1      
## [52] S7_0.2.2             distributional_0.9.0