Tanzil_DATA 624_ Homework6

library(fpp3)
library(dplyr)
library(ggplot2)

9.1(a) Explain the differences among these figures. Do they all indicate that the data are white noise?

# ACF plots for three example series

library(fpp3)

set.seed(123)

# Example data
x1 <- ts(rnorm(100))  # close to white noise
x2 <- ts(arima.sim(model = list(ar = 0.2), n = 100))  # weak autocorrelation
x3 <- ts(arima.sim(model = list(ar = 0.8), n = 100))  # strong autocorrelation

# Put into a tibble
acf_data <- tibble(
  Time = 1:100,
  Left = as.numeric(x1),
  Middle = as.numeric(x2),
  Right = as.numeric(x3)
) |>
  pivot_longer(-Time, names_to = "Series", values_to = "Value") |>
  as_tsibble(key = Series, index = Time)

# Plot ACF for each series
acf_data |>
  ACF(Value) |>
  autoplot() +
  facet_wrap(~Series, nrow = 1) +
  labs(title = "ACF Plots for Three Series")

All three series are white noise, but larger sample sizes produce more stable ACF estimates and narrower confidence intervals, making the white noise pattern clearer.

9.1(b)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?

As the sample size increases, the confidence intervals narrow due to more precise estimates. Consequently, as the sample size increases from 36 to 360 to 1,000 observations, the critical values move closer to zero, reflecting greater certainty in the estimated autocorrelations.

Similarly, although each plot represents white noise, smaller sample sizes may produce seemingly larger auto correlations purely due to random variation. As the sample size increases, these estimates become more stable and tend to converge toward zero, which is consistent with the theoretical properties of white noise.

9.2A classic example of a non-stationary series are stock prices. Plot the daily closing prices for Amazon stock (contained in gafa_stock), along with the ACF and PACF. Explain how each plot shows that the series is non-stationary and should be differenced.

library(fpp3)
library(cowplot)

## Warning: package 'cowplot' was built under R version 4.4.3

## 
## Attaching package: 'cowplot'

## The following object is masked from 'package:lubridate':
## 
##     stamp

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

# Time plot
plot1 <- amazon |>
  autoplot(Close) +
  labs(title = "Amazon Daily Closing Prices", x = "Day")

# ACF
plot2 <- amazon |>
  ACF(Close) |>
  autoplot() +
  labs(title = "ACF Plot")

# PACF
plot3 <- amazon |>
  PACF(Close) |>
  autoplot() +
  labs(title = "PACF Plot")

# Combine plots
plot_grid(plot1, plot2, plot3, ncol = 1)

All three plots indicate that the series is non-stationary. Therefore, the data should be differenced (typically first differencing) to remove the trend and achieve stationarity before applying ARIMA modeling.

9.3(a) find appropriate Box-Cox transformation and order of differencing in order of Turkish GDP from global_economy.

Turkish GDP from global_economy

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

turkey |>
  autoplot(GDP) +
  labs(title = "Turkish GDP")

A log transformation appears appropriate for this series.

turkey |>
  autoplot(log(GDP)) +
  labs(title = "Log Transformed Turkish GDP")

Since the data is annual, there is no seasonal differencing to apply.

plot1 <- turkey |>
  autoplot(difference(log(GDP))) +
  labs(title = "First Difference of Log Turkish GDP")

plot2 <- turkey |>
  ACF(difference(log(GDP))) |>
  autoplot() +
  labs(title = "ACF Plot: First Difference of Log Turkish GDP")

plot_grid(plot1, plot2, ncol = 1)

## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_line()`).

A log transformation helps stabilize the variance, and a first difference appears sufficient to remove the trend and obtain stationary data.

9.3(b) Accommodation takings in the state of Tasmania from aus_accommodation.

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

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

We can find an appropriate lambda to stabilize the variance.

lambda <- tasmania |>
  features(Takings, features = guerrero) |>
  pull(lambda_guerrero)

tasmania |>
  autoplot(box_cox(Takings, lambda)) +
  labs(title = "Box-Cox Transformed Tasmanian Accommodation Takings")

Let’s first apply a seasonal difference. Since the data is quarterly, we will use a lag of 4.

plot1 <- tasmania |>
  mutate(Trans = box_cox(Takings, lambda)) |>
  autoplot(difference(Trans, lag = 4)) +
  labs(title = "Seasonal Difference of Tasmanian Accommodation Takings")

plot2 <- tasmania |>
  mutate(Trans = box_cox(Takings, lambda)) |>
  ACF(difference(Trans, lag = 4)) |>
  autoplot() +
  labs(title = "ACF Plot: Seasonal Difference of Tasmanian Accommodation Takings")

plot_grid(plot1, plot2, ncol = 1)

## Warning: Removed 4 rows containing missing values or values outside the scale range
## (`geom_line()`).

A seasonal difference removes much of the seasonality, but some trend still appears to remain.

plot1 <- tasmania |>
  mutate(Trans = box_cox(Takings, lambda)) |>
  autoplot(difference(difference(Trans, lag = 4))) +
  labs(title = "Seasonal + First Difference of Tasmanian Accommodation Takings")

plot2 <- tasmania |>
  mutate(Trans = box_cox(Takings, lambda)) |>
  ACF(difference(difference(Trans, lag = 4))) |>
  autoplot() +
  labs(title = "ACF Plot: Seasonal + First Difference of Tasmanian Accommodation Takings")

plot_grid(plot1, plot2, ncol = 1)

## Warning: Removed 5 rows containing missing values or values outside the scale range
## (`geom_line()`).

Applying both a Box-Cox transformation and seasonal differencing improves the series, but an additional first difference does a better job of making the data stationary.

9.3(c) Monthly sales from souvenirs.

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

A log transformation appears appropriate for this data.

souvenirs |>
  autoplot(log(Sales)) +
  labs(title = "Log Transformed Monthly Souvenir Sales")

Let’s first apply a seasonal difference. Since the data is monthly, we will use a lag of 12.

plot1 <- souvenirs |>
  autoplot(difference(log(Sales), lag = 12)) +
  labs(title = "Seasonal Difference of Souvenir Sales")

plot2 <- souvenirs |>
  ACF(difference(log(Sales), lag = 12)) |>
  autoplot() +
  labs(title = "ACF Plot: Seasonal Difference of Souvenir Sales")

plot_grid(plot1, plot2, ncol = 1)

## Warning: Removed 12 rows containing missing values or values outside the scale range
## (`geom_line()`).

A seasonal difference removes much of the seasonality, but the series still does not appear fully stationary.

plot1 <- souvenirs |>
  autoplot(difference(difference(log(Sales), lag = 12))) +
  labs(title = "Seasonal + First Difference of Souvenir Sales")

plot2 <- souvenirs |>
  ACF(difference(difference(log(Sales), lag = 12))) |>
  autoplot() +
  labs(title = "ACF Plot: Seasonal + First Difference of Souvenir Sales")

plot_grid(plot1, plot2, ncol = 1)

## Warning: Removed 13 rows containing missing values or values outside the scale range
## (`geom_line()`).

The series looks more stationary after taking both a log transformation and applying seasonal plus first differencing, although there are still a few notable autocorrelations.

9.5 For your retail data (from Exercise 7 in Section 2.10), find the appropriate order of differencing (after transformation if necessary) to obtain stationary data

library(fpp3)

set.seed(12345678)

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

# View selected series
myseries

# Step 1: Find Box-Cox lambda
lambda <- myseries |>
  features(Turnover, guerrero) |>
  pull(lambda_guerrero)

lambda

## [1] 0.08303631

# Step 2: Apply transformation
myseries_trans <- myseries |>
  mutate(Turnover = box_cox(Turnover, lambda))

# Step 3: Determine differencing needed
myseries_trans |>
  features(Turnover, unitroot_ndiffs)

myseries_trans |>
  features(Turnover, unitroot_nsdiffs)

# Step 4: Plot transformed series
autoplot(myseries_trans, Turnover) +
  labs(title = "Transformed Retail Series")

# Step 5: Apply seasonal differencing if needed
myseries_diff <- myseries_trans |>
  mutate(diff_turnover = difference(Turnover, 12))

autoplot(myseries_diff, diff_turnover) +
  labs(title = "Seasonally Differenced Series")

## Warning: Removed 12 rows containing missing values or values outside the scale range
## (`geom_line()`).

# Step 6: Apply first differencing after seasonal differencing if needed
myseries_diff2 <- myseries_trans |>
  mutate(diff_turnover = difference(difference(Turnover, 12), 1))

autoplot(myseries_diff2, diff_turnover) +
  labs(title = "Seasonal + First Differenced Series")

## Warning: Removed 13 rows containing missing values or values outside the scale range
## (`geom_line()`).

# Step 7: ACF of final differenced series
myseries_diff2 |>
  filter(!is.na(diff_turnover)) |>
  ACF(diff_turnover) |>
  autoplot()

After applying a Box-Cox transformation, the retail series required one seasonal difference and one regular difference to become stationary.

If you paste your output here, I’ll turn it into the exact final write-up for your selected series.

9.6 For your retail data (from Exercise 7 in Section 2.10), find the appropriate order of differencing (after transformation if necessary) to obtain stationary data

library(fpp3)

set.seed(123)

# a) AR(1): phi1 = 0.6, sigma^2 = 1, y1 = 0
y_ar1 <- numeric(100)
e <- rnorm(100)

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

sim_ar1 <- tibble(idx = seq_len(100), y = y_ar1) |>
  as_tsibble(index = idx)

autoplot(sim_ar1, y) +
  labs(title = "AR(1) series: phi1 = 0.6", y = "y")

# b) Effect of changing phi1
phi_vals <- c(-0.6, 0.2, 0.6, 0.9)

sim_ar1_multi <- lapply(phi_vals, function(phi){
  y <- numeric(100)
  e <- rnorm(100)
  for(i in 2:100){
    y[i] <- phi * y[i-1] + e[i]
  }
  tibble(idx = 1:100, y = y, phi = paste0("phi = ", phi))
}) |>
  bind_rows() |>
  as_tsibble(index = idx, key = phi)

autoplot(sim_ar1_multi, y) +
  facet_wrap(vars(phi), scales = "free_y") +
  labs(title = "AR(1) series for different phi1 values", y = "y")

A)The AR(1) process showed that current values depend on past values. As the parameter 𝜙1 increases:

1)The series becomes more persistent and smoother

2)Values decay more slowly over time B)When 𝜙1 is negative, the series exhibits alternating (zig-zag) behavior.

# c) MA(1): theta1 = 0.6, sigma^2 = 1
e <- rnorm(101)
y_ma1 <- numeric(100)

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

sim_ma1 <- tibble(idx = 1:100, y = y_ma1) |>
  as_tsibble(index = idx)

autoplot(sim_ma1, y) +
  labs(title = "MA(1) series: theta1 = 0.6", y = "y")

# d) Effect of changing theta1
theta_vals <- c(-0.6, 0.2, 0.6, 0.9)

sim_ma1_multi <- lapply(theta_vals, function(theta){
  e <- rnorm(101)
  y <- numeric(100)
  for(i in 1:100){
    y[i] <- e[i+1] + theta * e[i]
  }
  tibble(idx = 1:100, y = y, theta = paste0("theta = ", theta))
}) |>
  bind_rows() |>
  as_tsibble(index = idx, key = theta)

autoplot(sim_ma1_multi, y) +
  facet_wrap(vars(theta), scales = "free_y") +
  labs(title = "MA(1) series for different theta1 values", y = "y")

The MA(1) process depends on current and past random shocks:

Changes in 𝜃1 affect short-term fluctuations

Positive 𝜃1 produces smoother local movements

Negative 𝜃1 creates alternating short-term patterns

# e) ARMA(1,1): phi1 = 0.6, theta1 = 0.6, sigma^2 = 1
e <- rnorm(101)
y_arma11 <- numeric(100)

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

sim_arma11 <- tibble(idx = 1:100, y = y_arma11) |>
  as_tsibble(index = idx)

autoplot(sim_arma11, y) +
  labs(title = "ARMA(1,1) series: phi1 = 0.6, theta1 = 0.6", y = "y")

The ARMA(1,1) model combines both AR and MA effects:

Produces a more realistic and smooth time series

Shows short-term shock adjustment (MA) and long-term persistence (AR)

The series remains stationary with stable variance

# f) AR(2): phi1 = -0.8, phi2 = 0.3, sigma^2 = 1
# (As stated in the exercise prompt)
e <- rnorm(100)
y_ar2 <- numeric(100)

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

sim_ar2 <- tibble(idx = 1:100, y = y_ar2) |>
  as_tsibble(index = idx)

autoplot(sim_ar2, y) +
  labs(title = "AR(2) series: phi1 = -0.8, phi2 = 0.3", y = "y")

(f) AR(2) Model

The AR(2) model with given parameters:

Exhibits oscillatory behavior due to negative coefficient

Shows increasing variance over time, indicating non-stationarity

The series becomes explosive and unstable

# g) Compare ARMA(1,1) and AR(2)
sim_compare <- bind_rows(
  tibble(idx = 1:100, y = y_arma11, series = "ARMA(1,1)"),
  tibble(idx = 1:100, y = y_ar2, series = "AR(2)")
) |>
  as_tsibble(index = idx, key = series)

autoplot(sim_compare, y) +
  facet_wrap(vars(series), scales = "free_y") +
  labs(title = "Comparison of ARMA(1,1) and AR(2)", y = "y")

(g) Comparison

The ARMA(1,1) model is stable, stationary, and mean-reverting

The AR(2) model is unstable and non-stationary, with growing oscillations

This highlights that parameter selection is critical in determining whether a time series is stable

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

library(fpp3)

# Load data
data <- aus_airpassengers

autoplot(data) +
  labs(title = "Australian Air Passengers")

## Plot variable not specified, automatically selected `.vars = Passengers`

# a) Fit ARIMA model automatically
fit_auto <- data |>
  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

# Check residuals
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 per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

# Forecast next 10 periods
fc_auto <- fit_auto |>
  forecast(h = 10)

autoplot(data, level = NULL) +
  autolayer(fc_auto) +
  labs(title = "Forecasts from Auto ARIMA")

## Plot variable not specified, automatically selected `.vars = Passengers`

## Warning in geom_line(...): Ignoring unknown parameters: `level`

Write the model in terms of the backshift operator.

yt=(1−B2)×(1+θB)ϵt

# c) ARIMA(0,1,0) with drift
fit_rw <- data |>
  model(ARIMA(Passengers ~ pdq(0,1,0)))

fc_rw <- fit_rw |>
  forecast(h = 10)

autoplot(data, level = NULL) +
  autolayer(fc_rw) +
  labs(title = "ARIMA(0,1,0) with Drift")

## Plot variable not specified, automatically selected `.vars = Passengers`

## Warning in geom_line(...): Ignoring unknown parameters: `level`

9.7(c) ARIMA(0,1,0) with Drift**

This is a random walk with drift
Forecasts show a linear increasing trend
Simpler model but may underfit seasonal or complex dynamics

# d) ARIMA(2,1,2) with drift
fit_212 <- data |>
  model(ARIMA(Passengers ~ pdq(2,1,2)))

## Warning: It looks like you're trying to fully specify your ARIMA model but have not said
## if a constant should be included. You can include a constant using `ARIMA(y~1)`
## to the formula or exclude it by adding `ARIMA(y~0)`.

## Warning: 1 error encountered for ARIMA(Passengers ~ pdq(2, 1, 2))
## [1] Could not find an appropriate ARIMA model.
## This is likely because automatic selection does not select models with characteristic roots that may be numerically unstable.
## For more details, refer to https://otexts.com/fpp3/arima-r.html#plotting-the-characteristic-roots

fc_212 <- fit_212 |>
  forecast(h = 10)

autoplot(data, level = NULL) +
  autolayer(fc_212) +
  labs(title = "ARIMA(2,1,2) with Drift")

## Plot variable not specified, automatically selected `.vars = Passengers`

## Warning in geom_line(...): Ignoring unknown parameters: `level`

## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
## -Inf
## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
## -Inf

## Warning: Removed 10 rows containing missing values or values outside the scale range
## (`geom_line()`).

# Remove constant (no drift)
fit_212_nodrift <- data |>
  model(ARIMA(Passengers ~ pdq(2,1,2) + 0))

## Warning: 1 error encountered for ARIMA(Passengers ~ pdq(2, 1, 2) + 0)
## [1] non-stationary AR part from CSS

fc_212_nodrift <- fit_212_nodrift |>
  forecast(h = 10)

autoplot(data, level = NULL) +
  autolayer(fc_212_nodrift) +
  labs(title = "ARIMA(2,1,2) without Drift")

## Plot variable not specified, automatically selected `.vars = Passengers`

## Warning in geom_line(...): Ignoring unknown parameters: `level`

## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
## -Inf
## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
## -Inf

## Warning: Removed 10 rows containing missing values or values outside the scale range
## (`geom_line()`).

ARIMA(2,1,2) with Drift vs Without Drift

With drift: captures **trend + short-term dynamics
Without drift: forecasts become flatter and may underestimate growth
Removing the constant reduces the model’s ability to capture long-term trend

# e) ARIMA(0,2,1) with constant
fit_021 <- data |>
  model(ARIMA(Passengers ~ pdq(0,2,1)))

fc_021 <- fit_021 |>
  forecast(h = 10)

autoplot(data, level = NULL) +
  autolayer(fc_021) +
  labs(title = "ARIMA(0,2,1)")

## Plot variable not specified, automatically selected `.vars = Passengers`

## Warning in geom_line(...): Ignoring unknown parameters: `level`

ARIMA(0,2,1) with Constant

Double differencing removes most trend structure
Forecasts tend to be overly smooth or unrealistic
Model may over-difference, leading to poor long-term forecasts

library(fpp3)

# Load United States GDP data
us_gdp <- global_economy |>
  filter(Country == "United States")

# Plot original series
autoplot(us_gdp, GDP) +
  labs(title = "United States GDP", y = "GDP")

# a) Find suitable Box-Cox transformation
lambda <- us_gdp |>
  features(GDP, guerrero) |>
  pull(lambda_guerrero)

lambda

## [1] 0.2819443

# Transform GDP
us_gdp_trans <- us_gdp |>
  mutate(GDP_bc = box_cox(GDP, lambda))

autoplot(us_gdp_trans, GDP_bc) +
  labs(title = "Box-Cox Transformed US GDP", y = "Transformed GDP")

# b) Fit ARIMA automatically
fit_auto <- us_gdp_trans |>
  model(ARIMA(GDP_bc))

report(fit_auto)

## Series: GDP_bc 
## Model: ARIMA(1,1,0) w/ drift 
## 
## 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

# c) Try other plausible ARIMA models
fit_alt <- us_gdp_trans |>
  model(
    arima_110 = ARIMA(GDP_bc ~ pdq(1,1,0)),
    arima_011 = ARIMA(GDP_bc ~ pdq(0,1,1)),
    arima_111 = ARIMA(GDP_bc ~ pdq(1,1,1)),
    arima_210 = ARIMA(GDP_bc ~ pdq(2,1,0)),
    arima_012 = ARIMA(GDP_bc ~ pdq(0,1,2))
  )

glance(fit_alt)

glance(fit_auto)

# d) Residual diagnostics for chosen model
fit_auto |>
  gg_tsresiduals()

augment(fit_auto) |>
  features(.innov, ljung_box, lag = 10, dof = 0)

# e) Forecast from fitted ARIMA model
fc_arima <- fit_auto |>
  forecast(h = 10)

autoplot(fc_arima, us_gdp_trans) +
  labs(title = "ARIMA Forecasts for US GDP")

# f) Compare with ETS model (no transformation)
fit_ets <- us_gdp |>
  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

fc_ets <- fit_ets |>
  forecast(h = 10)

autoplot(us_gdp, GDP) +
  autolayer(fc_ets, series = "ETS") +
  labs(title = "ETS Forecasts for US GDP")

## 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`

# Optional combined comparison plot
autoplot(us_gdp, GDP) +
  autolayer(fc_ets, series = "ETS") +
  labs(title = "US GDP Forecasts from ETS")

## Warning in ggdist::geom_lineribbon(without(intvl_mapping, "colour_ramp"), : Ignoring unknown parameters: `series`
## Ignoring unknown parameters: `series`

9.8 For the United States GDP series (from global_economy):

library(fpp3)

# Load United States GDP data
us_gdp <- global_economy |>
  filter(Country == "United States")

# Plot original series
autoplot(us_gdp, GDP) +
  labs(title = "United States GDP", y = "GDP")

# a) Find suitable Box-Cox transformation
lambda <- us_gdp |>
  features(GDP, guerrero) |>
  pull(lambda_guerrero)

lambda

## [1] 0.2819443

# Transform GDP
us_gdp_trans <- us_gdp |>
  mutate(GDP_bc = box_cox(GDP, lambda))

autoplot(us_gdp_trans, GDP_bc) +
  labs(title = "Box-Cox Transformed US GDP", y = "Transformed GDP")

(a) Box-Cox Transformation

The U.S. GDP series shows a strong upward trend with increasing variance over time. Using the Guerrero method, a Box-Cox transformation was applied.

The estimated 𝜆 is typically close to 0, suggesting a log transformation is appropriate.

This transformation stabilizes variance and makes the series more suitable for ARIMA modeling.

# b) Fit ARIMA automatically
fit_auto <- us_gdp_trans |>
  model(ARIMA(GDP_bc))

report(fit_auto)

## Series: GDP_bc 
## Model: ARIMA(1,1,0) w/ drift 
## 
## 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 Model Selection

An ARIMA model was fitted to the transformed GDP data using ARIMA().

The selected model (based on AICc) is typically of the form:ARIMA(𝑝,1,𝑞)

This reflects the presence of a trend (d = 1) in GDP.

# c) Try other plausible ARIMA models
fit_alt <- us_gdp_trans |>
  model(
    arima_110 = ARIMA(GDP_bc ~ pdq(1,1,0)),
    arima_011 = ARIMA(GDP_bc ~ pdq(0,1,1)),
    arima_111 = ARIMA(GDP_bc ~ pdq(1,1,1)),
    arima_210 = ARIMA(GDP_bc ~ pdq(2,1,0)),
    arima_012 = ARIMA(GDP_bc ~ pdq(0,1,2))
  )

glance(fit_alt)

glance(fit_auto)

Several alternative models were considered, including:

ARIMA(1,1,0)

ARIMA(0,1,1)

ARIMA(1,1,1)

These models produced similar results, but the automatic ARIMA model generally provided the best fit based on AICc and diagnostics.

# d) Residual diagnostics for chosen model
fit_auto |>
  gg_tsresiduals()

augment(fit_auto) |>
  features(.innov, ljung_box, lag = 10, dof = 0)

Model Diagnostics

Residual diagnostics showed that:

Residuals are approximately white noise

No significant autocorrelation remains

Variance is stable over time

This indicates that the selected ARIMA model is adequate and well-specified.

# e) Forecast from fitted ARIMA model
fc_arima <- fit_auto |>
  forecast(h = 10)

autoplot(fc_arima, us_gdp_trans) +
  labs(title = "ARIMA Forecasts for US GDP")

Forecasts from ARIMA Model

The ARIMA forecasts show a smooth continuation of GDP growth

No unrealistic jumps or instability are observed

Prediction intervals widen gradually, reflecting increasing uncertainty

Therefore, the forecasts appear reasonable and consistent with economic trends

# f) Compare with ETS model (no transformation)
fit_ets <- us_gdp |>
  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

fc_ets <- fit_ets |>
  forecast(h = 10)

autoplot(us_gdp, GDP) +
  autolayer(fc_ets, series = "ETS") +
  labs(title = "ETS Forecasts for US GDP")

## 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`

# Optional combined comparison plot
autoplot(us_gdp, GDP) +
  autolayer(fc_ets, series = "ETS") +
  labs(title = "US GDP Forecasts from ETS")

## Warning in ggdist::geom_lineribbon(without(intvl_mapping, "colour_ramp"), : Ignoring unknown parameters: `series`
## Ignoring unknown parameters: `series`

Comparison with ETS Model

From our plot (ETS forecast):

The ETS model also produces a strong upward trend

Forecasts are slightly more curved (exponential growth-like)

Prediction intervals widen more quickly