D624 Homework 6: ARIMA models

Exercise 9.1

Figure 9.32 (located here) shows the ACFs for 36 random numbers, 360 random numbers and 1,000 random numbers.

Part a

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

Since they are random numbers, they should be independent from each other, and so we would expect the Autocorrelation Function plots to be consistent with white noise, which they are except in one respect.

Each spike is the correlation between points with a given lag and the blue lines represent the 95% confidence interval that the correlation for a given lag is significant. Since this is random data we would expect one in 20 spikes to be outside the blue lines but none are. There may be issues with this random number generator even if, the data are behaving like white noise.

Notice as the amount of random numbers increases from 36 to 360 to 1,000 the blue 95% confidence intervals get closer to zero as with more data we have more evidence to support that there is no significant correlation between the data for different lags.

Part 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?

If the critical values, or spikes, were all exactly on the mean at zero I’m not sure if that’s perfect white noise or somehow randomness without the tell-tale traces of randomness. We expect the random numbers to have some correlation to each other at different lags and observe that as small fluctuations around zero. The fewer the amount of data points, the larger those spikes are going to be for the same reason that if you roll a die once it’s likely to get the highest and lowest possible numbers (1 and 6) but if you roll ten dice it’s incredibly unlikely to get the highest and lowest possible sums (10 and 60).

The autocorrelations are different to each other in each figure, even if the data in the smaller sets are reused in the larger sets, because there are more points to generate a critical value for at each lag.

Now I’m curious how brown noise would be characterised with an ACF plot.

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

The Partial Autocorrelation Function is the same as the ACF except it excludes the influence of the skipped over points of time when evaluating correlation at a given lag.

Since the PACF has a first value of basically 1 and then subsequent values at or within the confidence interval then each point is only related to the one immediately before it and then no others. After you account for the previous point you’ve already eliminated the influence of history on your current point and there’s no more influence in earlier plot. This can be seen in the ACF plot as well where we have the same starting value, basically one, and then since we’re not removing the intervening points when evaluating the lag we still show an autocorrelation across lags, however this is just a characteristic of random walks always starting from the point before it.

### Exercise 9.2

# Check out packages
library(fpp3)

# Select and show data, ACF and PACF plots
amazon <- gafa_stock |>
  filter(Symbol == "AMZN")
amazon |>
  gg_tsdisplay(Close, plot_type='partial')

Exercise 9.3

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

Part a

Turkish GDP from global_economy.

Below we find that one differencing is sufficient to obtain stationary data after the optimum box-cox transformation using a lambda of ~0.1572.

### Exercise 9.3 Part a

# Select Turkey's data
turkeygdp <- global_economy |>
  filter(Country == "Turkey")

# Find the right lambda for box-cox
lambda <- turkeygdp |>
  features(GDP, features = guerrero) |>
  pull(lambda_guerrero)

# Evidence stationarity with one difference
turkeygdp |>
  gg_tsdisplay(difference(box_cox(GDP, lambda)), plot_type='partial')

Part b

Accommodation takings in the state of Tasmania from aus_accommodation.

Below we find that one differencing is sufficient to obtain stationary data after the optimum box-cox transformation using a lambda of ~0.001819.

Since the lambda is essentially zero instead of doing a box-cox transformation we could just take the natural log of the Takings to ease communicating the methodology. (Alternative code using log instead of box-cox included as comment below)

### Exercise 9.3 Part b

# Select Tasmania's data
tasmaniatakings <- aus_accommodation |>
  filter(State == "Tasmania")

# Find the right lambda for box-cox
lambda <- tasmaniatakings |>
  features(Takings, features = guerrero) |>
  pull(lambda_guerrero)

# Evidence stationarity with one difference
tasmaniatakings |>
  gg_tsdisplay(difference(box_cox(Takings, lambda)), plot_type='partial')

# Alternative code using ln instead of box-cox since lambda is close to zero
#tasmaniatakings |>
#  gg_tsdisplay(difference(log(Takings)), plot_type='partial')

Part c

Monthly sales from souvenirs.

I thought I had to take two differences so that the data would be symmetric around zero but that’s not a condition for stationarity. Even though the mean is negative, it still has constant mean and variance.

To check I applied the Kwiatkowski-Phillips-Schmidt-Shin (KPSS) test to the souvenir sales data with zero, one, and two differences and got the following p-values:

Zero differences = 0.01 p-value
One differences = 0.10 p-value
Two differences = 0.10 p-value

KPSS test p-values get reported between 0.01 and 0.10 with larger and smaller values being reported as one of those two numbers. If the p-value is less than 0.05 we can reject the null-hypothesis that the data is stationary and so should try taking a difference to recheck stationarity.

As we can see from the p-values after only one differencing the data is stationary.

Note, similarly to part b, the lambda was close to zero at ~0.002118. This time we took the log of the data instead of doing box-cox.

### Exercise 9.3 Part c

# Check out packages
library(urca)

# Evidence stationarity with one difference
souvenirs |>
  gg_tsdisplay(difference(log(Sales), differences = 1), plot_type='partial')

# Apply the KPSS test to zero, one and two differences
kpss0 <- souvenirs |>
  features(Sales, unitroot_kpss)
kpss1 <- souvenirs |>
  features(difference(Sales), unitroot_kpss)
kpss2 <- souvenirs |>
  features(difference(Sales, differences = 2), unitroot_kpss)

Exercise 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

The minimum differencing for stationarity is one. We confirmed this with a KPSS test of 0.01 for zero differencing and 0.10 for one differencing.

Note, since the guerrero-lambda is close to zero at ~0.03543 we took the log of the data instead of box-cox.

### Exercise 9.5

# Select retail data
set.seed(20240915)
myseries <- aus_retail |>
  filter(`Series ID` == sample(aus_retail$`Series ID`,1))

# Lambda is close to zero so we're taking the log instead
lambda <- myseries |>
  features(Turnover, features = guerrero) |>
  pull(lambda_guerrero)

# Evidence stationarity with one difference
myseries |>
  gg_tsdisplay(difference(log(Turnover)), plot_type='partial')

# Confirming minimum differencing for stationarity is one
kpss0 <- myseries |>
  features(Turnover, unitroot_kpss)
kpss1 <- myseries |>
  features(difference(Turnover), unitroot_kpss)

Exercise 9.6

Simulate and plot some data from simple ARIMA models.

Part a

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.

-done, click Show button to the right to see code-

### Exercise 9.6 Part a

# AR(1) model with ϕ_1=0.6 and σ^2=1
y <- numeric(100)
set.seed(934289227)
e <- rnorm(100)
for(i in 2:100)
  y[i] <- 0.6*y[i-1] + e[i]
sim_0.6 <- tsibble(idx = seq_len(100), y = y, index = idx)

Part b

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

We can see that as phi approaches zero the data starts to look more like white noise and as phi approaches one the data starts to look like a random walk.

This is because as phi approaches zero we are reducing the influence of earlier states rapidly to the point of at zero there is no influence from earlier states and it’s truly independent states.

As phi approaches one we are increasing the influence of earlier states rapidly to the point of at one, we are now just adding the random error value to the preceding state.

### Exercise 9.6 Part b

# Check out packages
library(gridExtra)

# Generate three more AR(1) model with phis 0.05, 0.2 and 0.95
for(i in 2:100)
  y[i] <- 0.05*y[i-1] + e[i]
sim_0.05 <- tsibble(idx = seq_len(100), y = y, index = idx)
for(i in 2:100)
  y[i] <- 0.4*y[i-1] + e[i]
sim_0.4 <- tsibble(idx = seq_len(100), y = y, index = idx)
for(i in 2:100)
  y[i] <- .95*y[i-1] + e[i]
sim_0.95 <- tsibble(idx = seq_len(100), y = y, index = idx)

# Plot all four as a grid
phi_0.05 <- autoplot(sim_0.05) + 
  ggtitle("AR(1) model with ϕ_1=0.05 and σ^2=1") +
  xlab(NULL) +
  ylab(NULL)
phi_0.4 <- autoplot(sim_0.4) + 
  ggtitle("AR(1) model with ϕ_1=0.4 and σ^2=1") +
  xlab(NULL) +
  ylab(NULL)
phi_0.6 <- autoplot(sim_0.6) + 
  ggtitle("AR(1) model with ϕ_1=0.6 and σ^2=1") +
  xlab(NULL) +
  ylab(NULL)
phi_0.95 <- autoplot(sim_0.95) + 
  ggtitle("AR(1) model with ϕ_1=0.95 and σ^2=1") +
  xlab(NULL) +
  ylab(NULL)
grid.arrange(phi_0.05, phi_0.4, phi_0.6, phi_0.95, ncol=2)

Part c

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

-done, click Show button to the right to see code-

### Exercise 9.6 Part c

# MA(1) model with θ_1=0.6 and σ^2=1
y <- numeric(100)
set.seed(934289227)
e <- rnorm(100)
for(i in 2:100)
  y[i] <- e[i] + 0.6*e[i-1]
sim_ma_0.6 <- tsibble(idx = seq_len(100), y = y, index = idx)

Part d

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

Similarly to reducing phi in the AR(1) model close to zero, when you reduce theta in the MA(1) model close to 0, the data approaches white noise because previous errors are diminishingly included in the moving average.

And also, increasing theta towards 1 increases the auto correlation to previous errors, however, and I did not pick this up in my reading, it’s a negative autocorrelation, so instead of looking like a random walk, if the last value was positive, the next value will be negative unless the new error bumps it back over the line.

### Exercise 9.6 Part d

# Generate three more MA(1) model with thetas 0.05, 0.2 and 0.95
for(i in 2:100)
  y[i] <- e[i] + 0.05*e[i-1]
sim_ma_0.05 <- tsibble(idx = seq_len(100), y = y, index = idx)

for(i in 2:100)
  y[i] <- e[i] + 0.4*e[i-1]
sim_ma_0.4 <- tsibble(idx = seq_len(100), y = y, index = idx)

for(i in 2:100)
  y[i] <- e[i] + 0.95*e[i-1]
sim_ma_0.95 <- tsibble(idx = seq_len(100), y = y, index = idx)

# Plot all four as a grid
theta_0.05 <- autoplot(sim_ma_0.05) + 
  ggtitle("AR(1) model with θ_1=0.05 and σ^2=1") +
  xlab(NULL) +
  ylab(NULL)
theta_0.4 <- autoplot(sim_ma_0.4) + 
  ggtitle("AR(1) model with θ_1=0.4 and σ^2=1") +
  xlab(NULL) +
  ylab(NULL)
theta_0.6 <- autoplot(sim_ma_0.6) + 
  ggtitle("AR(1) model with θ_1=0.6 and σ^2=1") +
  xlab(NULL) +
  ylab(NULL)
theta_0.95 <- autoplot(sim_ma_0.95) + 
  ggtitle("AR(1) model with θ_1=0.95 and σ^2=1") +
  xlab(NULL) +
  ylab(NULL)
grid.arrange(theta_0.05, theta_0.4, theta_0.6, theta_0.95, ncol=2)

Part e

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

-done, click Show button to the right to see code-

### Exercise 9.6 Part e

# ARMA(1,1) model with ϕ_1=0.6, θ_1=0.6 and σ^2=1
y <- numeric(100)
set.seed(934289227)
e <- rnorm(100)
for(i in 2:100)
  y[i] <- 0.6*y[i-1] + 0.6*e[i-1] + e[i]
ARMA_1_1 <- tsibble(idx = seq_len(100), y = y, index = idx)

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

-done, click Show button to the right to see code-

### Exercise 9.6 Part f

# AR(2) model with ϕ_1=−0.8, ϕ_2=0.3 and σ^2=1
y <- numeric(100)
set.seed(934289227)
e <- rnorm(100)
for(i in 3:100)
  y[i] <- -0.8*y[i-1] + 0.3*y[i-2] + e[i]
AR_2 <- tsibble(idx = seq_len(100), y = y, index = idx)

Part g

Graph the latter two series and compare them.

Interesting. The ARMA(1,1) model looks like realistic data even though it was generated synthetically. It looks like the MA(1) term is softening the jagged edges of the AR(1) term.

The AR(2) model looks like a trumpet. The error terms normally distributed with a mean of zero and a standard deviation of one are not discernible because of the scale. The negative first phi causes the value to oscillate and the combination of the two phis causes the value to amplify with each state. I bet this is not invertible.

For an AR(2) model there are three constraints on the phis:
−1 < ϕ2 < 1
ϕ1 + ϕ2 < 1
ϕ2 − ϕ1 < 1

For ϕ1 = -0.8 and ϕ2 = 0.3, then the first and second constraints are met, but the third is not. That means if we were to increase -0.8 to greater than -0.7 or decrease 0.3 to less than 0.2 we should start to see normal behavior again.

### Exercise 9.6 Part g

# Plot parts e and f
plot_ARMA_1_1 <- autoplot(ARMA_1_1) + 
  ggtitle("ARMA(1,1) model with ϕ_1=0.6, θ_1=0.6 and σ^2=1") +
  xlab(NULL) +
  ylab(NULL)
plot_AR_2 <- autoplot(AR_2) + 
  ggtitle("AR(2) model with ϕ_1=−0.8, ϕ_2=0.3 and σ^2=1") +
  xlab(NULL) +
  ylab(NULL)
grid.arrange(plot_ARMA_1_1, plot_AR_2)

Exercise 9.7

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

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

ARIMA() selected an ARIMA(0,2,1) model.

The residuals look like white noise, though it looked like maybe the first half of the residuals were smaller than the second half, the critical values in the ACF plot are very close to 0.0 and the histogram looks a little too peaked and with too many outliers to be normally distributed. I tried running a box-cox transformation first but the three residual plots didn’t look different and the lambda was ~0.8927. Because of these two things I left the data untransformed for this exercise.

### Exercise 9.7 Part a

# Select passenger data
passengers <- aus_airpassengers |>
  filter_index(~ "2011") 

# Use ARIMA() to find a model
passengers_fit_a <- passengers |>
  model(ARIMA(Passengers))
passengers_fit_a

## # A mable: 1 x 1
##   `ARIMA(Passengers)`
##               <model>
## 1      <ARIMA(0,2,1)>

# Check residuals look like white noise
passengers_fit_a |>
  gg_tsresiduals()

# Plot forecasts for the next ten years
plota <- passengers_fit_a |>
  forecast(h=10) |>
  autoplot(passengers) +
  ggtitle("Forecast: ARIMA(0,2,1) model")
plota

Part b

Write the model in terms of the backshift operator.

Here is the model in backshift notation. We’re using the reported model coefficient below along with 2 for the differencing and 0 for the constant. Note the 1 below is the error term, not the constant.

\((1 - B)^2 y_t = (1 + -0.8863 B) \epsilon_t\)

### Exercise 9.7 Part b

# Obtain coefficients
report(passengers_fit_a)

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

Part c

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

The slope of the forecast is less steep than the one in part a. Also the prediction interval curves towards the center line while in part a they are straight lines like a lamp.

### Exercise 9.7 Part c

# Find an ARIMA(0,1,0) model with drift
passengers_fit_b <- passengers |>
  model(ARIMA(Passengers ~ pdq(0,1,0)))

# Plot forecasts for the next ten years
passengers_fit_b |>
  forecast(h=10) |>
  autoplot(passengers) +
  ggtitle("Forecast: ARIMA(0,1,0) model with drift")

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

The forecast center line and prediction interval boundaries now have a slight wave to them. The slope is higher than Part c’s. The prediction interval is the slimest of all three but start to curve out by year 10 so may have a larger prediction interval than the others by year 15 or 20.

When we remove the constant by adding “0 +” before the pdq() model doesn’t generate with the warning “non-stationary AR part from CSS”. (Note we could specify adding a constant with “1 +”.). My understanding is that by forcing the algorithm to look for a solution with four parameters, 0 + (2,1,2), it’s not able to find a model within the constraints.

### Exercise 9.7 Part d

# Find an ARIMA(2,1,2) model with drift
passengers_fit_c <- passengers |>
  model(ARIMA(Passengers ~ pdq(2,1,2)))

# Plot forecasts for the next ten years
passengers_fit_c |>
  forecast(h=10) |>
  autoplot(passengers) +
  ggtitle("Forecast: ARIMA(2,1,2) model with drift")

# Find an ARIMA(2,1,2) model with constant removed
passengers_fit_d <- passengers |>
  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

# Plot forecasts for the next ten years
passengers_fit_d |>
  forecast(h=10) |>
  autoplot(passengers) +
  ggtitle("Forecast: ARIMA(2,1,2) model with constant removed")

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

# Show model
report(passengers_fit_d)

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

Part e

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

The forecast starts to curve upwards and we get the following warning about having a “quadratic or higher order polynomial trend”.

### Exercise 9.7 Part e

# Find an ARIMA(0,2,1) model with a constant
passengers_fit_e <- passengers |>
  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.

# Plot forecasts for the next ten years
passengers_fit_e |>
  forecast(h=10) |>
  autoplot(passengers) +
  ggtitle("Forecast: ARIMA(0,2,a) model with a constant")

Exercise 9.8

For the United States GDP series (from global_economy):

Part a

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

The guerrero lambda is ~0.2819 so for ease in communicating methodology we’re taking out data to the quarter power which has a similar effect but is easier to explain.

### Exercise 9.8 Part a

# Select USA's data
usagdp <- global_economy |>
  filter(Country == "United States")

# Find the right lambda for box-cox
lambda <- usagdp |>
  features(GDP, features = guerrero) |>
  pull(lambda_guerrero)

# Show residuals with one difference
usagdp |>
  gg_tsdisplay(difference(GDP^0.25), plot_type='partial')

Part b

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

ARIMA() selected an ARIMA(1,1,0) model with drift or a constant.

\((1 - 0.4656 B) (1 - B) y_t = 11.5083 + \epsilon_t\)

### Exercise 9.8 Part b

# Use ARIMA() to find and report a model
usa_fit_1 <- usagdp |>
  model(ARIMA(GDP^0.25))
report(usa_fit_1)

## Series: GDP 
## Model: ARIMA(1,1,0) w/ drift 
## Transformation: GDP^0.25 
## 
## Coefficients:
##          ar1  constant
##       0.4656   11.5083
## s.e.  0.1191    0.9297
## 
## sigma^2 estimated as 52.46:  log likelihood=-192.84
## AIC=391.67   AICc=392.13   BIC=397.8

Part c

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

The minimum differencing for stationarity is one. We confirmed this with a KPSS test of 0.01 for zero differencing and 0.10 for one differencing.

Looking at the residuals in Exercise 9.8 Part a there is a significant spike at lag 1 in the PACF but none beyond lag 1 so I would choose to an ARIMA(1,1,0) model like what was selected, however we’ll try variations on that.

We tried six additional models and they all underperformed compared to the ARIMA() selected model (1,1,0), we also did an unbound search to find a higher order minimum but it returned the (1,1,0) model as well.

### Exercise 9.8 Part c

# Confirm differencing
kpss0 <- usagdp |>
  features(GDP^0.25, unitroot_kpss)
kpss1 <- usagdp |>
  features(difference(GDP^0.25), unitroot_kpss)

# Try several models
usa_fit_basket <- usagdp |>
  model(arima110 = ARIMA(GDP^0.25 ~ pdq(1,1,0)),
        arima011 = ARIMA(GDP^0.25 ~ pdq(0,1,1)),
        arima111 = ARIMA(GDP^0.25 ~ pdq(1,1,1)),
        arima211 = ARIMA(GDP^0.25 ~ pdq(2,1,1)),
        arima112 = ARIMA(GDP^0.25 ~ pdq(1,1,2)),
        arima212 = ARIMA(GDP^0.25 ~ pdq(2,1,2)),
        arima312 = ARIMA(GDP^0.25 ~ pdq(3,1,2)),
        stepwise = ARIMA(GDP^0.25),
        search = ARIMA(GDP^0.25, stepwise=FALSE))

# Show models
#usa_fit_basket |> pivot_longer(!Country, names_to = "Model name", values_to = "Orders")

# Display models ranked by AICc
glance(usa_fit_basket) |> arrange(AICc) |> select(.model:BIC)

## # A tibble: 9 × 6
##   .model   sigma2 log_lik   AIC  AICc   BIC
##   <chr>     <dbl>   <dbl> <dbl> <dbl> <dbl>
## 1 arima110   52.5   -193.  392.  392.  398.
## 2 stepwise   52.5   -193.  392.  392.  398.
## 3 search     52.5   -193.  392.  392.  398.
## 4 arima111   53.4   -193.  394.  394.  402.
## 5 arima011   54.8   -194.  394.  395.  400.
## 6 arima112   53.7   -192.  395.  396.  405.
## 7 arima211   53.8   -193.  395.  396.  405.
## 8 arima212   54.7   -192.  397.  399.  409.
## 9 arima312   56.5   -195.  401.  403.  413.

Part d

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

The best model is the ARIMA(1,1,0). Basically we did a manual step-wise search to confirm it.

Looking at the time series residual plots below, all autocorrelations are within the threshold limits which is consistent with white noise. (Note, the histogram looks skewed.)

Additionally we ran a ljung_box test and got a high p-value of ~0.7893, supporting that the residuals act like white noise and there’s no more autocorrelation to pull into the model. Note we used a lag of 10 for the ljung_box test because it’s not seasonal data, otherwise we would have doubled the number of periods in the season for the lag.

### Exercise 9.8 Part d

usa_fit_basket |>
  select(arima110) |>
  gg_tsresiduals()

# Ljung-box test
augment(usa_fit_basket) |>
  filter(.model=='arima110') |>
  features(.innov, ljung_box, lag = 10, dof = 3)

## # A tibble: 1 × 4
##   Country       .model   lb_stat lb_pvalue
##   <fct>         <chr>      <dbl>     <dbl>
## 1 United States arima110    3.92     0.789

Part e

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

The forecasts look reasonable.

### Exercise 9.8 Part e

# Forecast using best fit ARIMA model
arima_plot <- usa_fit_basket |>
  forecast(h=10) |>
  filter(.model=='arima110') |>
  autoplot(usagdp) +
  labs(title="10-year forecast: ARIMA model")
arima_plot

Part f

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

We have a flatter slope with the ETS model, so even though it’s ETS(M,A,N) without dampening it still looks more like a forecast than a hopecast.

The ARIMA prediction interval is narrower. While both the ETS and ARIMA model can’t predict changes in historical patterns, ARIMA models tend to have too narrow prediction intervals because they only account for variation in the errors and not variation in the parameter estimates and model order.

### Exercise 9.8 Part f

# Fit ETS model
ets_fit <- usagdp |>
  model(ETS(GDP))

# Plot 
ets_plot <- ets_fit |>
  forecast(h = 10) |>
  autoplot(usagdp) +
  labs(title="10-year forecast: ETS model")
ets_plot

References

These exercises come from ‘Forecasting: Principles and Practice’ by Rob J Hyndman and George Athanasopoulos, 3rd Ed
https://otexts.com/fpp3/

D624 Homework 6: ARIMA models

PK O’Flaherty

2024-10-13

Exercise 9.1

Part a

Part b

Exercise 9.2

Exercise 9.3

Part a

Part b

Part c

Exercise 9.5

Exercise 9.6

Part a

Part b

Part c

Part d

Part e

Part f

Part g

Exercise 9.7

Part a

Part b

Part c

Part d

Part e

Exercise 9.8

Part a

Part b

Part c

Part d

Part e

Part f

References