library(fpp3)
library(patchwork)
library(purrr)

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

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

set.seed(42)
wn_36  <- rnorm(36)
wn_360 <- rnorm(360)
wn_1000 <- rnorm(1000)

wn_36_ts <- tsibble(idx = 1:36, y = wn_36, index = idx)
wn_360_ts <- tsibble(idx = 1:360, y = wn_360, index = idx)
wn_1000_ts <- tsibble(idx = 1:1000, y = wn_1000, index = idx)

p1 <- wn_36_ts |> ACF(y, lag_max = 20) |> autoplot() + labs(title = "n = 36")
p2 <- wn_360_ts |> ACF(y, lag_max = 20) |> autoplot() + labs(title = "n = 360")
p3 <- wn_1000_ts |> ACF(y, lag_max = 20) |> autoplot() + labs(title = "n = 1000")

p1/p2/p3

The main difference between figures is sample size. With 36 points ACF looks messy and has bigger spikes. With 360 and 1000 it becomes more stable and closer to zero.

Yes, all of them still white noise.

b. Why are the critical values at different distances from the mean of zero? Why are the autocorrelations different in each figure when they each refer to white noise?

Critical values are at different distance from zero because they depend on sample size. They are calculated as \(\frac{1.96}{\sqrt{n}}\) so when number of observations is small significance bounds are wider, and once the number of observation gets larger confidence bounds get closer to zero.

The number 1.96 comes from normal distribution (about 95% interval), since autocorrelation estimates are approximately normal around 0.

Autocorrelations are different because they are just sample estimates. Even if all data is white noise, each sample is random, so ACF values will not be exactly the same. With bigger sample they become closer to zero and more stable.

  1. 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.
#gafa_stock
#unique(gafa_stock$Symbol)
amazon <- gafa_stock |> 
  filter(Symbol == 'AMZN') |> 
  select(Close)

amazon |> 
  autoplot(Close)+
  labs(title = "Amazon closing stock prices")

amazon |> 
  gg_tsdisplay(Close, plot_type = 'partial', lag_max = 30)

The time plot of Amazon closing prices shows a clear uprising trend over time, so the mean is not constant over time. This strongly indicate that series is non-stationary. Stock time series also appear like a random walk model.

The ACF also supports non-stationary. The autocorrelations remain large and positive for many lags and decay very slowly, which is a sign of non-stationarity. For a stationary series, the ACF expected to drop toward zero much faster.

The PACF also shows a very large spike at lag 1, followed by smaller values afterwards and it means that current value depends on the previous one. This indicates that the series are not stationary.

So overall, all plots all suggest that the Amazon stock price series is non-stationary and should be differenced before fitting an ARIMA model.

amazon |> 
  features(Close, unitroot_ndiffs)

As we see the test confirms that differnce can be applied

amazon |> 
  gg_tsdisplay(difference(Close), lag_max = 60, plot_type = 'partial')

After differencing, the series looks much more stable around zero, except after 2018 where variability strongly increases.

ACF does not drop very cleanly and we still see spikes even up to higher lags (around 30–40), so there is still some remaining dependence.

PACF also shows some spikes, but no clear strong pattern.

So differencing improved stationarity, but series is not perfectly white noise. There is still some structure of the data left after differencing.

Let’s see if second-order differencing is needed

amazon |> 
  features(difference(Close), unitroot_ndiffs)

The test shows that there is no need for second-order differencing.

I think that for this stock price series a random walk model is a reasonable choice.

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

a. Turkish GDP from global_economy.

#global_economy

turkey_gdp <- global_economy |> 
  filter(Country == 'Turkey') |> 
  select(GDP, Population) |> 
  mutate(GDP = GDP/1e9)
turkey_gdp |> 
  autoplot(GDP)+
  labs(title = 'Turkish GDP', y = 'GDP (billions)')

The Turkish GDP series has strong upward trend and variance also increases over time, so transformation likely is needed.

lambda_turkey <- turkey_gdp |> 
  features(GDP, box_cox, features = guerrero) |> 
  pull(lambda_guerrero)
lambda_turkey
## [1] 0.1572187
turkey_gdp |> 
  autoplot(box_cox(GDP, lambda_turkey))

The Box-Cox transformation (lambda = 0.16) helps to stabilize the variance and make the trend close to linear.

turkey_gdp |> 
  features(box_cox(GDP, lambda_turkey), unitroot_ndiffs)
turkey_gdp |> 
  features(difference(box_cox(GDP, lambda_turkey) |> difference()), unitroot_ndiffs)

The test shows that only first-order differencing is appropriate.

turkey_gdp |> 
  gg_tsdisplay(difference(box_cox(GDP, lambda_turkey)), plot_type = 'partial')

After first differencing, the series looks much more stable around a constant mean.

From ACF and PACF plots i don’t see the significant spikes to manually identify ARIMA parameters.

So for this series, Box-Cox transformation and first-order differencing is appropriate.

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

#aus_accommodation
#unique(aus_accommodation$State)
tasmania_accomodations <- aus_accommodation |> 
  filter(State == 'Tasmania') |> 
  select(Takings)
tasmania_accomodations |> 
  autoplot(Takings)

As we see series show clear seasonal pattern and some increasing variation over time, so log transformation is appropriate to stabilize variance.

Since data is quarterly, seasonal differencing with lag 4 is applied to remove seasonality.

tasmania_accomodations |> 
  gg_tsdisplay(difference(log(Takings), 4) |> 
                 difference(1),
               plot_type = 'partial')

After seasonal differencing, the series looks stable around zero and the seasonal pattern is mostly gone.

The Kwiatkowski-Phillips-Schmidt-Shin (KPSS) test is used to determine if differencing is required. The hypotheses are:

\[ \begin{aligned} H_0 &: \text{The data are stationary} \\ H_a &: \text{The data are not stationary (unit root exists)} \end{aligned} \]

tasmania_accomodations |> 
  mutate(Takings_log = log(Takings)) |> 
  features(
    difference(Takings_log, 4),
    list(
      kpss = unitroot_kpss,
      ndiffs = unitroot_ndiffs
    )
  )

The KPSS test with p-value >= 0.1 fails to reject the null hypothesis of stationarity, so no additional differencing is needed.

So log transformation and seasonal differencing is enough to make the series stationary.

c. Monthly sales from souvenirs.

p1 <- souvenirs |> 
  autoplot(Sales)+
  labs(title = 'Souvenirs sales')
p1

The series has strong upward trend and strong seasonality, and variation increases over time, so log transformation likely is needed.

lambda_souvenirs <- souvenirs |> 
  features(Sales, box_cox, feature = guerrero) |> 
  pull(lambda_guerrero)
lambda_souvenirs
## [1] 0.002118221

Lambda is almost 0 so I will apply log transformation.

p2 <- souvenirs |> 
  autoplot(log(Sales))+
  labs(title = 'Log transformed souvenir sales')
p1/p2

Log transformation helped to improve the variance, but since the series has trend and seasonal patterns seasonal difference and maybe first-order differencing are needed.

souvenirs |> 
  ACF(lag_max = 36) |> 
  autoplot()
## Response variable not specified, automatically selected `var = Sales`

souvenirs |> 
  mutate(Sales_log = log(Sales)) |> 
  features(
    difference(Sales_log),
    list(
      kpss = unitroot_kpss,
      ndiffs = unitroot_nsdiffs
    )
  )

ACF plot shows significant spikes at multiples of 12, and KPSS test shows that seasonal-differencing is needed.

p3 <- souvenirs |> 
  gg_tsdisplay(difference(log(Sales), 12), plot_type = 'partial')
p3
## 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()`).

The mean is still not stable and variance is still changing over time, so I think series is not fully stationary.
Let’s test if first-order differencing is needed

souvenirs |> 
  mutate(Sales_log = log(Sales)) |> 
  features(
    difference(Sales_log, 12),
    list(
      kpss = unitroot_kpss,
      ndiffs = unitroot_ndiffs
    )
  )

Even though KPSS test does not reject stationarity, I still think plot after seasonal differencing shows that series is not fully stationary, so one more difference is needed.

p4 <- souvenirs |> 
  gg_tsdisplay(difference(log(Sales), 12) |> difference(), plot_type = 'partial')
p4
## 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()`).

After applying both seasonal and first differencing, the series looks stable around zero and variance is more constant.

ACF and PACF show negative significant spikes at lags1, but they are not consistent and do not follow a clear pattern.

So there is no strong structure left, and the series can be considered stationary.

  1. 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.
supermarket <- aus_retail |> 
  filter(Industry == 'Supermarket and grocery stores') |> 
  summarise(Turnover = sum(Turnover))

#supermarket
p1 <- supermarket |> 
  autoplot(Turnover) +
  labs(title = 'Supermarket and retail sales')
p1

The supermarket and retail sales series shows strong upward trend and clear seasonality, so it is not stationary.

There is also a presence of increasing variance along increasing level, indicating heteroskedasticity. Log or Box-Cox transformation is needed.

lambda_supermarket <- supermarket |> 
  features(Turnover, box_cox, feature = guerrero) |> 
  pull(lambda_guerrero)
lambda_supermarket
## [1] 0.1842339
p2 <- supermarket |> 
  autoplot(log(Turnover))+
  labs(title = "Log transformed supermarket and retail sales")

p3 <- supermarket |> 
  ACF(box_cox(Turnover, lambda_supermarket)) |> 
  autoplot()
p1/p2

I applied log transformation because Box-Cox lambda is close to zero. From the plot above we can see log transformation improved variance of the data.

The series clearly have a rising trend and seasonal patterns.

supermarket |> 
  mutate(Turnover_bxcx = box_cox(Turnover, lambda_supermarket)) |> 
  features(Turnover_bxcx, unitroot_nsdiffs)

Unit root test confirms that it is appropiate to apply seasonal differencing to the series

supermarket |> 
  gg_tsdisplay(difference(box_cox(Turnover, lambda_supermarket), 12), 
               lag = 36, plot_type = 'partial')

Time plot still appear to have a slight downward trend, and it appears that average value of the differenced series changing over time. Also it appears like wandering random walk.

ACF plot clearly shows that autocorrelation coefficient remain large and positive for many lags (as we see all the way up to lag 36 and probably beyond). It is obviously is not dropping to near zero coefficients quickly, as one of the requirements of stationary series.

From all findings from plots above I can say that series are not stationary after seasonal differencing, and likely first-order differencing is required to obtain stationary series.

Let’s do KPSS tests using the unitroot_ndiffs() feature to see if we need another difference

supermarket |> 
  features(difference(box_cox(Turnover, lambda_supermarket), 12), unitroot_ndiffs)

As we see the test shows that series are not stationary and we apply first-order differencing as well

supermarket |> 
 # mutate(Turnover_log_sdiff = difference(log(Turnover), 12)) |> 
  gg_tsdisplay(difference(box_cox(Turnover, lambda_supermarket), 12) |> 
                 difference(), lag = 36,
               plot_type = 'partial')

Time plot shows that series became stationary with first order differencing applied after seasonal differencing.
While the ACF and PACF plots show many significant spikes, this does not necessarily mean the series is non-stationary, it suggests that the series likely have strong seasonal or autoregressive components that need to be modeled using ARIMA terms (p, q, P, Q).

supermarket |> 
  mutate(Turnover_bxcx = box_cox(Turnover, lambda_supermarket)) |> 
  features(
    difference(difference(Turnover_bxcx, 12)),
    list(
      kpss = unitroot_kpss,
      ndiffs = unitroot_ndiffs
    )
  )

The tests confirm that series is stationary now, and no further differencing is required.

  1. Simulate and plot some data from simple ARIMA models.

a. Use the following R code to generate data from an AR(1) model with \(\phi_1=0.6\) and \(\sigma^2 = 1\) . The process starts with \(y_1 = 0\) .

set.seed(123)
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)
sim |> 
  autoplot()
## Plot variable not specified, automatically selected `.vars = y`

b. Produce a time plot for the series. How does the plot change as you change \(\phi_1\) ?

In the code below I simulated series generated from AR(1) model with \(\phi\) values in a range of (0.1, 1) and \(\sigma^2 = 1\)

set.seed(123)

phi_values <- c(0.1, 0.3, 0.6, 0.9, 1)
e <- rnorm(100)

sim_list <- lapply(phi_values, function(phi) {
  y <- numeric(100)

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

  data.frame(
    idx = 1:100,
    y = y,
    phi = paste0("phi = ", phi)
  )
})

sim_all <- do.call(rbind, sim_list) |>
  as_tsibble(key = phi, index = idx)

sim_all |>
  autoplot(y) +
  facet_wrap(~phi, ncol = 1, scales = "free_y")

An AR model works by making each value depend on its previous value plus some random noise, so the series carries information from the past forward instead of being completely random. As \(\phi\) increases, the influence of past shocks decays more slowly, which creates stronger dependence between observations and results in smoother time series.

c. Write your own code to generate data from an MA(1) model with \(\theta = 0.6\) and \(\sigma^2 = 1\) .

set.seed(123)

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)

sim |> 
  autoplot(y)

theta_values <- c(1, 0.8, 0.6, 0.25, 0)

e <- rnorm(100)

sim_list <- lapply(theta_values, function(theta) {
  y[1] <- e[1]
  
  for(i in 2:100){
    y[i] <- e[i] + theta * e[i-1]
  }
  
  data.frame(
    idx = 1:100,
    y = y,
    theta = paste0("theta = ", theta)
  )
})

sim_all <- do.call(rbind, sim_list) |>
  as_tsibble(key = theta, index = idx)

sim_all |> 
  autoplot(y) +
  facet_wrap(~theta,ncol =1, scales = "free_y")

First, I would like to say that we assume that the MA(1) model is used to simulate the data, and it is not a forecast. The reason is because current error \(e_t\) is unknown in forecast and usually is set to zero, while in simulation it is generated as random noise.

Including \(\theta\) = 0 provides a baseline of pure white random noise. As \(\theta\) increases, the series shows slightly stronger dependence between consecutive observations, but the effect remains small and the series stay noisy, because the MA(1) model does not include errors (random shocks) beyond one lag.

We can see that AR(1) model becomes smoother as value of \(\phi\) gets closer to 1, while in MA(1) model the series do not get much smoothness as \(\theta\) gets closer to 1. This is because AR model depends on previous observations, so when \(\phi\) = 1 basically it just adds to previous observation small random error.
On the other hand MA model is opposite - it depends on current random error (which could be completely different than previous one) and adds to it only a portion (depending on \(\theta\)) of previous random error.

e. Generate data from an ARMA(1,1) model with \(\phi_1\) = 0.6, \(\theta_1\) = 0.6 and \(\sigma^2\) = 1.

set.seed(123)

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

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

sim <- tsibble(idx = 1:100, y = y, index = idx)

p1 <- sim |> 
  autoplot(y) +
  labs(title = bquote("ARMA(1,1) with " ~ phi[1] == 0.6 ~ ", " ~ theta[1] == 0.6 ~ " and " ~ sigma^2 == 1))
p1

set.seed(123)

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

# initialize first two values
y[1] <- e[1]
y[2] <- -0.8*y[1] + e[2]

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

sim <- tsibble(idx = 1:100, y = y, index = idx)

p2 <- sim |> 
  autoplot(y) +
  labs(title = bquote("AR(2) with " ~ phi[1] == -0.8 ~ ", " ~ phi[2] == 0.3 ~ " and " ~ sigma^2 == 1))
p1/p2

According chapters 9.3 and 9.4 we normally restrict autoregressive and moving average models to stationary data, in which case some constraints on the values of the parameters are required:

-1 < \(\phi_1\) < 1,

\(\phi_1\) + \(\phi_2\) < 1

\(\phi_2\) - \(\phi_1\) <1

\(\theta_2\) + \(\theta_1\) > -1,

\(\theta_1\) - \(\theta_2\) < 1 .

In case of ARMA(1, 1) model both \(\phi_1\) and \(\theta_1\) satisfy stationary requeirments therefore series appear to be stationary.

On the other hand, AR(2) model has parameters \(\phi_1\) = -0.8 and \(\phi_2\) = 0.3. Both parameters satisfy the requirement of being constrained between -1 and 1, their sum is < 1.
But when it comes to their difference :

\(\phi_2\) - \(\phi_1\) = 0.3 +0.8 = 1.1

and it is not < 1, and therefore do not satisfy stationary constraint requirement \(\phi_2\) - \(\phi_1\) <1.
That’s why AR(2) graph oscillate with growing amplitude.

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

#aus_airpassengers
airpassengers <- aus_airpassengers |> 
  filter(Year <= 2011) 

airpassengers |> 
  autoplot(Passengers)+
  labs(title = "Air passengers 1970 -2011 (millions")

The plot of time series reveals following:

The data is not stationary because is has uprising trend and it is appropriate to use Box-Cox or Log transformation.

lambda_airpass <- aus_airpassengers |> 
  filter(Year<=2011) |> 
  features(Passengers, box_cox, feature = guerrero) |> 
  pull(lambda_guerrero)
lambda_airpass
## [1] 0.8927658

The Box-Cox test relsulted in lambda = 0.89, and it is close to 1. So Box-Cox transformation is not going too have big effect on the data. Therefore I will not apply Box-Cox transformation.

airpassengers |> 
  gg_tsdisplay(difference(Passengers), plot_type = 'partial')

After applying first-order differencing, the series appear to be close to be stationary and I’m not sure if second-order is needed.

fit_airpass <- airpassengers |> 
  model(
    auto_reg = ARIMA(Passengers))

fit_airpass |> 
  report()
## Series: Passengers 
## Model: ARIMA(0,2,1) 
## 
## Coefficients:
##           ma1
##       -0.8756
## s.e.   0.0722
## 
## sigma^2 estimated as 4.671:  log likelihood=-87.8
## AIC=179.61   AICc=179.93   BIC=182.99

Automatically selected model is ARIMA(0, 2, 1):

fit_airpass |> accuracy()
fit_airpass |> 
  select(auto_reg) |> 
  gg_tsresiduals()
## Warning: `gg_tsresiduals()` was deprecated in feasts 0.4.2.
## ℹ Please use `ggtime::gg_tsresiduals()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

The residuals look randomly scattered around zero up to around 1988, and there are fluctuations of residuals after that. Especially negative one around 1989 and very large one in 2010. These residuals are visible as ouliers on residual distribution plot.

The ACF of residuals does not show strong autocorrelation, most values are within bounds.

It looks like the model didn’t fully explains the structure of the series.

augment(fit_airpass) |> 
  features(.innov, ljung_box, h = 10, dof = 1)

The Ljung-Box test p-value = 0 suggests some remaining autocorrelation. This confirms the model may not fully capture all the structure in the data and maybe missing some patterns.

fc_airpass <- fit_airpass |> 
  forecast(h = 10)
p1 <- fc_airpass |> 
  autoplot(aus_airpassengers)+
  labs(title = 'ARIMA(0,2,1) 10 years forecast')
p1

As we can see from the forecast, the ARIMA(0, 2, 1) model overestimate future values after 2011 because it captures the strong increase observed around 2009–2011 as acceleration in the series. Because the model uses second differencing, it assumes that this accelerated growth rate will continue, which creates always growing forecast (that could be unrealistic in economic terms).

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

ARIMA(0, 2, 1) doesn’t have AR model part, includes first and second-order differencing with MA(1) whith only one parameter backshift formula is:

\((1-B)^2 y_t = (1 + \theta_1 B) \varepsilon_t\)

where \((1 - B)\) is squared because of double differencing.

And if we insert \(\theta_1\) = -0.8756 the final backshift formula:

\((1-B)^2 y_t = (1 - 0.8756 B) \varepsilon_t\)

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

fit_airpass <- airpassengers |> 
  model(
    auto_reg = ARIMA(Passengers),
    arima_010 = ARIMA(Passengers~ pdq(0,1,0)),
    arima_212 = ARIMA(Passengers~ pdq(2,1,2)),
    arima_212_noconst = ARIMA(Passengers ~ 0 + pdq(2,1,2)),
    arima_021_const = 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.
## Warning: 1 error encountered for arima_212_noconst
## [1] non-stationary AR part from CSS
fit_airpass |>
  select(model = 'arima_010') |> 
  report()
## Series: Passengers 
## Model: ARIMA(0,1,0) w/ drift 
## 
## Coefficients:
##       constant
##         1.3669
## s.e.    0.3319
## 
## sigma^2 estimated as 4.629:  log likelihood=-89.08
## AIC=182.17   AICc=182.48   BIC=185.59
fc_010 <-  fit_airpass |> 
  select(model = 'arima_010') |> 
  forecast(h=10)
p2 <- fc_010 |> 
  autoplot(aus_airpassengers)+
  labs(title = 'ARIMA(0,1,0) 10 years forecast')
  
p1/p2

Compared to ARIMA(0,2,1), the ARIMA(0,1,0) with drift gives a more linear and stable forecast as we see it in the plots.

We can see the last foerecast values for ARIMA(0,2,1) is almost 90, while ARIMA(0,1,0) last forecast is litlle bit above 75. It also confirms that ARIMA(0,1,0) is less aggressively increasing.
It assumes a constant average increase over time, instead of accelerating growth.
Because of all that, this forecast looks more realistic than ARIMA(0,2,1).

So ARIMA(0,1,0) with drift gives a more stable forecast for this series.

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.

fc_212 <-  fit_airpass |> 
  select(model = 'arima_212') |> 
  forecast(h=10)
p3 <- fc_212 |> 
  autoplot(aus_airpassengers)+
  labs(title = 'ARIMA(2,1,2) 10 years forecast')
  
p1/p3

The main difference is that ARIMA(2,1,2) has narrower prediction intervals, which means it produces more confident forecasts. This happens because ARIMA(2,1,2) uses more parameters (AR and MA terms), so it explains more of the structure in the data.

The forecast values themselves look very similar to ARIMA(0,2,1), so both models capture a similar trend in the data.

fc_212 <- fit_airpass |>
  select(arima_212) |>
  forecast(h = 10)

fc_212_noconst <- fit_airpass |>
  select(arima_212_noconst) |>
  forecast(h = 10)
p6 <- fc_212 |>
  autoplot(aus_airpassengers) +
  labs(title = "ARIMA(2,1,2) with drift")

p7 <- fc_212_noconst |>
  autoplot(aus_airpassengers) +
  labs(title = "ARIMA(2,1,2) without constant")

p1/p2/p6/p7
## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
## -Inf
## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
## -Inf
## Warning: Removed 10 rows containing missing values or values outside the scale range
## (`geom_line()`).

When trying to remove the constant, the model fails because AR part becomes non-stationary, so it cannot be estimated. This shows that not all parameter combinations are valid, and removing the constant might break the model.

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

fit_airpass |>
  select(model = 'arima_021_const') |> 
  report()
## Series: Passengers 
## Model: ARIMA(0,2,1) w/ poly 
## 
## Coefficients:
##           ma1  constant
##       -1.0000    0.0721
## s.e.   0.0709    0.0260
## 
## sigma^2 estimated as 4.086:  log likelihood=-85.74
## AIC=177.48   AICc=178.15   BIC=182.55
fc_021_const <-  fit_airpass |> 
  select(model = 'arima_021_const') |> 
  forecast(h=10)
p4 <- fc_021_const |> 
  autoplot(aus_airpassengers)+
  labs(title = 'ARIMA(0,2,1) with the constant = 0.072 10 years forecast')
  
p1/p4

When the constant is added to ARIMA(0,2,1), the forecast becomes more curved and grows faster over time.

This happens because with second differencing (d = 2), the constant introduces a quadratic trend, so the model assumes accelerating growth.

Compared to the model without constant, the forecast is more aggressive and increases more quickly.

This can lead to unrealistic forecasts, especially in long term.

  1. For the United States GDP series (from global_economy):
  1. if necessary, find a suitable Box-Cox transformation for the data
usa_gdp <- global_economy |> 
  filter(Code == 'USA') |> 
  select(GDP, Population)

p1 <- usa_gdp |> 
  autoplot(GDP/1e9)+
  labs(title = 'US GDP (billions)', y = 'GDP (billions)')
p1

The US GDP series shows strong upward trend and increasing variance over time, so it is not stationary.

Because variance increases with the level, a Box-Cox transformation is appropriate.

lambda_usgdp <- usa_gdp |> 
  features(GDP, box_cox, feature = guerrero) |> 
  pull(lambda_guerrero)
lambda_usgdp
## [1] 0.2819443

The estimated lambda is around 0.28, which means the transformation helps stabilize the variance.

usa_gdp <- usa_gdp |> 
  mutate(GDP_bxcx = box_cox(GDP, lambda_usgdp))

p2 <- usa_gdp |> 
  autoplot(GDP_bxcx)+
  labs(title = 'US GDP Box-Cox transformed. Lambda = 0.282')
p1/p2

A Box-Cox transformation is applied with lambda = 0.2819 to stabilize the variance.

However, from the plot we can see that the transformation helped to make the series linear with the strong trend.

fit_usgdp <- usa_gdp |> 
  model(
    auto_bxcx = ARIMA(box_cox(GDP, lambda_usgdp),
                      stepwise = FALSE,
                      approx = FALSE)
  )
fit_usgdp |> 
  report()
## Series: GDP 
## Model: ARIMA(1,1,0) w/ drift 
## Transformation: box_cox(GDP, lambda_usgdp) 
## 
## 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

The automatic ARIMA selected on the Box-Cox transformed series is ARIMA(1,1,0) with drift.

usa_gdp |> 
  features(difference(GDP_bxcx), unitroot_ndiffs)
usa_gdp |> 
  gg_tsdisplay(difference(GDP_bxcx), plot_type = 'partial')

From the plots, PACF shows a strong spike at lag 1, which suggests an AR(1) component, so ARIMA(1,1,0) is reasonable.

On the other hand, ACF also shows a significant spike at lag 1, which suggests an MA(1) component, so ARIMA(0,1,1) is also reasonable. At the same time ACF plot is sinusoidal, and that’s kinda puts doubt in this model’s selection.

So both models are only my estimated candidates, and more flexible models like ARIMA(1,1,1) or ARIMA(1,2,1) may fit better.

Let’s see how different ARIMA models will behave with training set. For training set

fit_usgdp <- usa_gdp |> 
  model(
    auto_bxcx = ARIMA(box_cox(GDP, lambda_usgdp),
                      stepwise = FALSE,
                      approx = FALSE),
    arima_011 = ARIMA(box_cox(GDP, lambda_usgdp)~ pdq(0,1,1)),
    arima_010 = ARIMA(box_cox(GDP, lambda_usgdp)~ pdq(0,1,0)),
    arima_111 = ARIMA(box_cox(GDP, lambda_usgdp)~ pdq(1,1,1)),
    arima_021 = ARIMA(box_cox(GDP, lambda_usgdp)~ pdq(0,2,1)),
    arima_121 = ARIMA(box_cox(GDP, lambda_usgdp)~ pdq(1,2,1)),
    arima_120 = ARIMA(box_cox(GDP, lambda_usgdp)~ pdq(1,2,0))
  )
fit_usgdp |> 
  accuracy() |> 
  arrange(RMSE)

Let’s see how different ARIMA models will behave with training set. For training set I choose the series except last 5 years.

train <- usa_gdp |> 
  filter(Year <= max(Year) - 5)

fit_usgdp_train <- train |> 
  model(
    auto_bxcx = ARIMA(box_cox(GDP, lambda_usgdp),
                      stepwise = FALSE,
                      approx = FALSE),
    arima_011 = ARIMA(box_cox(GDP, lambda_usgdp)~ pdq(0,1,1)),
    arima_110 = ARIMA(box_cox(GDP, lambda_usgdp)~ pdq(1,1,0)),
    arima_010 = ARIMA(box_cox(GDP, lambda_usgdp)~ pdq(0,1,0)),
    arima_111 = ARIMA(box_cox(GDP, lambda_usgdp)~ pdq(1,1,1)),
    arima_021 = ARIMA(box_cox(GDP, lambda_usgdp)~ pdq(0,2,1)),
    arima_121 = ARIMA(box_cox(GDP, lambda_usgdp)~ pdq(1,2,1)),
    arima_120 = ARIMA(box_cox(GDP, lambda_usgdp)~ pdq(1,2,0))
  )
fit_usgdp_train |> 
  accuracy() |> 
  arrange(RMSE)

From the training results, ARIMA(0,1,1), ARIMA(1,1,1), ARIMA(1,1,0), and the automatic model ARIMA(1,1,0) with drift all have very similar RMSE values.

This includes the models ARIMA(0,1,1) and ARIMA(1,1,0) I manually selected, so there is no clear best model based on training data alone.

The differences are very small, which suggests that several models can explain the data almost equally well.

But training results are not enough to choose the final model, and we need to evaluate performance on test data.

fc_usgdp_test <- fit_usgdp_train |> 
  forecast(h = 5)

fc_usgdp_test |> 
  accuracy(usa_gdp) |> 
  arrange(RMSE)

From the test results, ARIMA(1,2,1) has the lowest RMSE, so it performs best for forecasting.

All other models have much higher errors, even though some of them performed much better and similarly on training data.

This shows that ARIMA(1,2,1) suggesting that double differencing and both AR(1) MA(1) structures are needed to capture the structure of the series.

So ARIMA(1,2,1) is chosen as the final model.

Let’s take a look at forecast plots of my final selected and ARIMA function selected models

p1 <- fc_usgdp_test |> 
  filter(.model == 'arima_121') |> 
  autoplot(usa_gdp) +
  labs(title = 'ARIMA(1, 2, 1) forecast for test data')

p2 <- fc_usgdp_test |> 
  filter(.model == 'auto_bxcx') |> 
  autoplot(usa_gdp)+
   labs(title = 'ARIMA(1, 1, 0) forecast for test data')
p2/p1

As we see auto-ARIMA(1,1,0) slightly overestimates forecasts, while ARIMA(1,2,1) forecast points are almost identical with actual test points.

fit_usgdp_train |>
  select(arima_121) |>
  gg_tsresiduals()+
  labs(title = 'Residuals plots of ARIMA(1,2,1)')

fit_usgdp_train |>
  select(arima_121) |>
  augment() |>
  features(.innov, ljung_box, lag = 10, dof = 2)

Resudual plots of ARIMA(1,2,1) show that residuals are mostly centerd around zero, but there are some large spikes, especially around 1980 and 2008–2010. These negative residuals are results of economic recessions during those years, and I think no ARIMA models can capture and predict these shocks. The ACF of residuals does not show strong autocorrelation, and most values are within the bounds.

The residual distribution looks roughly symmetric, but there are some outliers which are occured during recessions. The Ljung-Box test has a high p-value=0.87, so it suggests that the residuals behave like white noise.

So overall, the model captures most of the structure, but not all extreme shocks in the data.

fit_usgdp_train |>
  select(auto_bxcx) |>
  gg_tsresiduals()+
  labs(title = "Residuals plots of ARIMA(1,1,0) with the drift")

The residuals are centered around zero, but there are still some large spikes, especially around 1980 and 2008–2010.

The ACF mostly stays within bounds.

The residual distribution is roughly symmetric, but with some outliers (because of mentioned recession shocks).

Compared to ARIMA(1,2,1), this model leaves slightly more structure in the residuals.

So in my opinion ARIMA(1,1,0) with drift is reasonable, but not as good as ARIMA(1,2,1).

fit_usgdp_train |>
  select(auto_bxcx) |>
  augment() |>
  features(.innov, ljung_box, lag = 10, dof = 1)

The Ljung-Box test has a high p-value=0.97, so it suggests that the residuals behave like white noise and the model is also adequate.

However, compared to ARIMA(1,2,1), this model performs worse on the test data in RMSE terms, so it is not the best choice for forecasting.

fit_arima_ets <- train |> 
  model(
    arima_121_bxcx = ARIMA(box_cox(GDP, lambda_usgdp) ~ pdq(1,2,1)),
    ETS_auto = ETS(GDP)
  )
fit_arima_ets |> 
  select(model = 'ETS_auto') |> 
  report()
## Series: GDP 
## Model: ETS(M,A,N) 
##   Smoothing parameters:
##     alpha = 0.9998967 
##     beta  = 0.4982127 
## 
##   Initial states:
##          l[0]        b[0]
##  448093333334 64619909089
## 
##   sigma^2:  7e-04
## 
##      AIC     AICc      BIC 
## 2901.713 2902.989 2911.564

ETS(M,A,N) model has \(\alpha\) = 0.99 (almost 1) and and is basically acting like a “memoryless” forecaster here. With an \(\alpha\) of almost 1, the model is just reacting to the last data observation and completely ignoring the past history. It is pretty much equivalent to a Naive forecast with a steep trend \(b\) attached.
In my opinion damped ETS(M, Ad, N) would be a better choice to forecast since since economy can’t grove forever.

fit_arima_ets |> 
  forecast(h = 5) |> 
  accuracy(usa_gdp) |> 
  select(.model, RMSE, MAE, MAPE)

The ETS model has much higher RMSE and MAE compared to ARIMA(1,2,1).

This means ETS performs significantly worse for forecasting this series.

ARIMA(1,2,1) has much lower error and gives more accurate predictions.

So ARIMA(1,2,1) is clearly the better model for this data.

fc_arima_ets <- fit_arima_ets |> 
  forecast(h=5)

fc_arima_ets |> 
  autoplot(usa_gdp, level = NULL)+
  autolayer(usa_gdp, GDP, color = 'black')+
  labs(title = 'ARIMA(1,2,1) versus ETS 5-year forecast comparison')+
  guides(color = guide_legend(title = 'Model'))

From the plot, both ARIMA(1,2,1) and ETS produce very similar forecast paths, following the same upward trend.

However, ARIMA(1,2,1) performs much better based on test accuracy, with significantly lower RMSE as we saw earlier.

But even ARIMA(1,2,1) performs better based on test accuracy, both models assume that GDP will continue growing over time.

In reality, GDP cannot grow forever at the same rate all the time, so these models may produce forecasts that become unrealistic in the long run. I would say that both models should be used with caution for long-term predictions.