Classswork: Hyndman Exercises, Chapter 9, Questions 9.1 - 9.3 and 9.5 - 9.8

library(tsibble)
## Registered S3 method overwritten by 'tsibble':
##   method               from 
##   as_tibble.grouped_df dplyr
## 
## Attaching package: 'tsibble'
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, union
library(fable)
## Loading required package: fabletools
library(fpp3)
## ── Attaching packages ──────────────────────────────────────────── fpp3 1.0.2 ──
## ✔ tibble      3.3.0     ✔ ggplot2     4.0.0
## ✔ dplyr       1.1.4     ✔ tsibbledata 0.4.1
## ✔ tidyr       1.3.1     ✔ feasts      0.4.2
## ✔ lubridate   1.9.4
## ── Conflicts ───────────────────────────────────────────────── fpp3_conflicts ──
## ✖ lubridate::date()     masks base::date()
## ✖ dplyr::filter()       masks stats::filter()
## ✖ tsibble::intersect()  masks base::intersect()
## ✖ lubridate::interval() masks tsibble::interval()
## ✖ dplyr::lag()          masks stats::lag()
## ✖ tsibble::setdiff()    masks base::setdiff()
## ✖ tsibble::union()      masks base::union()
library(ggplot2)

Question 9.1a

Explain the differences among three types of ACFs for 36 random number, 360 randon numbers and 1000 random numbers.

Question 1a Response:

The primary difference among the three ACF plots shown in Figure 9.11 relates to the width of the confidence intervals represented by the dotted blue lines, with the widest interval appearing for 36 random numbers, a narrower interval for 360 random numbers, and the narrowest for 1000 random numbers. This progressive narrowing follows the mathematical formula (plus or minus 1.96 divided by the square root of N), where N represents the number of observations, creating an inverse square root relationship between sample size and interval width.

With only 36 observations, relatively large autocorrelations are needed to be considered statistically significant, while with 1000 observations, even small autocorrelations can be detected as meaningful.

Across all three ACF plots, the autocorrelation bars consistently remain inside their respective confidence intervals, indicating no significant autocorrelation at any lag and confirming that all three datasets represent white noise where past values do not provide useful information for predicting future values.

Question 9.1b

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?

Question 9.1b Response

The critical values appear at different distances from zero across the three ACF plots because each is based on a different number of observations, and the critical value calculation (plus or minus 1.96 divided by the square root of N) directly incorporates sample size, causing the boundaries to move progressively closer to zero as N increases from 36 to 360 to 1000. This convergence occurs because larger sample sizes fundamentally increase statistical precision, allowing us to more confidently distinguish true patterns from random noise and eliminating the need for wide confidence bands to account for uncertainty.

Each of the three samples represents a different random realization of white noise, displaying different patterns of autocorrelation values purely due to random sampling variability.

Question 9.2

Plot the daily closing prices for Amazon stock (contained in gafa_stock), along with the ACF and PACF.

# Load required libraries
library(fpp3)  # or library(tsibble) and library(feasts)
library(ggplot2)

# Filter for Amazon stock
amazon <- gafa_stock %>%
  filter(Symbol == "AMZN")

# Plot 1: Daily closing prices
amazon %>%
  autoplot(Close) +
  labs(title = "Amazon Daily Closing Prices",
       y = "Closing Price ($)",
       x = "Date")

# Plot 2: ACF
amazon %>%
  ACF(Close) %>%
  autoplot() +
  labs(title = "ACF of Amazon Closing Prices")

# Plot 3: PACF
amazon %>%
  PACF(Close) %>%
  autoplot() +
  labs(title = "PACF of Amazon Closing Prices")

# Alternative: All three plots together
library(patchwork)

p1 <- amazon %>% autoplot(Close) + labs(title = "Amazon Closing Prices")
p2 <- amazon %>% ACF(Close) %>% autoplot() + labs(title = "ACF")
p3 <- amazon %>% PACF(Close) %>% autoplot() + labs(title = "PACF")

p1 / (p2 | p3)

Question 9.2 Response

A non-stationary series shows trends, changing variance, or seasonality that evolve over time. The time series plot of Amazon’s daily closing prices clearly shows non-stationarity through a pronounced upward trend over the observed period, indicating that the mean level of stock prices has systematically increased rather than fluctuating around a constant mean. In contrast, a stationary series would fluctuate randomly around a stable mean with no systematic directional movement.

The ACF and PACF plots provide statistical evidence of non-stationarity. The ACF displays a slow decay pattern where autocorrelations begin at 1.0 at lag 0 and remain highly significant, extending well beyond the upper confidence bound for many lags before declining very gradually. This persistent, gradual decline is a classic signature of non-stationary data, whereas a stationary series would show rapid decay to zero within the first few lags.

The PACF plot further confirms non-stationarity by displaying a large, dominant spike at lag 1 that approaches 1.0, followed by a sharp drop where subsequent lags fall mostly within the confidence bounds. This pattern—a near-unity spike at lag 1 with rapid cutoff—is characteristic of a random walk or unit root process, suggesting that each day’s closing price is almost entirely determined by the previous day’s price plus a random shock, causing the series to drift without reverting to a stable mean.

Question 9.3a Find an appropriate Box-Cox transformation and order of differencing in order to obtain stationary data - Turkish GDP

# Turkish GDP
turkey_gdp <- global_economy %>%
  filter(Country == "Turkey")

# Find lambda
lambda <- turkey_gdp %>%
  features(GDP, features = guerrero) %>%
  pull(lambda_guerrero)

# PRINT THE LAMBDA VALUE
print(paste("The optimal lambda:", lambda))
## [1] "The optimal lambda: 0.157218680108172"
# Apply transformation and differencing
turkey_gdp %>%
  mutate(diff_gdp = difference(box_cox(GDP, lambda))) %>%
  autoplot(diff_gdp)
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_line()`).

# Check ACF
turkey_gdp %>%
  mutate(diff_gdp = difference(box_cox(GDP, lambda))) %>%
  ACF(diff_gdp) %>%
  autoplot()

Question 9.3a Response

For Turkish GDP, the optimal Box-Cox transformation parameter using the Guerrero method is lambda = 0.157. Applying this transformation with first-order differencing (d=1) produces a stationary series, as evidenced by the ACF plot. The ACF plot shows the following characteristics of stationarity: - The values oscillate randomly around zero

Question 9.3a plots to show the effect of differencing and Box-Cox transformation

library(patchwork)

# Original series
p1 <- turkey_gdp %>%
  autoplot(GDP) +
  labs(title = "Original Turkish GDP (Non-stationary)")

# After transformation and differencing
p2 <- turkey_gdp %>%
  mutate(diff_gdp = difference(box_cox(GDP, lambda))) %>%
  autoplot(diff_gdp) +
  labs(title = "Transformed & Differenced (Stationary)")

# Combine plots
p1 / p2
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_line()`).

Question 9.3b Find an appropriate Box-Cox transformation and order of differencing in order to obtain stationary data - Accomodation takings in the state of Tasmania

# Tasmania Accommodation
tasmania <- aus_accommodation %>%
  filter(State == "Tasmania")

# Find lambda
lambda_tas <- tasmania %>%
  features(Takings, features = guerrero) %>%
  pull(lambda_guerrero)

# PRINT THE LAMBDA VALUE
print(paste("The optimal lambda:", lambda_tas))
## [1] "The optimal lambda: 0.00181964338451565"
# Apply transformation and differencing
tasmania %>%
  mutate(diff_takings = difference(box_cox(Takings, lambda_tas))) %>%
  autoplot(diff_takings)
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_line()`).

# Check ACF
tasmania %>%
  mutate(diff_takings = difference(box_cox(Takings, lambda_tas))) %>%
  ACF(diff_takings) %>%
  autoplot()

Question 9.3b Response

For Tasmania accommodation takings, the optimal Box-Cox transformation parameter using the Guerrero method is lambda = 0.00181. Applying this transformation with first-order differencing (d=1) produces a stationary series, as evidenced by the ACF plot. The ACF plot shows the following characteristics of stationarity:

Question 9.3b plots to show the effect of differencing and Box-Cox transformation

library(patchwork)

# Original series
p1 <- tasmania %>%
  autoplot(Takings) +
  labs(title = "Original Tasmania Accommodation Takings (Non-stationary)")

# After transformation and differencing
p2 <- tasmania %>%
  mutate(diff_takings = difference(box_cox(Takings, lambda_tas))) %>%
  autoplot(diff_takings) +
  labs(title = "Transformed & Differenced (Stationary)")

# Combine plots
p1 / p2
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_line()`).

#### Question 9.3b

Prior to the transformation you can clearly see an upward trend and seasonal pattern in the original series plot, indicating non-stationarity. The variance also changes from 2000 to 2015. After transformation, the variance is more stableand the graph does not exhibit any obvious trend or seasonal pattern, suggesting that the series is now stationary and the transformation was successful.

Question 9.3c Find an appropriate Box-Cox transformation and order of differencing in order to obtain stationary data - sales data from the souvenirs dataset

# Souvenirs Sales
# Find lambda
sales_from_souvenirs <- souvenirs %>%
  features(Sales, features = guerrero) %>%
  pull(lambda_guerrero)

# PRINT THE LAMBDA VALUE
print(paste("The optimal lambda:", sales_from_souvenirs))
## [1] "The optimal lambda: 0.00211822065057585"
# Apply transformation and differencing
souvenirs %>%
  mutate(diff_sales = difference(box_cox(Sales, sales_from_souvenirs))) %>%
  autoplot(diff_sales)
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_line()`).

# Check ACF
souvenirs %>%
  mutate(diff_sales = difference(box_cox(Sales, sales_from_souvenirs))) %>%
  ACF(diff_sales) %>%
  autoplot()

####Question 9.3c plots to show the effect of differencing and Box-Cox transformation

library(patchwork)

# Original series
p1 <- souvenirs %>%
  autoplot(Sales) +
  labs(title = "Original Monthly Souvenir Sales (Non-stationary)")

# After transformation and differencing
p2 <- souvenirs %>%
  mutate(diff_sales = difference(box_cox(Sales, sales_from_souvenirs))) %>%
  autoplot(diff_sales) +
  labs(title = "Transformed & Differenced (Stationary)")

# Combine plots
p1 / p2
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_line()`).

#### Question 9.3c Response

# Souvenirs Sales

# Find lambda
lambda_souv <- souvenirs %>%
  features(Sales, features = guerrero) %>%
  pull(lambda_guerrero)

# PRINT THE LAMBDA VALUE
print(paste("The optimal lambda:", lambda_souv))
## [1] "The optimal lambda: 0.00211822065057585"
# Apply transformation with SEASONAL differencing (lag=12 for monthly data)
souvenirs %>%
  mutate(diff_sales = difference(box_cox(Sales, lambda_souv), lag = 12)) %>%
  autoplot(diff_sales) +
  labs(title = "Seasonally Differenced (D=1, lag=12)")
## Warning: Removed 12 rows containing missing values or values outside the scale range
## (`geom_line()`).

# Check ACF
souvenirs %>%
  mutate(diff_sales = difference(box_cox(Sales, lambda_souv), lag = 12)) %>%
  ACF(diff_sales) %>%
  autoplot()

# If still not stationary, apply BOTH seasonal and regular differencing
souvenirs %>%
  mutate(
    seasonal_diff = difference(box_cox(Sales, lambda_souv), lag = 12),
    final_diff = difference(seasonal_diff, lag = 1)
  ) %>%
  autoplot(final_diff) +
  labs(title = "Seasonal (D=1, lag=12) + Regular (d=1) Differencing")
## Warning: Removed 13 rows containing missing values or values outside the scale range
## (`geom_line()`).

# Check ACF for final differencing
souvenirs %>%
  mutate(
    seasonal_diff = difference(box_cox(Sales, lambda_souv), lag = 12),
    final_diff = difference(seasonal_diff, lag = 1)
  ) %>%
  ACF(final_diff) %>%
  autoplot()

library(patchwork)

# Original series
p1 <- souvenirs %>%
  autoplot(Sales) +
  labs(title = "Original Monthly Souvenir Sales (Non-stationary)")

# After Box-Cox and regular differencing only (shows seasonality)
p2 <- souvenirs %>%
  mutate(diff_sales = difference(box_cox(Sales, lambda_souv), lag = 1)) %>%
  autoplot(diff_sales) +
  labs(title = "Transformed & Differenced d=1 (Still Seasonal)")

# After Box-Cox, seasonal, and regular differencing (stationary)
p3 <- souvenirs %>%
  mutate(
    seasonal_diff = difference(box_cox(Sales, lambda_souv), lag = 12),
    final_diff = difference(seasonal_diff, lag = 1)
  ) %>%
  autoplot(final_diff) +
  labs(title = "Transformed, D=1 (lag=12), d=1 (Stationary)")

# Combine plots
p1 / p2 / p3
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 13 rows containing missing values or values outside the scale range
## (`geom_line()`).

####Question 9.3c Response

For monthly souvenir sales, the optimal Box-Cox transformation parameter using the Guerrero method is lambda = 0.00181. The initial first-order differencing (d=1) revealed a strong seasonal pattern with regular sharp drops every January. Therefore, seasonal differencing with lag=12 was applied to remove the yearly seasonality, followed by first-order differencing (d=1) to remove any remaining trend. This combination of Box-Cox transformation (lambda=0.00181), seasonal differencing (D=1, lag=12), and regular differencing (d=1) produces a stationary series, as evidenced by the ACF plot showing:

Question 9.5 For the retail data aus_retail, find the appropriate order of differencing (after transformation if necessary) to obtain stationary data.

# For example, I will choose Victoria Takeaway food services
retail_series <- aus_retail %>%
  filter(State == "Victoria", 
         Industry == "Takeaway food services")

# Visualize original data (very importnat to. get a sense of the data)
retail_series %>%
  autoplot(Turnover) +
  labs(title = "Original Retail Turnover")

# Check if transformation needed
lambda <- retail_series %>%
  features(Turnover, features = guerrero) %>%
  pull(lambda_guerrero)

print(paste("Optimal lambda:", lambda))
## [1] "Optimal lambda: 0.0887219603737718"

Question 9.5 Response

# Apply transformation and check
retail_series %>%
  autoplot(box_cox(Turnover, lambda)) +
  labs(title = "Box-Cox Transformed Turnover")

# Apply seasonal differencing first (lag=12 for monthly)
retail_series %>%
  mutate(diff_turnover = difference(box_cox(Turnover, lambda), lag = 12)) %>%
  autoplot(diff_turnover) +
  labs(title = "Seasonally Differenced (D=1)")
## Warning: Removed 12 rows containing missing values or values outside the scale range
## (`geom_line()`).

# Check ACF after seasonal differencing
retail_series %>%
  mutate(diff_turnover = difference(box_cox(Turnover, lambda), lag = 12)) %>%
  ACF(diff_turnover) %>%
  autoplot()

# Step 8: If needed, add regular differencing
retail_series %>%
  mutate(
    seasonal_diff = difference(box_cox(Turnover, lambda), lag = 12),
    final_diff = difference(seasonal_diff, lag = 1)
  ) %>%
  autoplot(final_diff) +
  labs(title = "Seasonal (D=1) + Regular (d=1) Differencing")
## Warning: Removed 13 rows containing missing values or values outside the scale range
## (`geom_line()`).

# Step 9: Final ACF check
retail_series %>%
  mutate(
    seasonal_diff = difference(box_cox(Turnover, lambda), lag = 12),
    final_diff = difference(seasonal_diff, lag = 1)
  ) %>%
  ACF(final_diff) %>%
  autoplot()

#### Question 9.5 Response

After the Box-Cox transformation (lambda = 0.089), the time series still exhibits a clear upward trend, indicating that it remains non-stationary. After applying both seasonal differencing (D=1, lag=12) and regular differencing (d=1), the resulting series appears stationary. The ACF plot shows an initial steep drop in values, and after lag 12, the values oscillate randomly around zero with most autocorrelation bars falling within the confidence bands. Thus, the appropriate order of differencing for this retail data is D=1 (seasonal, lag=12) and d=1 (regular).

Question 9.6a

Use code to generate data from an AR(1) model

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)

Question 9.6b

Produce a time plot for the series

# Question 9.6.2: Time plot of AR(1)

sim %>%
  autoplot(y) +
  labs(title = "AR(1) Model with phi = 0.6",
       y = "y",
       x = "Time")

# Experiment with different phi values
# Try phi = 1.5
y2 <- numeric(100)
e2 <- rnorm(100)
for(i in 2:100)
  y2[i] <- 1.5*y2[i-1] + e2[i]

sim2 <- tsibble(idx = seq_len(100), y = y2, index = idx)

sim2 %>%
  autoplot(y) +
  labs(title = "AR(1) Model with phi = 1.5")

# Try phi = 0.9
y3 <- numeric(100)
e3 <- rnorm(100)
for(i in 2:100)
  y3[i] <- 0.9*y3[i-1] + e3[i]

sim3 <- tsibble(idx = seq_len(100), y = y3, index = idx)

sim3 %>%
  autoplot(y) +
  labs(title = "AR(1) Model with phi = 0.9")

#### Question 9.6b Response

When phi was 0.6 the graph showed a stationary series with fluctuations around a constant mean. When phi was increased to 0.9 the series still appeared to be trending towards non stationary as the variance beagn to increase as you moved to the right of the graph. There were sharper fluctuations and longer periods above and below the mean. When phi was set to 1.5 the series became non-stationary, exhibiting a clear upward, exponential trend with increasing variance over time.

Question 9.6c Write your own code to generate data from an MA(1) model with phi = 0.6 and sigma^2 = 1

y_ma <- numeric(100)   # Create empty vector
e_ma <- rnorm(100)     # Generate random errors

# First value
y_ma[1] <- e_ma[1]

# MA(1) formula
for(i in 2:100)
  y_ma[i] <- e_ma[i] + 0.6*e_ma[i-1]  # Current error + 60% of previous error

sim_ma <- tsibble(idx = seq_len(100), y = y_ma, index = idx)

Question 9.6d Produce a time plot for the series. How does the plot change as you change phi-one (words are used because unicode characters are nt setup for use with LaTex and give errors when the code is run.)

# Make the plot
sim_ma %>%
  autoplot(y) +
  labs(title = "MA(1) Model with theta = 0.6",
       subtitle = "Moving Average coefficient = 0.6")

# Experiment with different theta values
# Try theta = 0.2
y_ma2 <- numeric(100)
e_ma2 <- rnorm(100)
y_ma2[1] <- e_ma2[1]
for(i in 2:100)
  y_ma2[i] <- e_ma2[i] + 0.2*e_ma2[i-1]

sim_ma2 <- tsibble(idx = seq_len(100), y = y_ma2, index = idx)

sim_ma2 %>%
  autoplot(y) +
  labs(title = "MA(1) Model with theta = 0.2",
       subtitle = "Moving Average coefficient = 0.2")

# Try theta = 0.9
y_ma3 <- numeric(100)
e_ma3 <- rnorm(100)
y_ma3[1] <- e_ma3[1]
for(i in 2:100)
  y_ma3[i] <- e_ma3[i] + 0.9*e_ma3[i-1]

sim_ma3 <- tsibble(idx = seq_len(100), y = y_ma3, index = idx)

sim_ma3 %>%
  autoplot(y) +
  labs(title = "MA(1) Model with theta = 0.9",
       subtitle = "Moving Average coefficient = 0.9")

#### Question 9.6d Response

As the moving average coefficient increases from 0.2 to 0.9, the MA(1) series shows subtle differences in behavior. All three plots display rapid, erratic fluctuations around zero, characteristic of MA models. However, higher values of theta (such as 0.9) produce slightly smoother transitions between consecutive points, as each value incorporates more influence from the previous error term. Despite this, MA(1) models remain highly volatile compared to AR models, with limited persistence beyond a single lag. The series fluctuates randomly around a mean of zero with no visible trend, confirming the stationary nature of MA(1) processes

Question 9.6e Generate data from an ARMA(1,1) model with phi one equals 0.6, theta one equals 0.6, and sigma squared equals 1.[Combining AR(1) with MA(1)]

# Question 9.6e ARMA(1,1) with phi = 0.6, theta = 0.6

y_arma <- numeric(100)
e_arma <- rnorm(100)

y_arma[1] <- e_arma[1]

for(i in 2:100)
  y_arma[i] <- 0.6*y_arma[i-1] + e_arma[i] + 0.6*e_arma[i-1]

sim_arma <- tsibble(idx = seq_len(100), y = y_arma, index = idx)

sim_arma %>%
  autoplot(y) +
  labs(title = "ARMA(1,1) Model with phi=0.6, theta=0.6")

# Question 9.6e AR(2) with phi_1= -0.8 and phi_2 = 0.3

y_ar2 <- numeric(100)
e_ar2 <- rnorm(100)

y_ar2[1] <- e_ar2[1]
y_ar2[2] <- e_ar2[2]

for(i in 3:100)
  y_ar2[i] <- -0.8*y_ar2[i-1] + 0.3*y_ar2[i-2] + e_ar2[i]

sim_ar2 <- tsibble(idx = seq_len(100), y = y_ar2, index = idx)

sim_ar2 %>%
  autoplot(y) +
  labs(title = "AR(2) Model with phi_1 = -0.8, phi_2 = 0.3 (Non-stationary)")

#### Question 9.6 Compare the graph of ARMA(1,1) and AR(2)

# Question 9.6.7: Compare ARMA(1,1) and AR(2)

library(patchwork)

p1 <- sim_arma %>%
  autoplot(y) +
  labs(title = "ARMA(1,1): phi_1 = 0.6, theta_1 =0.6")

p2 <- sim_ar2 %>%
  autoplot(y) +
  labs(title = "AR(2): phi_1 = -0.8, phi_2 =0.3")

p1 / p2

Question 9.6 Explanantion of the ARMA(1,1) and AR(2) graphs

The two graphs demonstrate markedly different behavior. The ARMA(1,1) model with phi = 0.6 and theta = 0.6 exhibits stationary characteristics, fluctuating randomly around zero with values ranging approximately between -5 and +4 throughout the entire series. The series shows moderate persistence but consistently reverts to its mean, displaying no trend or explosive behavior.

In contrast, the AR(2) model with phi = -0.8 and phi-two = 0.3 is clearly non-stationary and exhibits explosive behavior. The series remains relatively flat and stable for the first 60 observations but then displays increasingly large oscillations that grow exponentially in magnitude, reaching values beyond ±3000 by the end of the series. The negative phi coefficient creates alternating positive and negative swings, while the combination of parameters produces an unstable process that diverges from any fixed mean.

This comparison illustrates the critical importance of parameter selection in ARIMA models—the ARMA(1,1) parameters yield a well-behaved, predictable series suitable for forecasting, while the AR(2) parameters create an unstable system that would be inappropriate for practical applications.

Question 9.7 Consider aus_airpassengers, the total number of passengers (in millions) from Australian air carriers for the period 1970-2011.

library(fpp3)

# Load and view the data
aus_airpassengers
# Plot original data to understand it
aus_airpassengers %>%
  autoplot(Passengers) +
  labs(title = "Australian Air Passengers (1970-2011)",
       y = "Passengers (millions)",
       x = "Year")

#### Understanding the initial dataset

The data shows a clear upward trend indicating that it is non-stationary.The amount of international visitors to australia incresed by approximately 650% from 1970 to 2011.

Question 9.7a Use ARIMA() to find an appropriate ARIMA model.

# Fit automatic ARIMA model
fit_auto <- aus_airpassengers %>%
  model(auto = ARIMA(Passengers))

# Display the selected model
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
# Check residuals (should look like white noise)
fit_auto %>% gg_tsresiduals()
## Warning: `gg_tsresiduals()` was deprecated in feasts 0.4.2.
## ℹ Please use `ggtime::gg_tsresiduals()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

#Check of residuals for white noise
augment(fit_auto) %>%
  features(.innov, ljung_box, lag = 10)
# Plot forecasts for next 10 periods
fc_auto <- fit_auto %>%
  forecast(h = 10)

fc_auto %>%
  autoplot(aus_airpassengers) +
  labs(title = "Forecasts from Automatic ARIMA Model",
       y = "Passengers (millions)")

#### Question 9.7a Response

Question 9.7b Write the model in terms of the backshift operator (submitted in the Brightspace link).

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

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

# Fit ARIMA(0,1,0) with drift
fit_010_drift <- aus_airpassengers %>%
  model(arima010_drift = ARIMA(Passengers ~ 1 + pdq(0,1,0)))

# Report model
report(fit_010_drift)
## 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
# Forecast next 10 periods
fc_010_drift <- fit_010_drift %>%
  forecast(h = 10)

# Plot forecast from ARIMA(0,1,0) only
fc_010_drift %>%
  autoplot(aus_airpassengers) +
  labs(title = "ARIMA(0,1,0) with drift",
       y = "Passengers (millions)",
       x = "Year")

# Compare with part a (Automatic ARIMA)
fc_010_drift %>%
  autoplot(aus_airpassengers, level = NULL) +
  autolayer(fc_auto, level = NULL, color = "orange") +
  labs(title = "Comparison: ARIMA(0,1,0) with drift vs Automatic ARIMA",
       y = "Passengers (millions)",
       x = "Year",
       subtitle = "Blue = ARIMA(0,1,0) with drift, Red = Automatic ARIMA") +
  theme(legend.position = "bottom")

##### Question 9.7d Plot forecasts from an ARIMA(2,1,2) model with drift and compare these to parts a and c. Remove the constant and see what happens.

fit_212_drift <- aus_airpassengers %>%
  model(arima212_drift = ARIMA(Passengers ~ 1 + pdq(2,1,2)))

# Report model
report(fit_212_drift)
## 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
# Forecast
fc_212_drift <- fit_212_drift %>%
  forecast(h = 10)

# Compare all three models
fc_comparison <- bind_rows(
  fc_auto %>% mutate(Model = "Automatic ARIMA"),
  fc_010_drift %>% mutate(Model = "ARIMA(0,1,0) with drift"),
  fc_212_drift %>% mutate(Model = "ARIMA(2,1,2) with drift")
)

fc_comparison %>%
  autoplot(aus_airpassengers, level = NULL) +
  labs(title = "Comparison of Three ARIMA Models",
       y = "Passengers (millions)",
       color = "Model")

##### Question 9.7d removing the constsnt to see what happens

fit_212_no_drift <- aus_airpassengers %>%
  model(arima212_no_drift = ARIMA(Passengers ~ 0 + pdq(2,1,2)))
## Warning: 1 error encountered for arima212_no_drift
## [1] non-stationary AR part from CSS
# Report model
report(fit_212_no_drift)
## Series: Passengers 
## Model: NULL model 
## NULL model
# Forecast
fc_212_no_drift <- fit_212_no_drift %>%
  forecast(h = 10)

# Compare with and without drift
fc_212_compare <- bind_rows(
  fc_212_drift %>% mutate(Model = "ARIMA(2,1,2) WITH drift"),
  fc_212_no_drift %>% mutate(Model = "ARIMA(2,1,2) WITHOUT drift")
)

fc_212_compare %>%
  autoplot(aus_airpassengers, level = NULL) +
  labs(title = "ARIMA(2,1,2): With vs Without Drift",
       y = "Passengers (millions)",
       color = "Model")
## Warning: Removed 10 rows containing missing values or values outside the scale range
## (`geom_line()`).

Question 9.7e

fit_021_constant <- aus_airpassengers %>%
  model(arima021_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 model
report(fit_021_constant)
## 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
# Forecast
fc_021_constant <- fit_021_constant %>%
  forecast(h = 10)

# Plot forecast
fc_021_constant %>%
  autoplot(aus_airpassengers) +
  labs(title = "Forecasts from ARIMA(0,2,1) with constant",
       y = "Passengers (millions)",
       subtitle = "What happens with second-order differencing?")

#### Question 9.7e analysis

After plotting the model with a constant graph smoothens slightly but overall the data peaks at the same points.

Question 9.8

For the United States GDP series (from global_economy):

# Load required libraries
library(fpp3)
library(ggplot2)
library(patchwork)

# Load the data
us_gdp <- global_economy %>%
  filter(Country == "United States") %>%
  select(Year, GDP)

# Convert to tsibble if needed
us_gdp <- us_gdp %>%
  as_tsibble(index = Year)

Question 9.8a

cat("=== STEP 1: Box-Cox Transformation Analysis ===\n\n")
## === STEP 1: Box-Cox Transformation Analysis ===
# Visualize original data
p1 <- us_gdp %>%
  autoplot(GDP) +
  labs(title = "Original US GDP Series",
       y = "GDP (USD)", x = "Year") +
  theme_minimal()

# Find optimal lambda
lambda <- us_gdp %>%
  features(GDP, features = guerrero) %>%
  pull(lambda_guerrero)

cat(sprintf("Optimal lambda from Guerrero method: %.4f\n", lambda))
## Optimal lambda from Guerrero method: 0.2819
cat(sprintf("Interpretation: lambda ≈ %.2f suggests ", lambda))
## Interpretation: lambda ≈ 0.28 suggests
if (abs(lambda - 0) < 0.1) {
  cat("log transformation\n")
} else if (abs(lambda - 0.5) < 0.1) {
  cat("square root transformation\n")
} else if (abs(lambda - 1) < 0.1) {
  cat("no transformation needed\n")
} else {
  cat(sprintf("Box-Cox transformation with lambda = %.2f\n", lambda))
}
## Box-Cox transformation with lambda = 0.28
# Apply transformation
us_gdp <- us_gdp %>%
  mutate(GDP_transformed = box_cox(GDP, lambda))

# Visualize transformed data
p2 <- us_gdp %>%
  autoplot(GDP_transformed) +
  labs(title = sprintf("Box-Cox Transformed GDP (λ = %.3f)", lambda),
       y = "Transformed GDP", x = "Year") +
  theme_minimal()

# Compare variance stability
print(p1 / p2)

cat("\n")

Question 9.8b

cat("=== STEP 2: Automatic ARIMA Model Selection ===\n\n")
## === STEP 2: Automatic ARIMA Model Selection ===
# Fit automatic ARIMA
fit_auto <- us_gdp %>%
  model(auto_arima = ARIMA(box_cox(GDP, lambda)))

# Display model
cat("Automatically selected ARIMA model:\n")
## Automatically selected ARIMA model:
report(fit_auto)
## Series: GDP 
## Model: ARIMA(1,1,0) w/ drift 
## Transformation: box_cox(GDP, 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
cat("\n")
# Extract model specification
auto_spec <- fit_auto %>% 
  glance() %>%
  select(ar_roots, ma_roots)

Question 9.8c

Comparison of Alternative Models

# Fit multiple candidate models
fit_models <- us_gdp %>%
  model(
    auto_arima = ARIMA(box_cox(GDP, lambda)),
    arima_110 = ARIMA(box_cox(GDP, lambda) ~ pdq(1,1,0)),
    arima_011 = ARIMA(box_cox(GDP, lambda) ~ pdq(0,1,1)),
    arima_111 = ARIMA(box_cox(GDP, lambda) ~ pdq(1,1,1)),
    arima_210 = ARIMA(box_cox(GDP, lambda) ~ pdq(2,1,0)),
    arima_012 = ARIMA(box_cox(GDP, lambda) ~ pdq(0,1,2)),
    arima_211 = ARIMA(box_cox(GDP, lambda) ~ pdq(2,1,1)),
    arima_112 = ARIMA(box_cox(GDP, lambda) ~ pdq(1,1,2))
  )

# Compare models using information criteria
model_comparison <- fit_models %>%
  glance() %>%
  select(.model, AIC, AICc, BIC) %>%
  arrange(AICc)

cat("Model Comparison (ranked by AICc):\n")
## Model Comparison (ranked by AICc):
print(model_comparison, n = Inf)
## # A tibble: 8 × 4
##   .model       AIC  AICc   BIC
##   <chr>      <dbl> <dbl> <dbl>
## 1 auto_arima  657.  657.  663.
## 2 arima_110   657.  657.  663.
## 3 arima_011   659.  659.  665.
## 4 arima_111   659.  659.  667.
## 5 arima_210   659.  659.  667.
## 6 arima_012   659.  660.  667.
## 7 arima_112   660.  661.  670.
## 8 arima_211   660.  661.  671.
cat("\n")

Question 9.8d

Best Model Selection and Residual Diagnostics

# Select best model (lowest AICc)
best_model_name <- model_comparison$.model[1]
cat(sprintf("Selected best model: %s\n\n", best_model_name))
## Selected best model: auto_arima
# Extract best model for diagnostics
best_fit <- fit_models %>%
  select(all_of(best_model_name))

# Report on best model
cat("Best Model Details:\n")
## Best Model Details:
report(best_fit)
## Series: GDP 
## Model: ARIMA(1,1,0) w/ drift 
## Transformation: box_cox(GDP, 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
cat("\n")
# Get residuals
resids <- best_fit %>% residuals()

# Ljung-Box test
lb_test <- best_fit %>%
  augment() %>%
  features(.innov, ljung_box, lag = 10, dof = 0)

cat(sprintf("Ljung-Box test (lag=10): p-value = %.4f\n", lb_test$lb_pvalue))
## Ljung-Box test (lag=10): p-value = 0.9554
if (lb_test$lb_pvalue > 0.05) {
  cat("✓ Residuals appear to be white noise (p > 0.05)\n")
} else {
  cat("✗ Evidence of remaining autocorrelation (p < 0.05)\n")
}
## ✓ Residuals appear to be white noise (p > 0.05)
# Plot residual diagnostics
best_fit %>% gg_tsresiduals() +
  labs(title = sprintf("Residual Diagnostics: %s", best_model_name))

# Additional residual statistics
resid_stats <- resids %>%
  as_tibble() %>%
  summarise(
    Mean = mean(.resid, na.rm = TRUE),
    SD = sd(.resid, na.rm = TRUE),
    Min = min(.resid, na.rm = TRUE),
    Max = max(.resid, na.rm = TRUE)
  )

cat("\nResidual Statistics:\n")
## 
## Residual Statistics:
print(resid_stats)
## # A tibble: 1 × 4
##    Mean    SD   Min   Max
##   <dbl> <dbl> <dbl> <dbl>
## 1  1.54  72.7 -263.  138.
cat("2. MODEL SELECTION:\n")
## 2. MODEL SELECTION:
cat(sprintf("   - Best ARIMA model: %s\n", best_model_name))
##    - Best ARIMA model: auto_arima
cat("   - Selected based on minimum AICc criterion\n")
##    - Selected based on minimum AICc criterion
cat(sprintf("   - Ljung-Box test p-value: %.4f\n", lb_test$lb_pvalue))
##    - Ljung-Box test p-value: 0.9554
cat("   - Residuals show good properties for forecasting\n\n")
##    - Residuals show good properties for forecasting

Question 9.8e

Forecast Generation and Visualization

# Generate forecasts
forecast_horizon <- 10
fc_arima <- best_fit %>%
  forecast(h = forecast_horizon)

cat(sprintf("Generated %d-year ahead forecasts\n\n", forecast_horizon))
## Generated 10-year ahead forecasts
# Plot forecasts
p_forecast <- fc_arima %>%
  autoplot(us_gdp, level = c(80, 95)) +
  labs(
    title = sprintf("US GDP Forecasts: %s Model", best_model_name),
    subtitle = sprintf("Box-Cox transformation with λ = %.3f", lambda),
    y = "GDP (USD)",
    x = "Year"
  ) +
  theme_minimal() +
  theme(legend.position = "bottom")

print(p_forecast)

# Display point forecasts
cat("\nPoint Forecasts:\n")
## 
## Point Forecasts:
fc_table <- fc_arima %>%
  as_tibble() %>%
  select(Year, .mean) %>%
  mutate(.mean = format(.mean, big.mark = ",", scientific = FALSE))

print(fc_table, n = Inf)
## # A tibble: 10 × 2
##     Year .mean             
##    <dbl> <chr>             
##  1  2018 20,170,192,721,946
##  2  2019 20,968,649,069,689
##  3  2020 21,788,083,470,242
##  4  2021 22,629,555,325,733
##  5  2022 23,493,735,442,044
##  6  2023 24,381,134,356,366
##  7  2024 25,292,192,461,831
##  8  2025 26,227,318,951,165
##  9  2026 27,186,909,507,315
## 10  2027 28,171,354,538,189
cat("Key considerations:\n")
## Key considerations:
cat("1. Trend continuation: Do forecasts follow the historical trend?\n")
## 1. Trend continuation: Do forecasts follow the historical trend?
cat("2. Uncertainty intervals: Do they widen appropriately over time?\n")
## 2. Uncertainty intervals: Do they widen appropriately over time?
cat("3. Economic plausibility: Are growth rates realistic?\n")
## 3. Economic plausibility: Are growth rates realistic?
cat("4. Prediction intervals: Check 80% and 95% coverage\n\n")
## 4. Prediction intervals: Check 80% and 95% coverage
# Calculate implied growth rates
fc_growth <- fc_arima %>%
  as_tibble() %>%
  select(Year, .mean) %>%
  mutate(
    growth_rate = (.mean / lag(.mean) - 1) * 100,
    growth_rate = ifelse(row_number() == 1, NA, growth_rate)
  )

cat("Implied Annual Growth Rates (%):\n")
## Implied Annual Growth Rates (%):
print(fc_growth %>% filter(!is.na(growth_rate)), n = Inf)
## # A tibble: 9 × 3
##    Year   .mean growth_rate
##   <dbl>   <dbl>       <dbl>
## 1  2019 2.10e13        3.96
## 2  2020 2.18e13        3.91
## 3  2021 2.26e13        3.86
## 4  2022 2.35e13        3.82
## 5  2023 2.44e13        3.78
## 6  2024 2.53e13        3.74
## 7  2025 2.62e13        3.70
## 8  2026 2.72e13        3.66
## 9  2027 2.82e13        3.62
cat("\n")
cat("3. FORECAST QUALITY:\n")
## 3. FORECAST QUALITY:
cat("   - Forecasts follow historical trend patterns\n")
##    - Forecasts follow historical trend patterns
cat("   - Prediction intervals widen appropriately\n")
##    - Prediction intervals widen appropriately
cat("   - Growth rates appear economically plausible\n\n")
##    - Growth rates appear economically plausible

Question 9.8f

ETS Model Comparison (No Transformation

# Fit ETS model
fit_ets <- us_gdp %>%
  model(ets = ETS(GDP))

cat("ETS Model Details:\n")
## ETS Model Details:
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
cat("\n")
# Generate ETS forecasts
fc_ets <- fit_ets %>%
  forecast(h = forecast_horizon)

# Compare both models visually
p_comparison <- us_gdp %>%
  autoplot(GDP) +
  autolayer(fc_arima, series = "ARIMA", level = 95) +
  autolayer(fc_ets, series = "ETS", level = 95) +
  labs(
    title = "US GDP Forecasts: ARIMA vs ETS",
    subtitle = "95% prediction intervals shown",
    y = "GDP (USD)",
    x = "Year",
    colour = "Model"
  ) +
  theme_minimal() +
  theme(legend.position = "bottom")
## Warning in ggdist::geom_lineribbon(without(intvl_mapping, "colour_ramp"), :
## Ignoring unknown parameters: `series`
## Warning in geom_line(mapping = without(mapping, "shape"), data =
## unpack_data(object[single_row[["FALSE"]], : Ignoring unknown parameters:
## `series`
## Warning in ggdist::geom_lineribbon(without(intvl_mapping, "colour_ramp"), :
## Ignoring unknown parameters: `series`
## Warning in geom_line(mapping = without(mapping, "shape"), data =
## unpack_data(object[single_row[["FALSE"]], : Ignoring unknown parameters:
## `series`
## Scale for fill_ramp is already present.
## Adding another scale for fill_ramp, which will replace the existing scale.
print(p_comparison)
## Ignoring unknown labels:
## • colour : "Model"

# Accuracy comparison on training data
accuracy_train <- bind_rows(
  accuracy(best_fit) %>% mutate(Model = best_model_name),
  accuracy(fit_ets) %>% mutate(Model = "ETS")
) %>%
  select(Model, RMSE, MAE, MAPE, MASE)

cat("\n=== Training Set Accuracy Comparison ===\n")
## 
## === Training Set Accuracy Comparison ===
print(accuracy_train)
## # A tibble: 2 × 5
##   Model               RMSE           MAE  MAPE  MASE
##   <chr>              <dbl>         <dbl> <dbl> <dbl>
## 1 auto_arima 147232166284.  87539655369.  1.52 0.257
## 2 ETS        166906260796. 104653662551.  1.91 0.307
cat("\n")
# Compare forecasts numerically
cat("=== Point Forecast Comparison (First 5 years) ===\n")
## === Point Forecast Comparison (First 5 years) ===
fc_comparison <- bind_rows(
  fc_arima %>% as_tibble() %>% select(Year, .mean) %>% mutate(Model = "ARIMA"),
  fc_ets %>% as_tibble() %>% select(Year, .mean) %>% mutate(Model = "ETS")
) %>%
  pivot_wider(names_from = Model, values_from = .mean) %>%
  mutate(Difference = ARIMA - ETS,
         Pct_Diff = (Difference / ETS) * 100)

print(head(fc_comparison, 5))
## # A tibble: 5 × 5
##    Year   ARIMA     ETS    Difference Pct_Diff
##   <dbl>   <dbl>   <dbl>         <dbl>    <dbl>
## 1  2018 2.02e13 2.01e13 104467706074.    0.521
## 2  2019 2.10e13 2.07e13 227637138423.    1.10 
## 3  2020 2.18e13 2.14e13 371784623582.    1.74 
## 4  2021 2.26e13 2.21e13 537969563679.    2.44 
## 5  2022 2.35e13 2.28e13 726862764596.    3.19
cat("\n")
cat("4. ARIMA vs ETS:\n")
## 4. ARIMA vs ETS:
cat("   - Both models produce reasonable forecasts\n")
##    - Both models produce reasonable forecasts
cat("   - ARIMA with Box-Cox transformation captures trend better\n")
##    - ARIMA with Box-Cox transformation captures trend better
cat("   - ETS is simpler but may be less flexible for GDP data\n")
##    - ETS is simpler but may be less flexible for GDP data
cat("   - Choice depends on forecast horizon and interpretability needs\n\n")
##    - Choice depends on forecast horizon and interpretability needs