Homework6

Do the exercises 9.1, 9.2, 9.3, 9.5, 9.6, 9.7, 9.8 in Hyndman. Please submit both the Rpubs link as well as your .rmd file.

Problem 9.1

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

  1. Explain the differences among these figures. Do they all indicate that the data are white noise?
library(jpeg)

img <- readJPEG("Pic91.jpg")
plot(1:2, type='n')  
rasterImage(img, 1, 1, 2, 2)

Each plot shows the autocorrelation function (ACF) for a different sample size of random numbers. In the first plot, which uses only 36 random numbers, the ACF values jump around quite a bit. That’s expected because with such a small sample, there’s more random variation. In the second plot, with 360 numbers, the ACF looks more stable, though there are still a few small fluctuations. By the time we get to the third plot, with 1,000 random numbers, the ACF is almost completely flat every value is close to zero across all lags.

This happens because as the sample size increases, the random ups and downs start to cancel out, and the ACF gives a clearer picture of what’s really going on. In this case, all three plots represent white noise. The only reason they look different is because of sampling variability the randomness just evens out more as we collect more data.

  1. 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?

The critical values sit at different distances from zero because they depend on the sample size. The 95% confidence limits for the ACF are roughly ±1.96/√n. With fewer observations, √n is smaller, so the limits are wider; with more observations, they’re narrower. For example: with n=36 the limits are about ±0.33, with n=360 about ±0.10, and with n=1,000 about ±0.06.

  • n=36 => √n = 6 => ±1.96/6 = 0.327

  • n=360 => √n = 18.97 => ±1.96/18.97 = 0.103

  • n=1000 => √n = 31.62 => ±1.96/31.62 = 0.062

Even though all three plots show white noise where the true autocorrelations are zero; the sample ACFs look a little different because of random variation. With smaller samples, those random ups and downs are more noticeable, and some bars might even fall outside the confidence limits just by chance. When we have a larger sample, those random effects tend to cancel out, and the ACF looks flatter and closer to zero.

So, all three figures represent white noise; the differences just come from the sample size. Smaller samples look “noisier” while larger ones give a cleaner, more stable picture.

Problem 9.2

A 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.

I will start by loading the data and filtering for Amazon stock.

amazon_stock <- gafa_stock |> 
  filter(Symbol == "AMZN")

glimpse(amazon_stock)
## Rows: 1,258
## Columns: 8
## Key: Symbol [1]
## $ Symbol    <chr> "AMZN", "AMZN", "AMZN", "AMZN", "AMZN", "AMZN", "AMZN", "AMZ…
## $ Date      <date> 2014-01-02, 2014-01-03, 2014-01-06, 2014-01-07, 2014-01-08,…
## $ Open      <dbl> 398.80, 398.29, 395.85, 395.04, 398.47, 403.71, 402.53, 397.…
## $ High      <dbl> 399.36, 402.71, 397.00, 398.47, 403.00, 406.89, 403.76, 399.…
## $ Low       <dbl> 394.02, 396.22, 388.42, 394.29, 396.04, 398.44, 393.80, 388.…
## $ Close     <dbl> 397.97, 396.44, 393.63, 398.03, 401.92, 401.01, 397.66, 390.…
## $ Adj_Close <dbl> 397.97, 396.44, 393.63, 398.03, 401.92, 401.01, 397.66, 390.…
## $ Volume    <dbl> 2137800, 2210200, 3170600, 1916000, 2316500, 2103000, 267950…

I will plot the daily closing prices for Amazon stock.

amazon_stock |> 
  autoplot(Close) + labs(title = "Amazon Daily Closing Prices")

The plot of Amazon’s daily closing prices shows a clear upward trend over time, indicating that the series is non-stationary. A stationary series would have a constant mean and variance over time, but here we see that the prices are increasing for later have a reduction, which suggests that the mean is changing.

Next, I will plot the ACF and PACF for the Amazon stock closing prices.

amazon_stock |> 
  ACF(Close) |> 
  autoplot()+ labs(title = "ACF of Amazon Daily Closing Prices")

In the ACF plot, we see that the autocorrelations are very high at large lags and decrease slowly over time. This slow decay is a hallmark of non stationarity, as it indicates that past values have a long lasting influence on future values. In a stationary series, we would expect the autocorrelations to drop off more quickly.

amazon_stock |> 
  PACF(Close) |> 
  autoplot()+ labs(title = "PACF of Amazon Daily Closing Prices")

The PACF plot shows significant spikes at the first lag, which also suggests non stationarity.

The plot of Amazon’s daily closing prices clearly shows an upward trend over time. The values keep increasing rather than fluctuating around a constant average, which tells us the series isn’t stationary ; its mean changes over time.

In the ACF plot, the correlations stay positive for many lags and decline very slowly. That slow decay means today’s price is still highly related to prices from quite a while ago, another sign of non-stationarity.

The PACF plot shows one big spike at lag 1 and then drops off quickly, which is typical of a random walk pattern. Together, these plots make it clear that Amazon’s stock prices are non-stationary and should be differenced (usually once) to make the series stationary before building a forecasting model.

Problem 9.3

For the following series, find an appropriate Box-Cox transformation and order of differencing in order to obtain stationary data.

  1. Turkish GDP from global_economy.
turkish_gdp <- global_economy |> 
  filter(Country == "Turkey") |> 
  select(Year, GDP)
glimpse(turkish_gdp)
## Rows: 58
## Columns: 2
## $ Year <dbl> 1960, 1961, 1962, 1963, 1964, 1965, 1966, 1967, 1968, 1969, 1970,…
## $ GDP  <dbl> 13995067818, 8022222222, 8922222222, 10355555556, 11177777778, 11…

Plot the Turkish GDP data to visualize it.

turkish_gdp |> 
  autoplot(GDP) + labs(title = "Turkish GDP Over Time",
       y = "Gross Domestic Product",
       x = "Year")

The plot of Turkish GDP shows a clear upward trend over time, indicating that the series is non-stationary. To stabilize the variance, I will apply a Box-Cox transformation. The Box-Cox transformation can be found using the Guerrero method.

lambda <- turkish_gdp |> 
  features(GDP, features = guerrero) |> 
  pull(lambda_guerrero)

lambda
## [1] 0.1572187

The optimal Box–Cox transformation parameter (λ) is approximately 0.16 for the Turkish GDP data. This suggests that applying a Box–Cox transformation with this λ value will help stabilize the variance.

Note: When λ is close to 0, a logarithmic transformation (log GDP) is appropriate; when λ is close to 1, no transformation is needed.

Next, Box–Cox transformation will be applied.

turkish_gdp_transformed <- turkish_gdp |> 
  mutate(GDP_boxcox = (GDP^lambda - 1) / lambda)

turkish_gdp_transformed |>
  autoplot(GDP_boxcox) + labs(title = "Box-Cox Transformed Turkish GDP Over Time",
       y = "Box-Cox Transformed GDP",
       x = "Year")

The plot of the Box-Cox transformed Turkish GDP still shows an upward trend, indicating that the series is still non-stationary. I will apply first differencing to the transformed data.

turkish_gdp_diff <- turkish_gdp_transformed |> 
  mutate(GDP_diff = difference(GDP_boxcox))

autoplot(turkish_gdp_diff, GDP_diff) +
  labs(title = "First Difference of Box-Cox Transformed Turkish GDP",
       y = "Differenced GDP",
       x = "Year")

In the plot of the first difference of the Box-Cox transformed Turkish GDP, the series appears to fluctuate around a constant mean, suggesting that it is now stationary.

In conclusion Even after applying the Box Cox transformation (λ = 0.16), the Turkish GDP series still shows an upward trend, indicating that it remains non-stationary. The transformation stabilizes the variance but does not remove the trend in the data. Therefore, I applied first differencing to the transformed series. The differenced plot shows that the long term trend has been removed, and the data now fluctuate around a roughly constant mean, suggesting that the series is closer to being stationary.

  1. Accommodation takings in the state of Tasmania from aus_accommodation.
tas_accom <- aus_accommodation |> 
  filter(State == "Tasmania")

autoplot(tas_accom, Takings) +
  labs(title = "Accommodation Takings in Tasmania",
       y = "Takings ($ million)",
       x = "Quarter")

The plot of accommodation takings in Tasmania shows a clear upward trend and seasonal patterns, indicating that the series is non-stationary. To stabilize the variance, I will apply a Box-Cox transformation. The Box-Cox transformation can be found using the Guerrero method.

lambda_tas <- tas_accom |> 
  features(Takings, features = guerrero) |> 
  pull(lambda_guerrero)

lambda_tas
## [1] 0.001819643

The optimal Box–Cox transformation parameter (λ) for accommodation takings in Tasmania is approximately 0.0018, which is effectively zero. Therefore, a logarithmic transformation is appropriate to stabilize the variance.

Note: - λ ≈ 1 -> No transformation needed - λ between 0.3–0.7 -> Mild power transform - λ ≈ 0 -> Log transform

Next, I will apply the logarithmic transformation.

tas_accom <- tas_accom |> 
  mutate(Takings_log = log(Takings))

tas_accom |> 
  mutate(Takings_log_diff = difference(Takings_log, 4)) |> 
  autoplot(Takings_log_diff)

The plot of the log-transformed and seasonally differenced series shows that the large upward trend from the original data has been removed. The values now fluctuate around a fairly constant level, and the variance appears much more stable over time.

The plot of accommodation takings in Tasmania shows a strong upward trend and clear seasonal patterns, indicating that the series is non-stationary. To stabilize the variance, a Box–Cox transformation was applied using the Guerrero method. The optimal λ value was approximately 0.0018, which is effectively zero. This means a logarithmic transformation is the most suitable choice. After applying the log transform, the data were seasonally differenced (lag = 4, because the data are quarterly) to remove seasonality and then first-differenced to eliminate the remaining trend. The resulting series fluctuates around a constant mean with roughly constant variance, showing that it is now stationary.

  1. Monthly sales from souvenirs.
souvenirs <- souvenirs

glimpse(souvenirs)
## Rows: 84
## Columns: 2
## $ Month <mth> 1987 Jan, 1987 Feb, 1987 Mar, 1987 Apr, 1987 May, 1987 Jun, 1987…
## $ Sales <dbl> 1664.81, 2397.53, 2840.71, 3547.29, 3752.96, 3714.74, 4349.61, 3…
autoplot(souvenirs, Sales) +
  labs(title = "Monthly Souvenir Sales",
       y = "Sales ($)", x = "Year")

This plot of monthly souvenir sales shows a clear upward trend and seasonal patterns, indicating that the series is non-stationary. To stabilize the variance, I will apply a Box-Cox transformation. The Box-Cox transformation can be found using the Guerrero method.

lambda_souvenirs <- souvenirs |> 
  features(Sales, features = guerrero) |> 
  pull(lambda_guerrero)

lambda_souvenirs
## [1] 0.002118221

The optimal Box–Cox transformation parameter (λ) for monthly souvenir sales is approximately 0.0021, which is effectively zero. Therefore, a logarithmic transformation is appropriate to stabilize the variance. I will create two new columns , Sales_log_diff12 and Sales_log_diff12_1, to remove seasonality and trend.

souvenirs <- souvenirs |>
  mutate(Sales_log = log(Sales))

souvenirs <- souvenirs |>
  mutate(
    Sales_log_diff12 = difference(Sales_log, 12), 
    Sales_log_diff12_1 = difference(Sales_log_diff12, 1) 
  )

autoplot(souvenirs, Sales_log_diff12_1) +
  labs(title = "Seasonally and First Differenced Log Sales",
       y = "Differenced Log Sales", x = "Year")

The plot of the log-transformed, seasonally differenced, and first-differenced series shows that the upward trend and seasonal patterns present in the original data have been removed. The values now fluctuate around a fairly constant mean, and the variance appears stable over time, indicating that the series is now stationary.

Problem 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.

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


autoplot(myseries, Turnover) +
  labs(title = "Retail Turnover Time Series",
       y = "Turnover ($ million)", x = "Year")

The plot of retail turnover shows a clear upward trend and seasonal patterns, indicating that the series is non-stationary with strong seasonality.

Now I will apply a Box-Cox transformation to stabilize the variance. The Box-Cox transformation can be found using the Guerrero method.

lambda_retail <- myseries |>
  features(Turnover, features = guerrero) |>
  pull(lambda_guerrero)

lambda_retail
## [1] 0.08303631

The optimal Box–Cox transformation parameter (λ) for retail turnover is approximately 0.08, which is close to zero. Therefore, a logarithmic transformation is appropriate to stabilize the variance.

After stabilizing variance with the Box–Cox transformation, the next goal is to make the mean constant over time; that is what differencing does.

  • The number of differences required to make the series stationary.
unitroot_nsdiffs(myseries$Turnover_log)
## nsdiffs 
##       0

nsdiffs = 0 means no seasonal differencing is needed.

  • Non-seasonal differencing order
unitroot_ndiffs(myseries$Turnover_log)
## ndiffs 
##   -Inf

ndiffs = 1 means one regular (first) difference to remove the overall trend and make the mean constant.

myseries <- myseries |>
  mutate(
    Turnover_log = log(Turnover),
    Turnover_log_diff = difference(Turnover_log, 1)
  )

autoplot(myseries, Turnover_log_diff) +
  labs(title = "Log-Transformed and First-Differenced Retail Turnover",
       y = "Differenced Log Turnover", x = "Year")

After applying a log transformation and one first difference, the retail turnover series no longer shows the upward trend seen in the original data. The values now move up and down around a fairly constant average, and the size of those fluctuations looks steady over time. This means that a single regular difference was enough to make the series stationary.

Problem 9.6

Simulate and plot some data from simple ARIMA models.

  1. Use the following R code to generate data from an AR(1) model with ϕ1=0.6 and σ2=1.The process starts with y1=0.
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)

The script provided is creating 100 simulated observations from an AR(1) model where the autoregressive coefficient (ϕ1) is 0.6 and the variance of the white noise error term (σ2) is 1. Each subsequent value is generated based on the previous value plus a random error term drawn from a standard normal distribution.

b.Produce a time plot for the series. How does the plot change as you change ϕ1 ?

sim |> 
  autoplot(y) +
  labs(title = "Simulated AR(1) Process (ϕ1=0.6, σ2=1)",
       y = "y", x = "Time Index")

The plot of the simulated AR(1) process with ϕ1=0.6 and σ2=1 shows the series fluctuates around zero but tends to remember where it’s been for a short while. When one value increases, the next few values usually stay a bit higher before settling back down. This kind of short-term persistence is typical of an AR(1) model with ϕ = 0.6, where each observation is moderately influenced by the one before it.

Now I want to validate the autocorrelation

sim |> ACF(y) |> autoplot() +
  labs(title = "ACF of Simulated AR(1) Series")

The ACF plot shows a gradual decline in autocorrelation values as the lag increases, which is typical of an AR(1) process. The autocorrelation at lag 1 is positive and fairly strong, reflecting how each observation is influenced by the previous one. As the lag increases, the correlations steadily weaken, showing that the effect of past values fades over time. This pattern is consistent with the behavior expected from an AR(1) model with ϕ = 0.6.

c.Write your own code to generate data from an MA(1) model with ϕ1=0.6 and σ2=1

The data were generated using an MA(1) model defined by yt = et + θ1 * et-1, where et is white noise with mean 0 and variance σ2.

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

d.Produce a time plot for the series. How does the plot change as you change θ1 ?

sim_ma |> 
  autoplot(y) +
  labs(title = "Simulated MA(1) Process (θ1=0.6, σ2=1)",
       y = "y", x = "Time Index")

The plot of the simulated MA(1) process with θ₁ = 0.6 and σ² = 1 shows a series that fluctuates around zero but alternates more sharply than the AR(1) process. In this model, each value depends on both the current and previous random shocks. Because θ₁ is positive, when one value rises, the next often moves in the opposite direction, creating quick reversals between highs and lows.

Now I will validate the autocorrelation

sim_ma |> ACF(y) |> autoplot() +
  labs(title = "ACF of Simulated MA(1) Process (θ = 0.6)")

The ACF plot for the simulated MA(1) process with θ₁ = 0.6 shows a distinct pattern where the autocorrelation is significant at lag 1 and then drops to zero for all subsequent lags. This is characteristic of an MA(1) model, where only the first lag has a direct influence on the current value.

e.Generate data from an ARMA(1,1) model with ϕ1=0.6 , θ1=0.6 and σ2=1.

In this case, I will use ARMA(1,1) model defined by yt = ϕ1 * yt-1 + et + θ1 * et-1, where et is white noise with mean 0 and variance σ2.

y_arma <- numeric(100)
e_arma <- rnorm(100)
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 = "Simulated ARMA(1,1) Process (ϕ1=0.6, θ1=0.6, σ2=1)",
       y = "y", x = "Time Index")

The plot of the simulated ARMA(1,1) process with ϕ₁ = 0.6, θ₁ = 0.6, and σ² = 1 shows a series that blends the behavior of both AR and MA models. The values move up and down around zero, showing some short term persistence where high values tend to follow other highs. At the same time, sudden shifts appear because of random shocks from the MA component. Altogether, the series looks smoother than a pure MA(1) process but less steady than a simple AR(1), capturing both momentum and quick reversals.

f.Generate data from an AR(2) model with ϕ1=−0.8 ,ϕ2=0.3 and σ2=1. (Note that these parameters will give a non-stationary series.)

y_ar2 <- numeric(100)
e_ar2 <- rnorm(100)
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 = "Simulated AR(2) Process (ϕ1=-0.8, ϕ2=0.3, σ2=1)",
       y = "y", x = "Time Index")

The simulated AR(2) process with ϕ₁ = −0.8, ϕ₂ = 0.3, and σ² = 1 shows a series that moves up and down quite erratically without settling around a consistent average. Instead of staying stable, the values tend to drift and sometimes grow larger over time. This happens because the chosen parameters make the process non-stationary ; the effects of past values don’t fade quickly, causing the series to wander rather than return to a steady level.

g.Graph the latter two series and compare them.

sim_arma_plot <- sim_arma |> 
  autoplot(y) +
  labs(title = "Simulated ARMA(1,1) Process (ϕ1=0.6, θ1=0.6, σ2=1)",
       y = "y", x = "Time Index") +
  theme(plot.title = element_text(size = 9))  

sim_arma2_plot <- sim_ar2 |> 
  autoplot(y) +
  labs(title = "Simulated AR(2) Process (ϕ1=-0.8, ϕ2=0.3, σ2=1)",
       y = "y", x = "Time Index") +
  theme(plot.title = element_text(size = 9))  

grid.arrange(sim_arma_plot, sim_arma2_plot, ncol = 2, widths = c(1, 1))

The two plots side by side show clear differences in behavior between the ARMA(1,1) and AR(2) processes. The ARMA(1,1) series fluctuates around zero with a mix of short-term persistence and sudden shifts, reflecting both momentum and random shocks. In contrast, the AR(2) series appears more erratic and tends to drift over time without settling around a consistent average. This wandering behavior is due to the non-stationary nature of the AR(2) model with the given parameters, causing it to lack a stable mean. Overall, the ARMA(1,1) process looks more controlled and stable compared to the unpredictable and drifting AR(2) process.

Problem 9.7

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

  1. Use ARIMA() to find an appropriate ARIMA model. What model was selected. Check that the residuals look like white noise. Plot forecasts for the next 10 periods.

I will start running a time series plot

aus_airpassengers |> 
  autoplot(Passengers) +
  labs(title = "Australian Air Passengers (1970-2011)",
       y = "Passengers (millions)", x = "Year")

This plot of Australian air passengers from 1970 to 2011 shows a clear upward trend. The number of passengers increases over time, indicating growth in air travel.

Next, I will use the ARIMA() function to find an appropriate ARIMA model for the data.

fit <- aus_airpassengers |> 
  model(ARIMA(Passengers))
report(fit)
## 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

The function automatically finds an appropriate BOX-COX transformation if needed, it identifies the best differencing order to make the series stationary, and selects the AR and MA orders based on information criteria like AICc.

Model selected ARIMA(0,2,1) , which means: - No AR terms (p=0) - Two differences (d=2) to remove trend and make the series stationary - One MA term (q=1) to capture short-term correlations in the differenced series - No seasonal part, the model found that explicit seasonal differencing wasn’t necessary

Now I will check the residuals to see if they look like white noise.

gg_tsresiduals(fit)

augment(fit) |> 
  features(.innov, ljung_box, lag = 24, dof = 1)
## # A tibble: 1 × 3
##   .model            lb_stat lb_pvalue
##   <chr>               <dbl>     <dbl>
## 1 ARIMA(Passengers)    24.0     0.405

The residual diagnostics show no clear pattern or seasonality, and the Ljung–Box test (p = 0.405) confirms that the residuals are uncorrelated, indicating that the model provides a good fit.

Overall, the residuals behave like white noise, and the model is appropriate for forecasting future passenger numbers.

Now I will plot forecasts for the next 10 periods.

fit |> 
  forecast(h = 10) |> 
  autoplot(aus_airpassengers) +
  labs(title = "10-Period Forecast of Australian Air Passengers",
       y = "Passengers (millions)", x = "Year")

The time series plot of Australian air passengers from 1970 to 2011 shows a clear upward trend, reflecting steady growth in air travel over time. Using the ARIMA() function, the best-fitting model selected was ARIMA(0,2,1), which means two regular differences were needed to remove the strong trend, along with one moving average term to account for short-term fluctuations.

The residual diagnostics show that the residuals fluctuate randomly around zero, with no visible trend or seasonality. The Ljung–Box test (p = 0.405) confirms that the residuals behave like white noise, indicating the model fits the data well.

Forecasts for the next 10 years suggest that passenger numbers will continue to rise steadily. The prediction intervals widen over time, showing greater uncertainty for longer-term forecasts, but the overall trend points toward continued growth in Australian air travel.

  1. Write the model in terms of the backshift operator.

Backshift operator as also called the lag operator, denoted by B, is a mathematical tool used in time series analysis to simplify the representation of autoregressive (AR) and moving average (MA) models. It shifts a time series back by one period. For example, Byt = yt-1.

Previous fitting model was : - -0.8963 is the estimated MA(1) coefficient - d=2, two regular differences - p=0, no AR terms

Now I will write the model in terms of the backshift operator.

(1 - B)^2 yt = (1 - 0.8963B) et

This equation represents the ARIMA(0,2,1) model, where (1−𝐵)^2 indicates that the series has been differenced twice to achieve stationarity, and (1−0.8963B) represents the moving-average component with coefficient −0.8963.

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

In this case we will compare a random walk (ARIMA(0,1,0)) with drift to the ARIMA(0,2,1) model from part a.

fc_a <- fit |> 
  forecast(h = 10)

fit_rw <- aus_airpassengers |> 
  model(rw_drift = RW(Passengers ~ drift()))

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



autoplot(aus_airpassengers, level = NULL) +
  autolayer(fc_a, color = "blue",
            alpha = 0.3,  size = 1.1, 
            linetype = "solid", level = c(80,95)) +
  autolayer(fc_rw, color = "orange",
            alpha = 0.3,  size = 1.1, 
            linetype = "dashed", level = c(80,95)) +
  labs(
    title = "Forecast Comparison: ARIMA(0,2,1) vs ARIMA(0,1,0) with Drift",
    y = "Passengers (millions)",
    x = "Year",
    caption = "Blue: ARIMA(0,2,1)  |  Orange: ARIMA(0,1,0) with drift"
  ) +
  theme_minimal(base_size = 13)

The forecast comparison shows that both the ARIMA(0,2,1) model and the ARIMA(0,1,0) model with drift predict a steady increase in Australian air passengers over the next 10 years. The ARIMA(0,2,1) model (blue line) gives slightly higher forecasts and wider intervals, suggesting it captures more of the short-term ups and downs in the data. Meanwhile, the random walk with drift (orange line) produces a simpler, smoother trend that reflects a consistent average growth rate.

  1. 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.
library(fpp3)

# Fit ALL models together
fits <- aus_airpassengers |>
  model(
    arima_021  = ARIMA(Passengers),                 # part (a)
    rw_drift   = RW(Passengers ~ drift()),          # part (c)
    arima212_d = ARIMA(Passengers ~ 1 + pdq(2,1,2)),# (2,1,2) with drift
    arima212_0 = ARIMA(Passengers ~ 0 + pdq(2,1,2)) # (2,1,2) no constant
  )

# Forecast 10 periods
fc <- forecast(fits, h = 10)

# Nice, uncluttered plot: single 95% ribbon, transparent, clear legend
autoplot(fc, aus_airpassengers, level = 95) +
  scale_colour_manual(
    values = c(
      arima_021  = "#276EF1",   # blue
      rw_drift   = "#FF8C00",   # orange
      arima212_d = "#22A06B",   # green
      arima212_0 = "#8E44AD"    # purple
    ),
    name = "Model"
  ) +
  scale_fill_manual(
    values = c(
      arima_021  = "#276EF133",
      rw_drift   = "#FF8C0033",
      arima212_d = "#22A06B33",
      arima212_0 = "#8E44AD33"
    ),
    guide = "none"
  ) +
  guides(colour = guide_legend(override.aes = list(size = 1.2))) +
  labs(
    title   = "Forecast Comparison: ARIMA(0,2,1) vs RW+Drift vs ARIMA(2,1,2) (± constant)",
    y = "Passengers (millions)", x = "Year",
    caption = "Blue: ARIMA(0,2,1)  •  Orange: RW (0,1,0)+drift  •  Green: ARIMA(2,1,2)+drift  •  Purple: ARIMA(2,1,2) no constant"
  ) +
  theme_minimal(base_size = 13) +
  theme(plot.caption = element_text(margin = margin(t = 8)))

The ARIMA(2,1,2) model with drift (green) produces a similar upward trend to the ARIMA(0,2,1) and random walk with drift models but does a better job capturing short-term ups and downs in the data, thanks to its additional AR and MA components. When the drift term is removed (purple), the forecasts become noticeably flatter and stay closer to recent observations, since the model no longer assumes a steady long-term increase. Overall, the ARIMA(2,1,2) with drift strikes a good balance—capturing short-term variation while still reflecting the long-term growth in air travel—whereas the version without the constant is too conservative and tends to underestimate future passenger numbers.

  1. Plot forecasts from an ARIMA(0,2,1) model with a constant. What happens?

The model ARIMA(0,2,1) has p=0 (no auto regressive terms), d=2 (the data is differenced twice to remove trend), q=1 (one moving average term), and adding a constant introduces a drift term in the second differenced series.

Adding a constant to the ARIMA(0,2,1) model introduces a drift term that allows for ongoing growth in the forecasts. Instead of flattening out, the predictions now continue upward, reflecting a steady long term increase in passenger numbers.

fit_021_c <- aus_airpassengers |> 
  model(arima_021_c = ARIMA(Passengers ~ 1 + pdq(0,2,1)))
fc_021_c <- fit_021_c |> forecast(h = 10)

autoplot(fc_021_c, aus_airpassengers, level = 95) +
  labs(
    title = "Forecast from ARIMA(0,2,1) with Constant",
    y = "Passengers (millions)", x = "Year"
  ) +
  theme_minimal(base_size = 13)

The forecast from the ARIMA(0,2,1) model with a constant shows a steady upward trend in passenger numbers over the next 10 years. Including the constant term adds a drift component, allowing the model to represent consistent long-term growth rather than leveling off. As a result, the forecasts are higher compared to the version without a constant. The prediction intervals widen gradually, reflecting growing uncertainty over time. Overall, adding the constant enables the model to better capture the long-term upward trend seen in the historical data.

Problem 9.8

For the United States GDP series (from global_economy):

  1. if necessary, find a suitable Box-Cox transformation for the data;

I will start by loading the dataset

us_gdp <- global_economy |> 
  filter(Country == "United States") |> 
  select(Year, GDP)

glimpse(us_gdp)
## Rows: 58
## Columns: 2
## $ Year <dbl> 1960, 1961, 1962, 1963, 1964, 1965, 1966, 1967, 1968, 1969, 1970,…
## $ GDP  <dbl> 5.433000e+11, 5.633000e+11, 6.051000e+11, 6.386000e+11, 6.858000e…

I will plot the US GDP data to visualize it.

us_gdp |> 
  autoplot(GDP) + labs(title = "US GDP Over Time",
       y = "USD (current, billions)",
       x = "Year")

The plot of US GDP shows a clear upward trend over time, indicating that the series is non-stationary. To stabilize the variance, I will apply a Box-Cox transformation. The Box-Cox transformation can be found using the Guerrero method.

lambda_us <- us_gdp |> 
  features(GDP, features = guerrero) |> 
  pull(lambda_guerrero)
lambda_us
## [1] 0.2819443

The optimal Box–Cox transformation parameter (λ) is approximately 0.28 for the US GDP data. This suggests that applying a Box–Cox transformation with this λ value will help stabilize the variance.

us_gdp_transformed <- us_gdp |> 
  mutate(GDP_boxcox = box_cox(GDP, lambda_us))

us_gdp_transformed|>
  autoplot(GDP_boxcox) + labs(title = "Box-Cox Transformed US GDP Over Time",
       y = "Box-Cox Transformed GDP",
       x = "Year")

The United States GDP series shows a clear upward trend and increasing variability over time. Using the Guerrero method, the optimal Box–Cox transformation parameter was estimated as λ = 0.2819. Because this value is between 0 and 1, a Box–Cox transformation with this λ was applied to stabilize the variance. The transformed series displays more uniform variation over time and is ready for differencing to achieve stationarity.

  1. fit a suitable ARIMA model to the transformed data using ARIMA();

In this part I will use the ARIMA() function to fit a suitable ARIMA model to the Box-Cox transformed US GDP data.

fit_us <- us_gdp_transformed |> 
  model(ARIMA(GDP_boxcox))
report(fit_us)
## Series: GDP_boxcox 
## 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
gg_tsresiduals(fit_us)

augment(fit_us) |>
  features(.innov, ljung_box, lag = 24, dof = 3)
## # A tibble: 1 × 3
##   .model            lb_stat lb_pvalue
##   <chr>               <dbl>     <dbl>
## 1 ARIMA(GDP_boxcox)    10.4     0.973

After applying a Box–Cox transformation (λ = 0.2819) to stabilize the variance, the automatic ARIMA() function selected an ARIMA(1,1,0) model with drift for the transformed US GDP series. This model includes one autoregressive term and one order of differencing to remove the trend, along with a constant (drift) term representing steady long-term growth in GDP.

The residual diagnostics show that the model fits the data well. The Ljung–Box test statistic (10.43, p = 0.97) indicates no significant autocorrelation, meaning the residuals behave like white noise. Overall, the ARIMA(1,1,0) with drift model provides an excellent fit for the Box–Cox–transformed US GDP series and successfully captures its underlying dynamics.

Now I will visualize the forecasts

fit_us |> 
  forecast(h = 10) |> 
  autoplot(us_gdp_transformed) +
  labs(title = "Forecasts of US GDP (Box–Cox transformed)",
       y = "Transformed GDP", x = "Year")

The forecasts from the ARIMA(1,1,0) with drift model on the Box–Cox transformed US GDP series show a continued upward trend over the next 10 years. The prediction intervals widen over time, reflecting increasing uncertainty in longer-term forecasts. Overall, the model captures the steady growth pattern observed in the historical data.

  1. try some other plausible models by experimenting with the orders chosen;

For this part, I will experiment with several plausible ARIMA models by varying the orders of the autoregressive (p) and moving average (q) components while keeping the differencing order (d) fixed at 1, as suggested by the initial model. I will compare the models using the AICc criterion to identify the bestfitting one.

fit_us_alt <- us_gdp_transformed |>
  model(
    arima_011 = ARIMA(GDP_boxcox ~ pdq(0,1,1) + PDQ(0,0,0)),
    arima_111 = ARIMA(GDP_boxcox ~ pdq(1,1,1) + PDQ(0,0,0)),
    arima_210 = ARIMA(GDP_boxcox ~ pdq(2,1,0) + PDQ(0,0,0)),
    arima_211 = ARIMA(GDP_boxcox ~ pdq(2,1,1) + PDQ(0,0,0))
  ) 

fit_us_alt|> 
  glance() |> 
  arrange(AICc)
## # A tibble: 4 × 8
##   .model    sigma2 log_lik   AIC  AICc   BIC ar_roots  ma_roots 
##   <chr>      <dbl>   <dbl> <dbl> <dbl> <dbl> <list>    <list>   
## 1 arima_011  5689.   -326.  659.  659.  665. <cpl [0]> <cpl [1]>
## 2 arima_111  5580.   -325.  659.  659.  667. <cpl [1]> <cpl [1]>
## 3 arima_210  5580.   -325.  659.  659.  667. <cpl [2]> <cpl [0]>
## 4 arima_211  5647.   -325.  660.  661.  671. <cpl [2]> <cpl [1]>

To explore alternative specifications, several ARIMA models were tested on the Box–Cox–transformed US GDP data, including ARIMA(0,1,1), ARIMA(1,1,1), ARIMA(2,1,0), and ARIMA(2,1,1). Model comparison using AICc and BIC values showed that all models performed similarly, with AICc values ranging between 659 and 661.

The ARIMA(0,1,1) model had the lowest AICc (659.18), but the improvement over the simpler ARIMA(1,1,0) with drift model from part (b) was negligible. This confirms that the original model remains the most appropriate, offering a good balance between simplicity and accuracy while adequately capturing the dynamics of US GDP growth.

  1. choose what you think is the best model and check the residual diagnostics;
gg_tsresiduals(fit_us)

augment(fit_us) |>
  features(.innov, ljung_box, lag = 24, dof = 3)
## # A tibble: 1 × 3
##   .model            lb_stat lb_pvalue
##   <chr>               <dbl>     <dbl>
## 1 ARIMA(GDP_boxcox)    10.4     0.973

The residual diagnostics for the ARIMA(1,1,0) with drift model show that the residuals fluctuate randomly around zero, with no visible trend or autocorrelation. The ACF plot confirms that all autocorrelations lie within the 95% confidence bounds, and the histogram of residuals appears approximately normal.

The Ljung–Box test (p = 0.97) provides further evidence that the residuals are uncorrelated and behave like white noise. This indicates that the model successfully captures the main dynamics of the U.S. GDP series and is well-suited for forecasting. Based on the AICc and BIC comparisons, the ARIMA(1,1,0) with drift model was selected as the best fit for the Box–Cox–transformed U.S. GDP data, offering a good balance between simplicity and accuracy.

  1. produce forecasts of your fitted model. Do the forecasts look reasonable?
fc_us <- fit_us |> 
  forecast(h = 10)

autoplot(fc_us, us_gdp_transformed, level = 95) +
  labs(
    title = "10-Year Forecast of US GDP (Box–Cox Transformed)",
    y = "Transformed GDP", x = "Year"
  ) +
  theme_minimal(base_size = 13)

The forecasts from the ARIMA(1,1,0) with drift model show a continued upward trend in U.S. GDP over the next 10 years, consistent with the long-term growth observed in the historical data. The prediction intervals widen slightly as the forecast horizon increases, indicating greater uncertainty in longer-term projections.

  1. compare the results with what you would obtain using ETS() (with no transformation).

In this part I will use ETS() to fit an exponential smoothing state space model to the original (untransformed) US GDP data.

fit_ets_us <- us_gdp |> 
  model(ETS(GDP))

report(fit_ets_us)
## 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

An ETS(M,A,N) model was fitted to the untransformed U.S. GDP data. This model includes multiplicative errors and an additive trend, which effectively captures the steady growth pattern in GDP. The smoothing parameters (α = 0.999, β = 0.501) indicate that the model gives strong weight to recent observations while maintaining a moderately adaptive trend.

Now I will check the residual diagnostics to see if they look like white noise.

gg_tsresiduals(fit_ets_us)

The residual diagnostics for the ETS(M,A,N) model indicate a good fit. The residuals fluctuate randomly around zero, with no noticeable trend or correlation over time. The autocorrelation function confirms that all correlations fall within the 95% confidence limits, and the residuals are approximately normally distributed.

fc_ets_us <- fit_ets_us |> forecast(h = 10)

autoplot(fc_ets_us, us_gdp) +
  labs(
    title = "Forecasts of U.S. GDP using ETS(M,A,N) Model",
    y = "GDP (Billions of USD)",
    x = "Year"
  ) +
  theme_minimal(base_size = 13)

The forecasts from the ETS model show a smooth, consistent upward trajectory in GDP over the next 10 years, similar to the forecasts produced by the ARIMA(1,1,0) with drift model. However, the ETS forecasts appear slightly smoother because the model directly estimates the level and trend, while the ARIMA model captures changes through differencing.

Overall, both models produce plausible and realistic forecasts. The ARIMA model provides a slightly more flexible fit by accounting for short-term autocorrelation, whereas the ETS(M,A,N) model offers a simpler, trend-driven view that emphasizes long-term growth stability.