Homework 6

library(fpp3)
library(tidyverse)
library(patchwork)

---

Exercise 9.1 – ACF of White Noise Series

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

Part a – Explain the differences; do they all indicate white noise?

All three plots display white noise: no autocorrelation is significantly different from zero at any lag (all bars fall within the blue dashed critical-value bounds). However, the plots look different because the critical value bounds narrow as sample size grows.

The critical bounds are approximately ±1.96/√T. As T increases the bounds shrink:

Sample size (T)	Approx. bound
36	±0.327
360	±0.103
1000	±0.062

With only 36 observations the wide bands make it hard to detect genuine autocorrelation, whereas with 1,000 observations very small departures from zero would appear significant.

Part b – Why are critical values at different distances and autocorrelations different?

Different critical values: The bounds equal ±1.96/√T. As T grows, the standard error of each sample autocorrelation shrinks, so the bands narrow. This is a purely statistical effect – larger samples give more precise estimates.

Different autocorrelation values: Each plot is a different realisation of a white noise process. Sample autocorrelations are random variables. Even when the true autocorrelation is zero, a finite sample will produce non-zero estimates due to sampling variability. Larger samples produce estimates closer to the true value of zero, which is why the bars in the 1,000-observation plot are generally shorter than those in the 36-observation plot.

---

Exercise 9.2 – Amazon Stock Prices (Non-Stationarity)

Plot daily closing prices for Amazon, plus ACF and PACF. Explain why each shows non-stationarity.

amazon <- gafa_stock |>
  filter(Symbol == "AMZN") |>
  mutate(trading_day = row_number()) |>
  as_tsibble(index = trading_day)

p1 <- amazon |>
  autoplot(Close) +
  labs(title = "Amazon Daily Closing Price",
       x = "Trading Day", y = "Closing Price (USD)")

p2 <- amazon |>
  ACF(Close, lag_max = 50) |>
  autoplot() +
  labs(title = "ACF – Amazon Close Price")

p3 <- amazon |>
  PACF(Close, lag_max = 50) |>
  autoplot() +
  labs(title = "PACF – Amazon Close Price")

p1 / (p2 | p3)

Time plot: The series shows a clear long-term upward trend with no evidence of mean-reversion. A stationary series fluctuates around a constant mean, which this clearly does not.

ACF: The autocorrelations decay extremely slowly from a value close to 1 and remain large for many lags. This “long memory” signature is the classic hallmark of a non-stationary (unit-root) series. A stationary series would show ACF values that drop quickly to near zero.

PACF: There is a single very large spike at lag 1 (close to 1), with all subsequent partial autocorrelations near zero. This is consistent with a random walk / ARIMA(0,1,0) process, confirming that one round of differencing is needed to induce stationarity.

---

Exercise 9.3 – Box-Cox Transformation and Differencing for Stationarity

Part a – Turkish GDP (global_economy)

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

# Choose Box-Cox lambda
lambda_turkey <- turkey_gdp |>
  features(GDP, features = guerrero) |>
  pull(lambda_guerrero)
cat("Guerrero lambda for Turkish GDP:", round(lambda_turkey, 4), "\n")

## Guerrero lambda for Turkish GDP: 0.1572

p_raw <- turkey_gdp |>
  autoplot(GDP) + labs(title = "Turkish GDP (raw)")

p_bc <- turkey_gdp |>
  autoplot(box_cox(GDP, lambda_turkey)) +
  labs(title = paste0("Turkish GDP (Box-Cox λ=", round(lambda_turkey, 3), ")"))

p_diff <- turkey_gdp |>
  autoplot(difference(box_cox(GDP, lambda_turkey))) +
  labs(title = "First difference of Box-Cox transformed GDP")

(p_raw | p_bc) / p_diff

# Formal unit-root test on differenced series
turkey_gdp |>
  mutate(gdp_tr = box_cox(GDP, lambda_turkey)) |>
  features(difference(gdp_tr), unitroot_kpss)

## # A tibble: 1 × 3
##   Country kpss_stat kpss_pvalue
##   <fct>       <dbl>       <dbl>
## 1 Turkey     0.0889         0.1

Finding: A Box-Cox transformation with the Guerrero-selected λ stabilises variance. One round of regular differencing (d = 1) removes the trend and produces a stationary series, as confirmed by the KPSS test (p-value > 0.05 indicates stationarity).

Part b – Accommodation Takings in Tasmania (aus_accommodation)

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

lambda_tas <- tasmania |>
  features(Takings, features = guerrero) |>
  pull(lambda_guerrero)
cat("Guerrero lambda for Tasmania takings:", round(lambda_tas, 4), "\n")

## Guerrero lambda for Tasmania takings: 0.0018

tasmania |>
  autoplot(Takings) +
  labs(title = "Tasmania Accommodation Takings (raw)")

tasmania |>
  autoplot(box_cox(Takings, lambda_tas)) +
  labs(title = paste0("Tasmania Takings (Box-Cox λ=", round(lambda_tas, 3), ")"))

# Seasonal differencing first, then regular
tasmania |>
  mutate(tak_tr = box_cox(Takings, lambda_tas)) |>
  features(tak_tr, unitroot_nsdiffs)

## # A tibble: 1 × 2
##   State    nsdiffs
##   <chr>      <int>
## 1 Tasmania       1

tasmania |>
  mutate(tak_tr  = box_cox(Takings, lambda_tas),
         tak_sd  = difference(tak_tr, lag = 4)) |>
  features(tak_sd, unitroot_ndiffs)

## # A tibble: 1 × 2
##   State    ndiffs
##   <chr>     <int>
## 1 Tasmania      0

tasmania |>
  mutate(tak_tr = box_cox(Takings, lambda_tas),
         tak_sd = difference(tak_tr, lag = 4),
         tak_d  = difference(tak_sd)) |>
  autoplot() +
  labs(title = "Tasmania Takings – seasonally and first-differenced (Box-Cox)")

Finding: The series has strong seasonality (quarterly data). The Guerrero λ stabilises variance. unitroot_nsdiffs() recommends one seasonal difference (D = 1, lag = 4), and unitroot_ndiffs() on the seasonally differenced series recommends one additional regular difference (d = 1), giving ARIMA(p, 1, q)(P, 1, Q)[4].

Part c – Monthly Souvenir Sales (souvenirs)

lambda_souv <- souvenirs |>
  features(Sales, features = guerrero) |>
  pull(lambda_guerrero)
cat("Guerrero lambda for souvenirs:", round(lambda_souv, 4), "\n")

## Guerrero lambda for souvenirs: 0.0021

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

souvenirs |>
  autoplot(box_cox(Sales, lambda_souv)) +
  labs(title = paste0("Souvenir Sales (Box-Cox λ=", round(lambda_souv, 3), ")"))

# Check seasonal differencing
souvenirs |>
  mutate(s_tr = box_cox(Sales, lambda_souv)) |>
  features(s_tr, unitroot_nsdiffs)

## # A tibble: 1 × 1
##   nsdiffs
##     <int>
## 1       1

souvenirs |>
  mutate(s_tr = box_cox(Sales, lambda_souv),
         s_sd = difference(s_tr, lag = 12)) |>
  features(s_sd, unitroot_ndiffs)

## # A tibble: 1 × 1
##   ndiffs
##    <int>
## 1      0

souvenirs |>
  mutate(s_tr = box_cox(Sales, lambda_souv),
         s_sd = difference(s_tr, lag = 12)) |>
  autoplot(s_sd) +
  labs(title = "Souvenir Sales – seasonally differenced (Box-Cox)")

Finding: The Guerrero-selected λ is close to 0, suggesting a log transformation is appropriate (it stabilises the growing seasonal variance). One seasonal difference (D = 1, lag = 12) is sufficient; no additional regular differencing is needed.

---

Exercise 9.5 – Retail Data: Order of Differencing

Using the retail series from Exercise 7 of Section 2.10.

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

myseries |> autoplot(Turnover) +
  labs(title = "Retail Turnover – selected series")

# Box-Cox
lambda_retail <- myseries |>
  features(Turnover, features = guerrero) |>
  pull(lambda_guerrero)
cat("Guerrero lambda for retail:", round(lambda_retail, 4), "\n")

## Guerrero lambda for retail: 0.083

myseries |>
  autoplot(box_cox(Turnover, lambda_retail)) +
  labs(title = paste0("Retail Turnover – Box-Cox (λ=", round(lambda_retail, 3), ")"))

# Number of seasonal differences
myseries |>
  mutate(tr = box_cox(Turnover, lambda_retail)) |>
  features(tr, unitroot_nsdiffs)

## # A tibble: 1 × 3
##   State              Industry                                            nsdiffs
##   <chr>              <chr>                                                 <int>
## 1 Northern Territory Clothing, footwear and personal accessory retailing       1

# Number of regular differences after seasonal
myseries |>
  mutate(tr  = box_cox(Turnover, lambda_retail),
         str = difference(tr, lag = 12)) |>
  features(str, unitroot_ndiffs)

## # A tibble: 1 × 3
##   State              Industry                                            ndiffs
##   <chr>              <chr>                                                <int>
## 1 Northern Territory Clothing, footwear and personal accessory retailing      0

# Final stationary series
myseries |>
  mutate(tr   = box_cox(Turnover, lambda_retail),
         str  = difference(tr, lag = 12),
         strd = difference(str)) |>
  autoplot(strd) +
  labs(title = "Retail Turnover – Box-Cox + seasonal diff + regular diff")

Finding: After applying the Box-Cox transformation, one seasonal difference and one regular difference yield a stationary series. This suggests a seasonal ARIMA model with D = 1 and d = 1 on the transformed scale.

---

Exercise 9.6 – Simulating ARIMA Models

Part a & b – AR(1) with φ₁ = 0.6

set.seed(42)

# AR(1) with phi1 = 0.6
sim_ar1 <- function(phi1, n = 100) {
  y <- numeric(n)
  e <- rnorm(n)
  for (i in 2:n) y[i] <- phi1 * y[i-1] + e[i]
  tsibble(idx = seq_len(n), y = y, index = idx)
}

sim06  <- sim_ar1(0.6)
sim09  <- sim_ar1(0.9)
simneg <- sim_ar1(-0.6)

p1 <- sim06  |> autoplot(y) + labs(title = "AR(1): φ₁ = 0.6")
p2 <- sim09  |> autoplot(y) + labs(title = "AR(1): φ₁ = 0.9")
p3 <- simneg |> autoplot(y) + labs(title = "AR(1): φ₁ = -0.6")

p1 / p2 / p3

Effect of φ₁: As φ₁ increases toward 1 the series becomes smoother and more persistent (mean-reversion is slower). A large positive φ₁ (e.g., 0.9) produces long, sweeping runs above or below the mean. A negative φ₁ (e.g., -0.6) causes the series to oscillate rapidly around the mean.

Part c & d – MA(1) with θ₁ = 0.6

sim_ma1 <- function(theta1, n = 100) {
  y <- numeric(n)
  e <- rnorm(n + 1)
  for (i in 1:n) y[i] <- e[i+1] + theta1 * e[i]
  tsibble(idx = seq_len(n), y = y, index = idx)
}

ma06  <- sim_ma1(0.6)
ma09  <- sim_ma1(0.9)
maneg <- sim_ma1(-0.6)

q1 <- ma06  |> autoplot(y) + labs(title = "MA(1): θ₁ = 0.6")
q2 <- ma09  |> autoplot(y) + labs(title = "MA(1): θ₁ = 0.9")
q3 <- maneg |> autoplot(y) + labs(title = "MA(1): θ₁ = -0.6")

q1 / q2 / q3

Effect of θ₁: All MA(1) series are stationary regardless of θ₁. A large positive θ₁ produces positively correlated successive values (smoother), while a negative θ₁ creates negatively correlated successive values (more jagged / oscillatory). The overall amplitude looks similar across values of θ₁.

Part e – ARMA(1,1)

set.seed(42)
n <- 100
y_arma <- numeric(n)
e <- rnorm(n + 1)
for (i in 2:n)
  y_arma[i] <- 0.6 * y_arma[i-1] + e[i] + 0.6 * e[i-1]

arma_sim <- tsibble(idx = seq_len(n), y = y_arma, index = idx)
arma_sim |> autoplot(y) +
  labs(title = "ARMA(1,1): φ₁ = 0.6, θ₁ = 0.6")

Part f – AR(2) with non-stationary parameters

set.seed(42)
y_ar2 <- numeric(100)
e2 <- rnorm(100)
for (i in 3:100)
  y_ar2[i] <- -0.8 * y_ar2[i-1] + 0.3 * y_ar2[i-2] + e2[i]

ar2_sim <- tsibble(idx = seq_len(100), y = y_ar2, index = idx)
ar2_sim |> autoplot(y) +
  labs(title = "AR(2): φ₁ = -0.8, φ₂ = 0.3  (non-stationary region)")

Part g – Comparing ARMA(1,1) and AR(2)

p_arma <- arma_sim |> autoplot(y) +
  labs(title = "ARMA(1,1): φ₁=0.6, θ₁=0.6")
p_ar2  <- ar2_sim  |> autoplot(y) +
  labs(title = "AR(2): φ₁=-0.8, φ₂=0.3")

p_arma | p_ar2

Comparison:

• The ARMA(1,1) series is stationary; it oscillates around zero with moderate autocorrelation and no explosive behaviour.

• The AR(2) series with φ₁ = −0.8, φ₂ = 0.3 lies outside the stationarity triangle (the roots of 1 + 0.8B − 0.3B² = 0 lie inside the unit circle). This produces a rapidly oscillating, non-stationary series that can grow without bound over time.

---

Exercise 9.7 – Australian Air Passengers

air <- aus_airpassengers |>
  as_tsibble()

air |> autoplot(Passengers) +
  labs(title = "Total Australian Air Passengers (millions)")

Part a – Automatic ARIMA

fit_auto <- air |>
  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

# Residuals time plot
augment(fit_auto) |>
  autoplot(.innov) +
  labs(title = "Residuals – Auto ARIMA")

# ACF of residuals
augment(fit_auto) |>
  ACF(.innov) |>
  autoplot() +
  labs(title = "ACF of Residuals – Auto ARIMA")

# Ljung-Box test
augment(fit_auto) |>
  features(.innov, ljung_box, lag = 10, dof = 2)

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

# Forecasts
fit_auto |>
  forecast(h = 10) |>
  autoplot(air) +
  labs(title = "Auto ARIMA Forecasts – Next 10 Periods")

Selected model: ARIMA() selects an appropriate model based on AICc minimisation. Check the report() output for the exact specification. The residuals should look like white noise (no significant ACF spikes, roughly normal distribution).

Part b – Backshift operator notation

For an ARIMA(p, d, q) model, the general form using the backshift operator B is:

φ(B)(1 − B)^d y_t = c + θ(B) ε_t

where φ(B) = 1 − φ₁B − … − φ_pB^p and θ(B) = 1 + θ₁B + … + θ_qB^q.

(Substitute the specific p, d, q and coefficient values from report() above into this template.)

Part c – ARIMA(0,1,0) with drift

fit_rw <- air |>
  model(ARIMA(Passengers ~ 1 + pdq(0,1,0)))

fit_rw |>
  forecast(h = 10) |>
  autoplot(air) +
  labs(title = "ARIMA(0,1,0) with Drift – Forecasts")

The ARIMA(0,1,0) with drift is simply a random walk with drift. Forecasts are a straight line with constant upward slope equal to the average period- to-period increment. Compared to the auto-selected model, the prediction intervals may be wider or narrower depending on how much the AR/MA terms reduce residual variance.

Part d – ARIMA(2,1,2) with and without constant

fit_212 <- air |>
  model(
    with_const    = ARIMA(Passengers ~ 1 + pdq(2,1,2)),
    without_const = ARIMA(Passengers ~ 0 + pdq(2,1,2))
  )

fit_212 |>
  forecast(h = 10) |>
  autoplot(air) +
  facet_wrap(~.model, ncol = 1) +
  labs(title = "ARIMA(2,1,2): with vs without constant")

With constant: forecasts include a deterministic drift term. Without constant: forecasts revert toward the long-run mean of differences (zero), so the level forecasts eventually flatten out relative to the trend-bearing version.

Part e – ARIMA(0,2,1) with constant

fit_021 <- air |>
  model(ARIMA(Passengers ~ 1 + pdq(0,2,1)))

fit_021 |>
  forecast(h = 10) |>
  autoplot(air) +
  labs(title = "ARIMA(0,2,1) with Constant")

What happens: Including a constant in a twice-differenced model (d = 2) introduces a quadratic trend into the forecasts, which can cause them to curve upward quite steeply. This is often undesirable and is why ARIMA() typically suppresses the constant when d ≥ 2.

---

Exercise 9.8 – US GDP

Part a – Box-Cox transformation

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

us_gdp |> autoplot(GDP) +
  labs(title = "United States GDP (raw)")

lambda_us <- us_gdp |>
  features(GDP, features = guerrero) |>
  pull(lambda_guerrero)
cat("Guerrero lambda:", round(lambda_us, 4), "\n")

## Guerrero lambda: 0.2819

us_gdp |>
  autoplot(box_cox(GDP, lambda_us)) +
  labs(title = paste0("US GDP – Box-Cox (λ=", round(lambda_us, 3), ")"))

Part b – Fit ARIMA to transformed data

fit_us <- us_gdp |>
  model(arima = ARIMA(box_cox(GDP, lambda_us)))

report(fit_us)

## Series: GDP 
## Model: ARIMA(1,1,0) w/ drift 
## Transformation: box_cox(GDP, lambda_us) 
## 
## Coefficients:
##          ar1  constant
##       0.4586  118.1822
## s.e.  0.1198    9.5047
## 
## sigma^2 estimated as 5479:  log likelihood=-325.32
## AIC=656.65   AICc=657.1   BIC=662.78

Part c – Try alternative models

fit_us_cands <- us_gdp |>
  model(
    auto       = ARIMA(box_cox(GDP, lambda_us)),
    arima011   = ARIMA(box_cox(GDP, lambda_us) ~ pdq(0,1,1)),
    arima110   = ARIMA(box_cox(GDP, lambda_us) ~ pdq(1,1,0)),
    arima111   = ARIMA(box_cox(GDP, lambda_us) ~ pdq(1,1,1)),
    arima012   = ARIMA(box_cox(GDP, lambda_us) ~ pdq(0,1,2))
  )

glance(fit_us_cands) |>
  select(.model, AIC, AICc, BIC) |>
  arrange(AICc)

## # A tibble: 5 × 4
##   .model     AIC  AICc   BIC
##   <chr>    <dbl> <dbl> <dbl>
## 1 auto      657.  657.  663.
## 2 arima110  657.  657.  663.
## 3 arima011  659.  659.  665.
## 4 arima111  659.  659.  667.
## 5 arima012  659.  660.  667.

Part d – Best model and residual diagnostics

# Residual plot
augment(fit_us) |>
  autoplot(.innov) +
  labs(title = "Residuals – Best US GDP ARIMA")

# ACF of residuals
augment(fit_us) |>
  ACF(.innov) |>
  autoplot() +
  labs(title = "ACF of Residuals")

# Ljung-Box test
augment(fit_us) |>
  features(.innov, ljung_box, lag = 10, dof = 2)

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

Interpretation: If the Ljung-Box p-value is > 0.05, the residuals are consistent with white noise. Check the ACF of residuals – no significant spikes should remain.

Part e – Forecasts

fit_us |>
  forecast(h = 10) |>
  autoplot(us_gdp) +
  labs(title = "US GDP – ARIMA Forecasts (10 years ahead, Box-Cox scale)")

The forecasts follow the established trend. The prediction intervals widen over the forecast horizon, reflecting increasing uncertainty.

Part f – Compare with ETS (no transformation)

fit_comparison <- us_gdp |>
  model(
    ARIMA = ARIMA(box_cox(GDP, lambda_us)),
    ETS   = ETS(GDP)
  )

fit_comparison |>
  forecast(h = 10) |>
  autoplot(us_gdp) +
  facet_wrap(~.model, ncol = 1) +
  labs(title = "US GDP: ARIMA (Box-Cox) vs ETS (no transformation)")

glance(fit_comparison) |> select(.model, AIC, AICc, BIC)

## # A tibble: 2 × 4
##   .model   AIC  AICc   BIC
##   <chr>  <dbl> <dbl> <dbl>
## 1 ARIMA   657.  657.  663.
## 2 ETS    3191. 3192. 3201.

Comparison: ETS typically selects an exponential smoothing model that captures the trend. Both approaches should produce broadly similar forecasts for a series with a smooth upward trend. Differences arise in: (1) how transformation affects variance stabilisation, (2) whether the uncertainty fan widens differently. Comparing AICc values (on the same scale) or cross-validation accuracy would inform the better model choice.

---

End of Chapter 9 ARIMA Exercises