Introduction

The following RMD contains answers to exercises in chapter 9 of the Forecasting: Principles and Practice textbook for CUNY SPS DATA 624 Spring 2025 https://otexts.com/fpp3/arima-exercises.html. This chapter focuses on ARIMA models. The 9.11 Exercises answered here include 9.1, 9.2, 9.3, 9.5, 9.6, 9.7 and 9.8.

Load Libraries

library(fpp3)
library(tidyverse)
library(dplyr)
library(ggplot2)
library(tidyr)
library(urca)

Q1

Prompt: Figure 9.32 [see textbook] 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?

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

a.

The graphs all indicate that the data are white noise as the three of them all contain data within the significance limit boundaries (blue lines). The differences are in their distance from the 0.0 center. The left ACF is for 36 numbers and contains the highest distances. The middle graph has data for 360 numbers and has about 1/3 the variance of the left. The right ACF graph has 1,000 numbers and contains the least variance. This is because larger sample sizes for this data mean random variation decreases.

b.

The critical values are at different distances from the mean of zero across graphs because each graph shows data for different sample sizes. The formula for the critical values is +-1.96/rootT where T is the time series length. As the graphs of data increase in sample size, the T increases and the critical values shrink. The autocorrelations are different for the different sample sizes, with the higher sample sizes having less autocorrelation for large samples and values nearing zero for the white noise expected.

Q2

Prompt: 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.

# Save the data for gafa_stock
amzn_df <- gafa_stock |>
  filter(Symbol == "AMZN") |>
  select("Date", "Close")

# Graph the time series, acf, and pacf
amzn_df |>
  gg_tsdisplay(Close, plot_type='partial') +
  labs(title = "Amazon Closing Stock Prices")
## Warning: Provided data has an irregular interval, results should be treated with caution. Computing ACF by observation.
## Provided data has an irregular interval, results should be treated with caution. Computing ACF by observation.

A non-stationary dataset will have no trend, constant variance, and no seasonal or cyclical trends. In this data, there is an upward trend as seen on the first visualization, the time series. The lag[1] ACF plot shows that the r1 is showing an autocorrelation across the data that is not decreasing quickly. The PACF graph has a large spike at lag 1, which also shows the data is non-stationary. The data should be first order differenced (difference(1)). This will show the change between one observation and the next, effectively making the dataset stationary. The seasonal differencing is not necessary here, as the seasonal componenent is not one of the factors contributing to the stationary vs non-stationary issue.

Q3

Prompt: 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.

  2. Accommodation takings in the state of Tasmania from aus_accommodation.

  3. Monthly sales from souvenirs.

a.

# Save the data to turkish_df
turkish_df <- global_economy |>
  filter(Country == "Turkey") |>
  select("Year", "GDP")

# Visualize the data
turkish_plot <- turkish_df |>
  autoplot(GDP) +
  labs(title = 'Turkish GDP')
turkish_plot

The data shows an upward trend and different variances across time.

# Find the lambda
turkish_lambda <- turkish_df |>
  features(GDP, features = guerrero) |>
  pull(lambda_guerrero)

# Do the box-cox transformation with the calculated
trans_turkish_df <- turkish_df |>
  mutate(GDP_trans = box_cox(GDP, lambda = turkish_lambda)) |> 
  mutate(GDP_diff = difference(GDP_trans)) 

trans_turkish_df |>
  autoplot(GDP_diff) +
  labs(title = "Turkish GDP, Differenced and Transformed")
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_line()`).

The transformed data does look more centered around 0. Let’s visualize the before and after time series, acf, and pacf graphs.

# Show original plots
turkish_df |>
  gg_tsdisplay(plot_type='partial') +
  labs(title = "Original Turkish GDP")
## Plot variable not specified, automatically selected `y = GDP`

# Show transformed plots
trans_turkish_df |>
  gg_tsdisplay(difference(box_cox(GDP, turkish_lambda)), plot_type='partial') +
  labs(title = "Differenced Turkish GDP with lambda = 1.57")
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_point()`).

This data looks great transformed. Before the transformation, the data has an upward trend, the acf is out of the boundaries for more than half of the data, and the pacf is quite high at lag 1. After the transformation, the data has a horizontal trend, and the acf and pacf both stay within their boundaries. The first order differencing and box-cox transformation helped to turn the non-stationary series into a stationary one.

b.

# Save the data to tas_df
tas_df <- aus_accommodation |>
  filter(State == "Tasmania") |>
  select("Date", "Takings")

# Visualize the data
tas_plot <- tas_df |>
  autoplot(Takings) +
  labs(title = 'Tasmanian Takings')
tas_plot

The data has an upward trend and some serious seasonality by quarters. Variance is also increasing over time so a box-cox transformation will be used.

# Find the lambda
tas_lambda <- tas_df |>
  features(Takings, features = guerrero) |>
  pull(lambda_guerrero)

# Do the box-cox transformation with the calculated lambda
trans_tas_df <- tas_df |>
  mutate(Takings_trans = box_cox(Takings, lambda = tas_lambda)) |> 
  mutate(Takings_diff = difference(Takings_trans)) 

# Show the graph
trans_tas_df |>
  autoplot(Takings_diff) +
  labs(title = "Tasmanian Takings, Differenced and Transformed")
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_line()`).

There is still seasonality here. Lets difference by quarter.

# Do the box-cox transformation with the calculated lambda and quarterly lag
trans_tas_df <- tas_df |>
  mutate(Takings_trans = box_cox(Takings, lambda = tas_lambda)) |> 
  mutate(Takings_diff = difference(Takings_trans, lag = 4))

# Show the graph
trans_tas_df |>
  autoplot(Takings_diff) +
  labs(title = "Tasmanian Takings, Differenced and Transformed")
## Warning: Removed 4 rows containing missing values or values outside the scale range
## (`geom_line()`).

Now this data looks to be corrected for seasonality, variance, and trend. Let’s compare the before and afters.

# Show original plots
tas_df |>
  gg_tsdisplay(plot_type='partial') +
  labs(title = "Original Tasmanian Takings")
## Plot variable not specified, automatically selected `y = Takings`

# Show transformed plots
trans_tas_df |>
  gg_tsdisplay(difference(box_cox(Takings, tas_lambda), 4), plot_type='partial') +
  labs(title = "Differenced Tasmanian Takings with lambda = 0.002")
## Warning: Removed 4 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 4 rows containing missing values or values outside the scale range
## (`geom_point()`).

This data looks better transformed. Before the transformation, the data has an upward trend with seasonality, the acf is out of the boundaries for more than half of the data, and the pacf is out of bounds for 4 lags. After the transformation, the data has a horizontal trend, and the acf and pacf both stay within their boundaries more than before the transformations. The first order differencing and box-cox transformation helped to turn the non-stationary series into a stationary one.

c.

# Save the data to souv_df
souv_df <- souvenirs

# Visualize the data
souv_plot <- souv_df |>
  autoplot(Sales) +
  labs(title = 'Monthly Souvenir Sales')
souv_plot

The data shows a small upward trend, yearly seasonality, and increasing variances across time.

# Find the lambda
souv_lambda <- souv_df |>
  features(Sales, features = guerrero) |>
  pull(lambda_guerrero)

# Do the box-cox transformation with the calculated
trans_souv_df <- souv_df |>
  mutate(Sales_trans = box_cox(Sales, lambda = souv_lambda)) |> 
  mutate(Sales_diff = difference(Sales_trans, lag = 12)) 

# Display transformed data
trans_souv_df |>
  autoplot(Sales_diff) +
  labs(title = "Souvenir Sales, Differenced and Transformed")
## Warning: Removed 12 rows containing missing values or values outside the scale range
## (`geom_line()`).

The transformed data does look more centered around 0. Let’s visualize the before and after time series, acf, and pacf graphs.

# Show original plots
souv_df |>
  gg_tsdisplay(plot_type='partial') +
  labs(title = "Original Souvenir Sales")
## Plot variable not specified, automatically selected `y = Sales`

# Show transformed plots
trans_souv_df |>
  gg_tsdisplay(difference(box_cox(Sales, souv_lambda), 12), plot_type='partial') +
  labs(title = "Differenced Turkish GDP with lambda = 0.002")
## Warning: Removed 12 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 12 rows containing missing values or values outside the scale range
## (`geom_point()`).

This data looks better transformed. Before the transformation, the data has an upward trend with seasonality and the pacf is out of bounds for 3 lags. After the transformation, the data has a horizontal trend, and the acf and pacf both stay within their boundaries more than before the transformations. The first order differencing and box-cox transformation helped to turn the non-stationary series into a stationary one.

Q5

Prompt: 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.

# Code from 2.10 exercise 7
set.seed(64)
myseries <- aus_retail |>
  filter(`Series ID` == sample(aus_retail$`Series ID`,1))
head(myseries)
## # A tsibble: 6 x 5 [1M]
## # Key:       State, Industry [1]
##   State    Industry                         `Series ID`    Month Turnover
##   <chr>    <chr>                            <chr>          <mth>    <dbl>
## 1 Victoria Other specialised food retailing A3349799R   1982 Apr     34.9
## 2 Victoria Other specialised food retailing A3349799R   1982 May     34.6
## 3 Victoria Other specialised food retailing A3349799R   1982 Jun     34.6
## 4 Victoria Other specialised food retailing A3349799R   1982 Jul     35.2
## 5 Victoria Other specialised food retailing A3349799R   1982 Aug     33.8
## 6 Victoria Other specialised food retailing A3349799R   1982 Sep     35.4
# Visualize data
autoplot(myseries, Turnover)

The data has an upward trend, seasonality, and variance than increasese then decreases. That last one makes this a bit of a challenge. Let’s visualize the acf and pcf graphs.

# Visualize the acf and pcf data
myseries |>
  gg_tsdisplay(Turnover, plot_type='partial') +
  labs(title = "Australian Retail Turnover")

The new info from the acf graph shows we should do a box-cox transformation. A stationary lag plot will not have values close to 1 for the entire dataset.

# Find the lambda
myseries_lambda <- myseries |>
  features(Turnover, features = guerrero) |>
  pull(lambda_guerrero)

# Do the box-cox transformation with the calculated
trans_myseries <- myseries |>
  mutate(Turnover_trans = box_cox(Turnover, lambda = myseries_lambda)) |> 
  mutate(Turnover_diff = difference(Turnover_trans, lag = 12)) 

# Display transformed data
trans_myseries |>
  autoplot(Turnover_diff) +
  labs(title = "Australian Retail Turnover, Differenced and Transformed")
## Warning: Removed 12 rows containing missing values or values outside the scale range
## (`geom_line()`).

The transformed data does look more centered around 0. Let’s visualize the acf and pacf graphs.

# Show transformed plots
trans_myseries |>
  gg_tsdisplay(difference(box_cox(Turnover, myseries_lambda), 12), plot_type='partial') +
  labs(title = "Differenced Australian Retail Turnover with lambda = -0.355")
## Warning: Removed 12 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 12 rows containing missing values or values outside the scale range
## (`geom_point()`).

It looks like this could use another differencing to make it stationary.

# Do the box-cox transformation with a second differencing and calculated lambda
trans_myseries <- myseries |>
  mutate(Turnover_trans = box_cox(Turnover, lambda = myseries_lambda)) |> 
  mutate(Turnover_diff1 = difference(Turnover_trans, lag = 12)) |>
  mutate(Turnover_diff2 = difference(Turnover_diff1, lag = 1))  # Second differencing

# Display transformed data
trans_myseries |>
  gg_tsdisplay(difference(difference(box_cox(Turnover, myseries_lambda), 12), 1), 
               plot_type = 'partial') +
  labs(title = "Twice Differenced Australian Retail Turnover with lambda = -0.355")
## Warning: Removed 13 rows 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_point()`).

Alright, that brought on its own issues. The lag being so far out of bounds at 12 for acf and pacf indicate that there is an issue with the seasonality adjustments. Let’s go with the original transformation. Now, compare with the original.

# Visualize the original
myseries |>
  gg_tsdisplay(Turnover, plot_type='partial') +
  labs(title = "Australian Retail Turnover")

# Visualize the transformation
trans_myseries |>
  gg_tsdisplay(difference(box_cox(Turnover, myseries_lambda), 12), plot_type='partial') +
  labs(title = "Differenced Australian Retail Turnover with lambda = -0.355")
## Warning: Removed 12 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 12 rows containing missing values or values outside the scale range
## (`geom_point()`).

This data looks better transformed. Before the transformation, the data has an upward trend with seasonality, the acf is out of bounds for all lags, and the pacf is out of bounds for 7 lags. After the transformation, the data has a horizontal trend, and the acf stays within their boundaries more than before the transformations. The first order differencing and box-cox transformation helped to turn the non-stationary series into a stationary one.

Q6

Prompt: Simulate and plot some data from simple ARIMA models.

a.

Use the following R code to generate data from an AR(1) model with phi1 = 0.6 and sigma2 = 1. The process starts with y1 = 0.

# Given code
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)

b.

Produce a time plot for the series. How does the plot change as you change phi1?

# Get phi1 = 1 data
for(i in 2:100)
  y[i] <- 1*y[i-1] + e[i]

# Visualize phi1 = 1
tsibble(idx = seq_len(100), y = y, index = idx) |>
  autoplot(y) +
  labs(title = "AR(1) model with phi1 = 1")

# Get phi1 = 0.9 data
for(i in 2:100)
  y[i] <- 0.9*y[i-1] + e[i]

# Visualize phi1 = 1
tsibble(idx = seq_len(100), y = y, index = idx) |>
  autoplot(y) +
  labs(title = "AR(1) model with phi1 = 0.9")

# Plot the AR(1) model
sim |>
  autoplot(y) +
  labs(title = "AR(1) model with phi1 = 0.6")

# Get phi1 = 0 data
for(i in 2:100)
  y[i] <- 0*y[i-1] + e[i]

# Visualize phi1 = 0
tsibble(idx = seq_len(100), y = y, index = idx) |>
  autoplot(y) +
  labs(title = "AR(1) model with phi1 = 0")

# Get phi1 = -0.6 data
for(i in 2:100)
  y[i] <- -0.6*y[i-1] + e[i]

# Visualize phi1 = -0.6
tsibble(idx = seq_len(100), y = y, index = idx) |>
  autoplot(y) +
  labs(title = "AR(1) model with phi1 = -0.6")

# Get phi1 = -0.9 data
for(i in 2:100)
  y[i] <- -0.9*y[i-1] + e[i]

# Visualize phi1 = -0.9
tsibble(idx = seq_len(100), y = y, index = idx) |>
  autoplot(y) +
  labs(title = "AR(1) model with phi1 = -0.9")

As phi1 decreases, so the variance increases. The oscillations get much closer together for the negative values of phi1 compared to their positive counterparts.

c. 

Write your own code to generate data from an MA(1) model with OMEGA1 = 0.6 and sigma2 = 1.

# Generate data from an MA(1) model with specified omega1 and sigma2
for(i in 2:100)
  y[i] <- e[i] + 0.6 * e[i-1]

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

d. 

Produce a time plot for the series. How does the plot change as you change THETA1?

# Get theta1 = 1 data
for(i in 2:100)
  y[i] <- e[i] + 1 * e[i-1]

# Visualize theta1 = 1
tsibble(idx = seq_len(100), y = y, index = idx) |>
  autoplot(y) +
  labs(title = "MA(1) model with theta1 = 1")

# Get theta1 = 0.9 data
for(i in 2:100)
  y[i] <- e[i] + 0.9 * e[i-1]

# Visualize theta1 = 0.9
tsibble(idx = seq_len(100), y = y, index = idx) |>
  autoplot(y) +
  labs(title = "MA(1) model with theta1 = 0.9")

# Plot the AR(1) model
sim |>
  autoplot(y) +
  labs(title = "MA(1) model with theta1 = 0.6")

# Get theta1 = 0 data
for(i in 2:100)
  y[i] <- e[i] + 0 * e[i-1]

# Visualize theta1 = 0
tsibble(idx = seq_len(100), y = y, index = idx) |>
  autoplot(y) +
  labs(title = "MA(1) model with theta1 = 0")

# Get theta1 = -0.6 data
for(i in 2:100)
  y[i] <- e[i] + -0.6 * e[i-1]

# Visualize theta1 = -0.6
tsibble(idx = seq_len(100), y = y, index = idx) |>
  autoplot(y) +
  labs(title = "MA(1) model with theta1 = -0.6")

# Get theta1 = -0.9 data
for(i in 2:100)
  y[i] <- e[i] + -0.9 * e[i-1]

# Visualize theta1 = -0.9
tsibble(idx = seq_len(100), y = y, index = idx) |>
  autoplot(y) +
  labs(title = "MA(1) model with theta1 = -0.9")

Like when phi increases, increased thetas mean less condensed oscillations of values. Also, there seems to be a decrease in total range as theta increases.

e.

Generate data from an ARMA(1,1) model with phi1 = 0.6, THETA1 = 0.6 and sigma2 = 1.

# Simulate numeric data for y
set.seed(64)
y <- numeric(100)

# Using phi1 = 0.6, THETA1 = 0.6 and sigma2 = 1, generate ARMA data
for(i in 2:100)
  y[i] <- e[i] + 0.6 * e[i-1] + 0.6 * y[i-1]

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

f. 

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

# Simulate numeric data for y
set.seed(64)
y <- numeric(100)

# Using phi1 = -0.8, phi2 = 0.3 and sigma2 = 1, generate AR(2) data
for(i in 3:100)
  y[i] <- -0.8 * y[i-1] + 0.3 * y[i-2] + e[i]

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

g.

Graph the latter two series and compare them.

# Graph first series
arma1_data |>
  autoplot(y) +
  labs(title = "ARMA(1,1) with phi1 = 0.6, THETA1 = 0.6 and sigma2 = 1")

# Graph second series
ar2_data |>
  autoplot(y) +
  labs(title = "AR(2) with phi1 = -0.8, phi2 = 0.3 and sigma2 = 1")

The ARMA(1,1) process has values that go up and down but stay close to zero, changing in a random way without getting too extreme. On the other hand, the AR(2) process keeps growing bigger and bouncing up and down more and more. This happens because of the way it’s set up, making it unstable, while the ARMA(1,1) process stays more steady.

Q7

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

a.

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.

# Visualize data
aus_airpassengers |>
  autoplot(Passengers) +
  labs(title = "Australian Air Passengers")

# Fit a model to find the values
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
# Show the model
fit |>
  gg_tsresiduals() +
  labs(title = "Residuals for ARIMA(0,2,1)")

# Forecast 
fit |>
  forecast(h = 10) |>
  autoplot(aus_airpassengers, level = NULL) +
  labs(title = "ARIMA(0,2,1) Forecast for Australian Air Passengers")

The ARIMA(0,2,1) model that was chosen does a good job of showing the increase in Australian air passengers. The leftover differences look random, like white noise. The 10-period forecast predicts more growth, but it gets harder to know exactly what will happen as time goes on.

b.

Write the model in terms of the backshift operator.

(1-B)^2 y[subt] = (1 - 0.8963B)e[subt]

c. 

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

# Get ARIMA values for 0,1,0
fit_drift <- aus_airpassengers |>
  model(ARIMA(Passengers ~ 1 + pdq(0,1,0)))
report(fit_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
# Fit the forecast
fc_drift <- fit_drift |>
  forecast(h = 10)

# Plot
fc_drift |>
  autoplot(aus_airpassengers, level = NULL) +
  labs(title = "ARIMA(0,1,0) with Drift Forecast for Australian Air Passengers")

# Generate forecasts for both models
fc_arima <- forecast(fit, h = 10) |> mutate(Model = "ARIMA(0,2,1)")
fc_drift <- forecast(fit_drift, h = 10) |> mutate(Model = "ARIMA(0,1,0) with Drift")

# Combine the forecasts
fc_comparison <- bind_rows(fc_arima, fc_drift)

# Plot the forecast comparison
fc_comparison |> 
  autoplot(aus_airpassengers, level = NULL) +
  labs(title = "Forecast Comparison: ARIMA(0,2,1) vs ARIMA(0,1,0) with Drift")

Compared to a, the drift graph has a lower slope.

d. 

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 ARIMA(2,1,2) with drift
fit_arima212_drift <- aus_airpassengers |>
  model(ARIMA(Passengers ~ 1 + pdq(2,1,2)))

# Forecast for the next 10 periods
fc_arima212_drift <- fit_arima212_drift |>
  forecast(h = 10)

# Plot
fc_arima212_drift |>
  autoplot(aus_airpassengers, level = NULL) +
  labs(title = "ARIMA(2,1,2) with Drift")

# Generate forecasts for each model
fc_arima_021 <- forecast(fit, h = 10) |>
  mutate(Model = "ARIMA(0,2,1)")
fc_arima_010_drift <- forecast(fit_drift, h = 10) |> 
  mutate(Model = "ARIMA(0,1,0) with Drift")
fc_arima_212_drift <- forecast(fit_arima212_drift, h = 10) |> 
  mutate(Model = "ARIMA(2,1,2) with Drift")

# Combine the forecasts
fc_comparison <- bind_rows(fc_arima_021, fc_arima_010_drift, fc_arima_212_drift)

# Plot the forecast comparison
fc_comparison |> 
  autoplot(aus_airpassengers, level = NULL) +
  labs(title = "Forecast Comparison")

# Take out the constant
fc_arima_no_const <- aus_airpassengers |>
  model(ARIMA(Passengers ~ 0 + pdq(2,1,2)))
## Warning: 1 error encountered for ARIMA(Passengers ~ 0 + pdq(2, 1, 2))
## [1] non-stationary AR part from CSS
# Show report for no constant
report(fc_arima_no_const)
## Series: Passengers 
## Model: NULL model 
## NULL model

The side-by-side comparison helps us see how each model predicts the future. The red and green lines, ARIMA(0,2,1) and ARIMA(0,1,0) with drift respectively, look almost the same. The blue line, ARIMA(2,1,2) with drift, has a much higher slope than the other two.

Without the constant, the model fails.

e.

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

# Add a constatn to the ARIMA(0,2,1)
fit_arima021_constant <- aus_airpassengers |>
  model(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.
# See the report
report(fit_arima021_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_arima021_constant <- fit_arima021_constant |>
  forecast(h = 10)

# Plot the forecast
fc_arima021_constant |>
  autoplot(aus_airpassengers, level = NULL) +
  labs(title = "ARIMA(0,2,1) with Constant Forecast")

# Forecast much more
fc_more_arima021_constant <- fit_arima021_constant |>
  forecast(h = 100)

# Plot the longer forecast
fc_more_arima021_constant |>
  autoplot(aus_airpassengers, level = NULL) +
  labs(title = "ARIMA(0,2,1) with Constant Longer Forecast")

The forecast shows a quadratic increase. In the shorter forecast, h=10 and it is difficult to see that the graph does not show a linear forecast. The second graph has h=100 and makes it easy to see the exponential increase in the forecast.

Q8

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

a.

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

us_gdp <- global_economy |>
  filter(Code == 'USA')

us_gdp |>
  gg_tsdisplay(GDP, plot_type='partial') +
  labs(title = "United States GDP")

# Find the lambda
us_gdp_lambda <- us_gdp |>
  features(GDP, features = guerrero) |>
  pull(lambda_guerrero)

# Visualize
us_gdp |>
  autoplot(box_cox(GDP, us_gdp_lambda))

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

# Fit the model
fit_us_gdp <- us_gdp |>
  model(arima = ARIMA(box_cox(GDP, us_gdp_lambda), stepwise = FALSE, approx = FALSE))

# Show the report 
report(fit_us_gdp)
## Series: GDP 
## Model: ARIMA(1,1,0) w/ drift 
## Transformation: box_cox(GDP, us_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

c. 

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

# Fit multiple models 
us_gdp_fit_multiple <- us_gdp |>
  model(arima110 = ARIMA(box_cox(GDP,us_gdp_lambda) ~ pdq(1,1,0)),
        arima111 = ARIMA(box_cox(GDP,us_gdp_lambda) ~ pdq(1,1,1)),
        arima112 = ARIMA(box_cox(GDP,us_gdp_lambda) ~ pdq(1,1,2)),
        arima113 = ARIMA(box_cox(GDP,us_gdp_lambda) ~ pdq(1,1,3)),
        arima120 = ARIMA(box_cox(GDP,us_gdp_lambda) ~ pdq(1,2,0)),
        arima121 = ARIMA(box_cox(GDP,us_gdp_lambda) ~ pdq(1,2,1)),
        arima210 = ARIMA(box_cox(GDP,us_gdp_lambda) ~ pdq(2,1,0)),
        arima212 = ARIMA(box_cox(GDP,us_gdp_lambda) ~ pdq(2,1,2)))

glance(us_gdp_fit_multiple) |>
  arrange(AICc) |>
  select(.model:BIC)
## # A tibble: 8 × 6
##   .model   sigma2 log_lik   AIC  AICc   BIC
##   <chr>     <dbl>   <dbl> <dbl> <dbl> <dbl>
## 1 arima121  5761.   -321.  648.  649.  655.
## 2 arima120  6780.   -326.  656.  656.  660.
## 3 arima110  5479.   -325.  657.  657.  663.
## 4 arima111  5580.   -325.  659.  659.  667.
## 5 arima210  5580.   -325.  659.  659.  667.
## 6 arima112  5630.   -325.  660.  661.  670.
## 7 arima113  5732.   -325.  662.  664.  674.
## 8 arima212  5734.   -325.  662.  664.  674.

d. 

choose what you think is the best model and check the residual diagnostics;

# The lowest AICc value is ARIMA(1,2,0)
us_gdp_fit_multiple |>
  select(arima120) |>
  gg_tsresiduals() +
  ggtitle("Residuals for ARIMA(1,2,0)")

e.

produce forecasts of your fitted model. Do the forecasts look reasonable?

us_gdp_fit_multiple |>
  forecast(h=30) |>
  filter(.model=='arima120') |>
  autoplot(us_gdp, level = NULL) +
  labs(title = "ARIMA(1,1,0) forecast")

f. 

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

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

fit_ets |>
  forecast(h=30) |>
  autoplot(us_gdp, level = NULL) +
  labs(title = "ETS Forecast")

Both the ARIMA(1,1,0) and ETS models show the GDP going up. But the ETS model has a linear prediction where the ARIMA(1,1,0) has a quadratic forecast.