DATA 624 Homework 6

library(fpp2)

## Warning: package 'fpp2' was built under R version 4.5.2

## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo

## ── Attaching packages ──────────────────────────────────────────── fpp2 2.5.1 ──

## ✔ ggplot2   4.0.0      ✔ fma       2.5   
## ✔ forecast  8.24.0     ✔ expsmooth 2.3

##

library(fpp3)

## Registered S3 method overwritten by 'tsibble':
##   method               from 
##   as_tibble.grouped_df dplyr

## ── Attaching packages ──────────────────────────────────────────── fpp3 1.0.2 ──

## ✔ tibble      3.2.1     ✔ tsibble     1.1.6
## ✔ dplyr       1.1.4     ✔ tsibbledata 0.4.1
## ✔ tidyr       1.3.1     ✔ feasts      0.4.2
## ✔ lubridate   1.9.4     ✔ fable       0.5.0

## Warning: package 'fable' was built under R version 4.5.2

## ── Conflicts ───────────────────────────────────────────────── fpp3_conflicts ──
## ✖ lubridate::date()    masks base::date()
## ✖ dplyr::filter()      masks stats::filter()
## ✖ tsibble::intersect() masks base::intersect()
## ✖ tsibble::interval()  masks lubridate::interval()
## ✖ dplyr::lag()         masks stats::lag()
## ✖ tsibble::setdiff()   masks base::setdiff()
## ✖ tsibble::union()     masks base::union()

## 
## Attaching package: 'fpp3'

## The following object is masked from 'package:fpp2':
## 
##     insurance

library(forecast)
library(readxl)
library(fable)

Introduction

In this report, I will be answering questions 1, 2, 3, 5, 6, 7, and 8 in the 9th chapter from the Hyndman book.

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

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

The three ACF plots look different mainly because of the different sample sizes. The series with 36 observations looks much noisier, with larger mins and maxs, since small samples produce more random variation in the sample autocorrelations. As the sample size increases to 360 and then 1,000 observations, the ACF values stay closer to zero and the plots look more stable. Even though the plots look different, all three are consistent with white noise since none of them show a clear pattern (increasing, decreasing, or seasonal structure). In fact, the autocorrelations fluctuate randomly around zero, which is what is expected from a white noise process.

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 are at different distances from zero because they depend on the sample size. The approximate critical values for white noise are + or - 1.96/√T, where T is the number of observations. When the sample size is small, there is more sampling variability. As the sample size increases, the bounds become narrower because the estimates become more precise. Each plot comes from a different random sample, which makes the white noise representations to be different. In theory, white noise has zero autocorrelation at all lags, but in practice the estimated autocorrelations will vary slightly due to random sampling. Larger samples produce estimates that are closer to zero.

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

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

autoplot(amzn, Close) +
  labs(title = "Amazon Daily Closing Prices")

ggAcf(amzn$Close) +
  labs(title = "ACF of Amazon Closing Prices")

ggPacf(amzn$Close) +
  labs(title = "PACF of Amazon Closing Prices")

The plot of Amazon’s daily closing prices has a strong upward trend over most of the time, with changes in level over time. This means the series does not fluctuate around a constant mean, so it is non-stationarity. The ACF plot reinforces this. The autocorrelations are extremely high and positive for many lags, and they decrease really slowly. The PACF shows one very large spike at lag 1, followed by much smaller values after that. This tells us that the series has strong persistence and should be differenced before fitting an ARIMA model.

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

Turkish GDP from global_economy.

turkey_gdp <- global_economy |>
  filter(Country == "Turkey") |>
  select(Year, GDP)

autoplot(turkey_gdp, GDP) +
  labs(title = "Turkey GDP",
       x = "Year",
       y = "GDP")

The plot shows that the size of the fluctuations increases as the level of the GDP increases, so the variance increases. I will use the Guerrero method to find the optimal \(\lambda\) for the box-cox transformation.

lambda_turkey <- turkey_gdp |>
  features(GDP, guerrero) |>
  pull(lambda_guerrero)

lambda_turkey

## [1] 0.1572187

The \(\lambda\) is around 0.157, which is somewhat close to 0, so I will perform a log transformation.

turkey_gdp_t <- turkey_gdp |>
  mutate(GDP_trans = box_cox(GDP, lambda_turkey))

autoplot(turkey_gdp_t, GDP_trans) +
  labs(title = "Transformed Turkey GDP",
       x = "Year",
       y = "Transformed GDP")

turkey_gdp_t |>
  features(GDP_trans, unitroot_ndiffs)

## # A tibble: 1 × 1
##   ndiffs
##    <int>
## 1      1

After the transformation, the series needs one regular difference to become stationary. So the appropriate choice is a Box-Cox transformation with \(\lambda = 0.157\) and first differencing.

Accommodation takings in the state of Tasmania from aus_accommodation.

tas_accomodation <- aus_accommodation |>
  filter(State == "Tasmania") |>
  select(Date, Takings)

autoplot(tas_accomodation, Takings) +
  labs(title = "Accommodation in Tasmania",
       y = "Takings")

The series shows a clear upward trend, strong seasonality, and increasing variability.

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

lambda_tas

## [1] 0.001819643

The Box-Cox parameter (\(\lambda = 0.0018\)) suggests a log transformation.

tas_accomodation |>
  features(Takings, unitroot_ndiffs)

## # A tibble: 1 × 1
##   ndiffs
##    <int>
## 1      1

tas_accomodation |>
  features(Takings, unitroot_nsdiffs)

## # A tibble: 1 × 1
##   nsdiffs
##     <int>
## 1       1

After transformation, one regular difference and one seasonal difference are required to obtain stationarity.

Monthly sales from souvenirs.

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

The series shows a clear upward trend, strong seasonality, and increasing variability.

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

lambda_souv

## [1] 0.002118221

The Box-Cox parameter (\(\lambda = 0.0021\)) suggests a log transformation.

souvenirs |>
  features(Sales, unitroot_ndiffs)

## # A tibble: 1 × 1
##   ndiffs
##    <int>
## 1      1

souvenirs |>
  features(Sales, unitroot_nsdiffs)

## # A tibble: 1 × 1
##   nsdiffs
##     <int>
## 1       1

After transformation, one regular difference and one seasonal difference are required to obtain stationarity.

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

retaildata_raw <- "https://raw.githubusercontent.com/JDO-MSDS/Data-624/refs/heads/main/Homework%205/retail%20-%20Data1.csv"
retaildata <- read.csv(retaildata_raw, skip = 1)

myts <- ts(retaildata[,"A3349399C"],
           frequency = 12,
           start = c(1982,4))

autoplot(myts) +
  labs(title = "Monthly Australian Clothing Retail Sales",
       y = "Retail Turnover")

Trasformation needed?

lambda <- BoxCox.lambda(myts)
lambda

## [1] 0.02074707

Since lambda is 0.02074707, which is very close to zero, I will apply a log transformation.

myts_trans <- log(myts)
autoplot(myts_trans) +
  labs(title = "Log Transformed Clothing Retail Series",
       y = "log(Clothing Retail Turnover)")

The original retail series shows an upward trend and strong seasonal patterns. The seasonal extremes become larger as the level of the series increases, indicating that the variance grows over time. Because of this, I considered a Box-Cox transformation. The estimated parameter was 0.0207, which is very close to zero. This suggests that a log transformation is appropriate.

After applying the log transformation, the seasonal fluctuations seem a little more stable and the variance is more consistent over time. However, the transformed series still shows the same general trend and seasonality, but the scale is more suitable for further time series analysis.

ndiffs(myts_trans)

## [1] 1

nsdiffs(myts_trans)

## [1] 1

Since ndiffs = 1 and nsdiffs = 1, we need one seasonal difference (lag 12) and one regular difference.

myts_diff <- diff(myts_trans, lag = 12)
myts_diff2 <- diff(myts_diff)

autoplot(myts_diff2) +
   labs(title = "Differenced Log-Transformed Series",
       x = "Time",
       y = "Differenced log(Turnover)")

ACF of differenced series

ggAcf(myts_diff2) +
  labs(title = "ACF of Differenced Log-Transformed Series")

The transformed series still shows both trend and seasonality, so we need difference to make it stationary. Since the data is monthly but with a strong yearly seasonal pattern, one seasonal difference with lag 12 makes sense to remove the upward trend.

After applying both one seasonal and one regular difference, the series fluctuates around a roughly constant mean. Also, now there is no strong seasonal pattern anymore. So, the appropriate differencing order is one regular difference and one seasonal difference at lag 12.

Question 6: Simulate and plot some data from simple ARIMA models.

Use the following R code to generate data from an AR(1) model with \(ϕ_{1}=0.6\) and \(σ^2=1\). The process starts with \(y_1=0\).

set.seed(1234)

y <- ts(numeric(100))
e <- rnorm(100)

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

autoplot(y) +
  labs(
    title = "Simulated AR(1) Series",
    x = "Time",
    y = "y"
  )

The simulated series shows persistence over time, meaning that observations tend to be related to their recent past values. Because ϕ_{1}=0.6, the series tends to move smoothly and values stay above or below zero for short periods before changing direction.

Produce a time plot for the series. How does the plot change as you change \(ϕ_{1}\)?

simulate_ar1 <- function(phi, n = 100){
  y <- ts(numeric(n))
  e <- rnorm(n)
  
  for(i in 2:n)
    y[i] <- phi*y[i-1] + e[i]
  
  y
}

y1 <- simulate_ar1(0.2)
y2 <- simulate_ar1(0.5)
y3 <- simulate_ar1(0.8)

autoplot(cbind(phi_0.2 = y1,
               phi_0.5 = y2,
               phi_0.8 = y3),
         facets = TRUE) +
  labs(title = "AR1 Processes with Different Values of phi",
       x = "Time",
       y = "y")

By looking at the output plots for different ϕ, we can see that when ϕ₁ is small (0.2), the series seems like white noise, it looks random since it is constantly changing direction. As ϕ₁ increases, the series is smoother and with less changing in direction, meaning that values tend to remain above or below the mean for longer periods. For ϕ₁=0.8, the series has long swings (trends) and a strong dependence on past values.

Write your own code to generate data from an MA(1) model with \(θ_{1}=0.6\) and \(σ^2=1\).

MA(1) model: \(y_{t}= ε_{t} + θ_{1}ε_{t-1}\) where \(θ_{1}=0.6\) and \(σ^2=1\)

set.seed(1234)

y_ma1 <- ts(numeric(100))
e <- rnorm(100)

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

head(y_ma1, 10)

## Time Series:
## Start = 1 
## End = 10 
## Frequency = 1 
##  [1]  0.0000000 -0.4468102  1.2508987 -1.6950330 -0.9782939  0.7635307
##  [7] -0.2711064 -0.8914758 -0.8924311 -1.2287090

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

autoplot(y_ma1) +
  labs(title = "Simulated MA(1) Series (theta = 0.6)",
       x = "Time",
       y = "y")

simulate_ma1 <- function(theta, n = 100){
  y <- ts(numeric(n))
  e <- rnorm(n)
  
  for(i in 2:n){
    y[i] <- e[i] + theta * e[i-1]
  }
  
  y
}

set.seed(1234)

y1 <- simulate_ma1(0.2)
y2 <- simulate_ma1(0.5)
y3 <- simulate_ma1(0.8)

autoplot(cbind(theta_0.2 = y1,
               theta_0.5 = y2,
               theta_0.8 = y3),
         facets = TRUE) +
  labs(title = "MA(1) Processes with Different Values of theta",
       x = "Time",
       y = "y")

When θ₁ is small, the series behaves similarly to white noise with little a clear dependence between observations. As θ₁ increases, neighboring observations are more related to each other (short memory). However, unlike in the AR(1) process, the effect of extremes does not persist for long periods, that’s why the series still changes direction frequently.

Generate data from an ARMA(1,1) model with \(ϕ_{1}=0.6\), \(θ_{1}=0.6\) and \(σ_{2}=1\).

ARMA(1,1) model: \(y_{t} = ϕ_{1}y_{t-1} + ε_{t} + θ_{1}ε_{t-1}\)

set.seed(1234)

y_arma <- ts(numeric(100))
e <- rnorm(100)

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

# time plot
autoplot(y_arma) +
  labs(title = "Simulated ARMA(1,1) Series (phi = 0.6, theta = 0.6)",
       x = "Time",
       y = "y")

The ARMA(1,1) series combines characteristics of both AR and MA processes. The autoregressive component makes the series more persistence, while the MA component adds short-term effects from recent shocks/behavior. So, the series seems smoother than white noise but still has clear fluctuations over time.

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

AR(2) model: \(y_{t} = ϕ_{1}y_{t-1} + ϕ_{2}y_{t-2} + ε_{t}\) with \(ϕ_{1}=-0.8\) , \(ϕ_{2}=0.3\) and \(σ_{2}=1\)

set.seed(1234)

y_ar2 <- ts(numeric(100))
e <- rnorm(100)

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

# time plot
autoplot(y_ar2) +
  labs(title = "Simulated AR(2) Series (phi1 = -0.8, phi2 = 0.3)",
       x = "Time",
       y = "y")

The simulated AR(2) series has an oscillating pattern where the values repeatedly switch between positive and negative. Over time, the oscillations become larger (more amplitude), causing the series to diverge awayfrom zero, indicating that the process is non-stationary.

Graph the latter two series and compare them.

autoplot(cbind(ARMA11 = y_arma, AR2 = y_ar2), facets = TRUE) +
  labs(title = "Comparison of Simulated ARMA(1,1) and AR(2) Series",
       x = "Time",
       y = "y")

The two series have pretty difference behavior. The ARMA(1,1) series fluctuates around zero and remains relatively stable over time. The size of the fluctuations is pretty consistent, which is typical of a stationary time series.

In contrast, the AR(2) series shows oscillations that grow larger over time. The values constantly change between positive and negative, but the size of these swings increases with time. This causes the series to diverge away from zero, creating the growing behavior seen in the plot.

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

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.

fit_aus_airpassengers <- auto.arima(aus_airpassengers)
fit_aus_airpassengers

## Series: aus_airpassengers 
## ARIMA(0,2,1) 
## 
## Coefficients:
##           ma1
##       -0.8963
## s.e.   0.0594
## 
## sigma^2 = 4.308:  log likelihood = -97.02
## AIC=198.04   AICc=198.32   BIC=201.65

The auto.arima() function selected an ARIMA(0,2,1) model. Basically this means that the series needed two differences to achieve stationarity, which is makes sense considering the strong upward trend in the original data. The model has a moving average term (MA(1)). The estimated MA coefficient is −0.8963, so statistically significant, which means that there is a strong relationship between current and past shocks.

Checking Residuals

checkresiduals(fit_aus_airpassengers)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(0,2,1)
## Q* = 6.5774, df = 8, p-value = 0.5828
## 
## Model df: 1.   Total lags used: 9

We can confirm that the residuals look like white noise and that there is no significant autocorrelation remaining. So, the model properly captures the structure of the series.

forecast_fit_aus <- forecast(fit_aus_airpassengers, h = 10)

autoplot(forecast_fit_aus) +
  labs(title = "Forecasts of Australian Air Passengers",
       x = "Year",
       y = "Passengers [millions]")

Forecasts for the next 10 periods show an upward trend. Naturally, the prediction intervals are larger with time (more uncertainty).

In the end, the ARIMA(0,2,1) model provides a good fit to the data and generates reasonable forecasts that continue the previous growth pattern and still accounting for uncertainty.

Write the model in terms of the backshift operator.

For the fitted ARIMA(0,2,1) model, the backshift form is

\[(1-B)^2 y_t = (1+\theta_1 B)\varepsilon_t\] the \((1-B)^2 y_t\) is the second order differencing, the \((1+\theta_1 B)\) is the first order movinf average, and the \(\varepsilon_t\) is the white noise error term.

Since the estimated MA coefficient is \(\theta_1 = -0.8963\), the model becomes

\[(1-B)^2 y_t = (1-0.8963B)\varepsilon_t\]

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

ARIMA(0,1,0) with drift

y <- as.ts(aus_airpassengers)

# fit the 2 models
fit_021 <- Arima(y, order = c(0,2,1))
fit_010_drift <- Arima(y, order = c(0,1,0), include.drift = TRUE)

# forecast both
fc_021 <- forecast(fit_021, h = 10)
fc_010 <- forecast(fit_010_drift, h = 10)

# plot both
autoplot(y) +
  autolayer(fc_021$mean, series = "ARIMA(0,2,1)") +
  autolayer(fc_010$mean, series = "ARIMA(0,1,0) with drift") +
  labs(title = "Forecast Comparison",
       y = "Passengers [millions]") +
  guides(colour = guide_legend(title = "Model"))

Both models show upward trends, but ARIMA(0,2,1) increases slightly faster (higher slope), producing higher forecasts, while ARIMA(0,1,0) with drift follows a more steady linear growth.

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.

# with drift
fit_212_drift <- Arima(y, order = c(2,1,2), include.drift = TRUE)
fc_212_drift <- forecast(fit_212_drift, h = 10)

autoplot(fc_212_drift) +
  labs(title = "Forecasts from ARIMA(2,1,2) with Drift",
       y = "Passengers [millions]")

# without drift
fit_212_nodrift <- Arima(y, order = c(2,1,2), include.drift = FALSE, method = "ML")
fc_212_nodrift <- forecast(fit_212_nodrift, h = 10)

autoplot(fc_212_nodrift) +
  labs(title = "Forecasts from ARIMA(2,1,2) without Drift",
       y = "Passengers [millions]")

The ARIMA(2,1,2) model with drift shows steady growth similar to part (c), while the version without drift produces steeper, faster-growing forecasts. Removing the drift leads to more aggressive predictions. Removing drift actually mades forecasts steeper.

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

fit_021_const <- Arima(y, order = c(0,2,1), include.drift = TRUE)

## Warning in Arima(y, order = c(0, 2, 1), include.drift = TRUE): No drift term
## fitted as the order of difference is 2 or more.

fc_021_const <- forecast(fit_021_const, h = 10)

autoplot(fc_021_const) +
  labs(title = "ARIMA(0,2,1) with Drift",
       y = "Passengers [millions]")

When I tried to include a constant in the ARIMA(0,2,1) model, R returns a warning saying that no drift term is fitted because the order of differencing is 2. So, the constant is ignored and the model remains the same.

Question 8: For the usgdp series:

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

# plot to see behavior
autoplot(usgdp) +
  labs(title = "US GDP Time Series",
       x = "Year",
       y = "GDP")

# box cox parameter
lambda <- BoxCox.lambda(usgdp)
lambda

## [1] 0.366352

The plot of the usgdp series shows a strong upward trend and increasing variability over time, suggesting that a transformation may be useful to stabilize the variance.

The Box-Cox transformation parameter was estimated using BoxCox.lambda(), which gave lambda = 0.366. Since this value is between 0 and 1, a Box-Cox transformation is appropriate to reduce the increasing variance in the series.

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

# using auto arima with previous calculated lambda
fit_usgdp <- auto.arima(usgdp, lambda = lambda)
fit_usgdp

## Series: usgdp 
## ARIMA(2,1,0) with drift 
## Box Cox transformation: lambda= 0.366352 
## 
## Coefficients:
##          ar1     ar2   drift
##       0.2795  0.1208  0.1829
## s.e.  0.0647  0.0648  0.0202
## 
## sigma^2 = 0.03518:  log likelihood = 61.56
## AIC=-115.11   AICc=-114.94   BIC=-101.26

The auto.arima() function selected ARIMA(2,1,0) with drift and got the following coefficients: Estimated coefficients:

AR(1) = 0.2795
AR(2) = 0.1208
Drift = 0.1829

This means that the GDP series required one difference (d = 1) to remove the trend, and the remaining dynamics are captured by two autoregressive terms. The drift term represents the long-term upward growth in GDP.

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

Trying ARIMA(1,1,0), ARIMA(0,1,1), and ARIMA(1,1,1)

fit1 <- Arima(usgdp, order = c(1,1,0), lambda = lambda)
fit2 <- Arima(usgdp, order = c(0,1,1), lambda = lambda)
fit3 <- Arima(usgdp, order = c(1,1,1), lambda = lambda)

fit1

## Series: usgdp 
## ARIMA(1,1,0) 
## Box Cox transformation: lambda= 0.366352 
## 
## Coefficients:
##          ar1
##       0.6326
## s.e.  0.0504
## 
## sigma^2 = 0.04384:  log likelihood = 34.39
## AIC=-64.78   AICc=-64.73   BIC=-57.85

fit2

## Series: usgdp 
## ARIMA(0,1,1) 
## Box Cox transformation: lambda= 0.366352 
## 
## Coefficients:
##          ma1
##       0.4133
## s.e.  0.0453
## 
## sigma^2 = 0.05519:  log likelihood = 7.4
## AIC=-10.79   AICc=-10.74   BIC=-3.87

fit3

## Series: usgdp 
## ARIMA(1,1,1) 
## Box Cox transformation: lambda= 0.366352 
## 
## Coefficients:
##          ar1      ma1
##       0.9112  -0.5333
## s.e.  0.0461   0.1146
## 
## sigma^2 = 0.03979:  log likelihood = 46.21
## AIC=-86.43   AICc=-86.32   BIC=-76.04

I tested ARIMA(1,1,0), ARIMA(0,1,1), and ARIMA(1,1,1).

The resulting AIC values were:

ARIMA(1,1,0): AIC = −64.78
ARIMA(0,1,1): AIC = −10.79
ARIMA(1,1,1): AIC = −86.43

ARIMA(1,1,1) provided the best fit of the tested models, as it had the lowest AIC. However, the model selected earlier by auto.arima() (ARIMA(2,1,0) with drift - AIC = −115.11) still provides a better fit overall because it has a substantially lower AIC.

So, the ARIMA(2,1,0) with drift still is the preferred model for the transformed GDP series.

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

As explained above, the best model is ARIMA(2,1,0) since it has the lowest AIC.

Residuals for ARIMA(2,1,0)

checkresiduals(fit_usgdp)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(2,1,0) with drift
## Q* = 6.5772, df = 6, p-value = 0.3617
## 
## Model df: 2.   Total lags used: 8

The residual time plot shows that the residuals fluctuate randomly around zero with no obvious patterns. The histogram of residuals appears approximately normally distributed. Also, the residual ACF shows no significant spikes outside the confidence bounds, which means that there is little remaining autocorrelation.

The Ljung–box test gives a p-value of 0.3617 (>0.05). Sp, we fail to reject the null hypothesis that the residuals are independent.

The residual diagnostics tells us that the ARIMA(2,1,0) with drift model provides an adequate fit to the transformed GDP series.

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

usgdp_fc <- forecast(fit_usgdp, h = 10)
autoplot(usgdp_fc)

The forecasts show a continuation of the upward trend in US GDP. The prediction intervals are larger as the forecast horizon increases, showing a greater uncertainty in longer-term forecasts, as expected. The forecasts appear reasonable because they follow the long-term growth trend present in the historical GDP series and considering the greater uncertainty over time.

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

fit_ets <- ets(usgdp)
fit_ets

## ETS(A,A,N) 
## 
## Call:
## ets(y = usgdp)
## 
##   Smoothing parameters:
##     alpha = 0.9999 
##     beta  = 0.278 
## 
##   Initial states:
##     l = 1557.4589 
##     b = 18.6862 
## 
##   sigma:  41.8895
## 
##      AIC     AICc      BIC 
## 3072.303 3072.562 3089.643

usgdp_ets_fc <- forecast(fit_ets, h = 10)
autoplot(usgdp_ets_fc)

Using ets(usgdp), the selected model was ETS(A,A,N) (additive error, additive trend, no seasonality). The ETS forecasts also show a continued upward trend in US GDP, similar to the forecasts from the ARIMA(2,1,0) with drift model. Both models seem tocapture the strong growth pattern in the data, and the forecasts appear very similar. The main difference is that the ARIMA model used a Box-Cox transformation and differencing, while the ETS model directly models the level and trend components of the series.