Exercise Week 9: Solutions

library(fpp3)

fpp3 9.11, Ex6

Simulate and plot some data from simple ARIMA models.

  1. 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.
ar1 <- function(phi, n = 100L) {
  y <- numeric(n)
  e <- rnorm(n)
  for (i in 2:n) {
    y[i] <- phi * y[i - 1] + e[i]
  }
  tsibble(idx = seq_len(n), y = y, index = idx)
}
  1. Produce a time plot for the series. How does the plot change as you change \phi_1?

Some examples of changing \phi_1

ar1(0.6)
# A tsibble: 100 x 2 [1]
     idx      y
   <int>  <dbl>
 1     1  0    
 2     2 -1.17 
 3     3  0.323
 4     4 -1.01 
 5     5  1.29 
 6     6  1.69 
 7     7  1.67 
 8     8  0.897
 9     9  2.95 
10    10  1.77 
# ℹ 90 more rows
ar1(0.6) |> autoplot(y) + labs(title = expression(paste(phi, "=0.6")))

ar1(0.95) |> autoplot(y) + labs(title = expression(paste(phi, "=0.95")))

ar1(0.05) |> autoplot(y) + labs(title = expression(paste(phi, "=0.05")))

ar1(-0.65) |> autoplot(y) + labs(title = expression(paste(phi, "=-0.65")))

  • With \phi_1 close to 1, the series wanders, with long excursions away from zero (it is close to a random walk).
  • With \phi_1 close to 0, the series looks like white noise.
  • With \phi_1 < 0, the series oscillates rapidly around zero, as consecutive observations tend to have opposite signs.
  1. Write your own code to generate data from an MA(1) model with \theta_{1} = 0.6 and \sigma^2 = 1.
ma1 <- function(theta, n = 100L) {
  y <- numeric(n)
  e <- rnorm(n)
  for (i in 2:n) {
    y[i] <- theta * e[i - 1] + e[i]
  }
  tsibble(idx = seq_len(n), y = y, index = idx)
}
  1. Produce a time plot for the series. How does the plot change as you change \theta_1?
ma1(0.6) |> autoplot(y) + labs(title = expression(paste(theta, "=0.6")))

ma1(0.95) |> autoplot(y) + labs(title = expression(paste(theta, "=0.95")))

ma1(0.05) |> autoplot(y) + labs(title = expression(paste(theta, "=0.05")))

ma1(-0.65) |> autoplot(y) + labs(title = expression(paste(theta, "=-0.65")))

The MA(1) series all look fairly similar to white noise. With \theta_1 > 0, the series is a little smoother (consecutive observations are positively correlated), and with \theta_1 < 0 it is a little more jagged. The variance of the series, (1 + \theta_1^2)\sigma^2, increases with |\theta_1|.

  1. Generate data from an ARMA(1,1) model with \phi_{1} = 0.6, \theta_{1} = 0.6 and \sigma^2 = 1.
arma11 <- function(phi, theta, n = 100) {
  y <- numeric(n)
  e <- rnorm(n)
  for (i in 2:n) {
    y[i] <- phi * y[i - 1] + theta * e[i - 1] + e[i]
  }
  tsibble(idx = seq_len(n), y = y, index = idx)
}
arma11(0.6, 0.6) |> autoplot(y)

  1. Generate data from an AR(2) model with \phi_{1} = -0.8, \phi_{2} = 0.3 and \sigma^2 = 1. (Note that these parameters will give a non-stationary series.)
ar2 <- function(phi1, phi2, n = 100) {
  y <- numeric(n)
  e <- rnorm(n)
  for (i in 3:n) {
    y[i] <- phi1 * y[i - 1] + phi2 * y[i - 2] + e[i]
  }
  tsibble(idx = seq_len(n), y = y, index = idx)
}
ar2(-0.8, 0.3) |> autoplot(y)

  1. Graph the latter two series and compare them.

See graphs above. The ARMA(1,1) series is stationary, and wanders around zero with positive autocorrelation. The non-stationarity of the AR(2) process has led to increasing oscillations.

fpp3 9.11, Ex7

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

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

fit <- aus_airpassengers |>
  model(arima = 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
fit |> gg_tsresiduals()

fit |> forecast(h = 10) |> autoplot(aus_airpassengers)

An ARIMA(0,2,1) model was selected. The residuals look like white noise: there are no significant autocorrelations. The histogram suggests the residuals may not be normally distributed, but that only affects the prediction intervals.

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

(1 - B)^2y_t = (1 + \theta B)\varepsilon_t where \varepsilon_t \sim \text{N}(0,\sigma^2), \theta = -0.90 and \sigma^2 = 4.31.

  1. Plot forecasts from an ARIMA(0,1,0) model with drift and compare these to part a.
aus_airpassengers |>
  model(arima = ARIMA(Passengers ~ 1 + pdq(0, 1, 0))) |>
  forecast(h = 10) |>
  autoplot(aus_airpassengers)

  • Both contain increasing trends, but the ARIMA(0,2,1) model has an implicit trend due to the double-differencing, while the ARIMA(0,1,0) with drift models the trend directly via the drift coefficient.
  • The ARIMA(0,2,1) forecasts have a steeper slope (about 2.2 million passengers per year, compared with 1.4 million for the model with drift), as they follow the more recent growth.
  • The intervals are narrower when there are fewer differences.
  1. Plot forecasts from an ARIMA(2,1,2) model with drift and compare these to parts a and c. Remove the constant and see what happens.
aus_airpassengers |>
  model(arima = ARIMA(Passengers ~ 1 + pdq(2, 1, 2))) |>
  forecast(h = 10) |>
  autoplot(aus_airpassengers)

aus_airpassengers |>
  model(arima = ARIMA(Passengers ~ 0 + pdq(2, 1, 2)))
Warning: 1 error encountered for arima
[1] non-stationary AR part from CSS
# A mable: 1 x 1
         arima
       <model>
1 <NULL model>
  • The point forecasts from the ARIMA(2,1,2) model with drift are very similar to those from the ARIMA(0,1,0) model with drift, apart from some small oscillations in the first few years. The prediction intervals are slightly wider.
  • Removing the constant makes the estimation fail, with a “non-stationary AR part from CSS” error. Without the drift, the AR terms have to account for the trend, which pushes the AR polynomial towards non-stationarity.
  1. Plot forecasts from an ARIMA(0,2,1) model with a constant. What happens?
aus_airpassengers |>
  model(arima = ARIMA(Passengers ~ 1 + pdq(0, 2, 1))) |>
  forecast(h = 10) |>
  autoplot(aus_airpassengers)
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.

The forecast trend is now quadratic, and there is a warning that this is generally a bad idea. A constant with two differences implies a quadratic trend in the forecasts, which is rarely sensible.

fpp3 9.11, Ex8

For the United States GDP series (from global_economy):

  1. if necessary, find a suitable Box-Cox transformation for the data;
us_economy <- global_economy |>
  filter(Code == "USA")
us_economy |>
  autoplot(GDP)

lambda <- us_economy |>
  features(GDP, features = guerrero) |>
  pull(lambda_guerrero)
lambda
[1] 0.2819443
us_economy |>
  autoplot(box_cox(GDP, lambda))

The Guerrero method gives \lambda = 0.28. The transformed series is closer to linear, so a Box-Cox transformation may be useful here.

  1. fit a suitable ARIMA model to the transformed data using ARIMA();
fit <- us_economy |>
  model(ARIMA(box_cox(GDP, lambda)))
report(fit)
Series: GDP 
Model: ARIMA(1,1,0) w/ drift 
Transformation: box_cox(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

An ARIMA(1,1,0) model with drift was selected.

  1. try some other plausible models by experimenting with the orders chosen;
fit <- us_economy |>
  model(
    arima010 = ARIMA(box_cox(GDP, lambda) ~ 1 + pdq(0, 1, 0)),
    arima011 = ARIMA(box_cox(GDP, lambda) ~ 1 + pdq(0, 1, 1)),
    arima012 = ARIMA(box_cox(GDP, lambda) ~ 1 + pdq(0, 1, 2)),
    arima013 = ARIMA(box_cox(GDP, lambda) ~ 1 + pdq(0, 1, 3)),
    arima110 = ARIMA(box_cox(GDP, lambda) ~ 1 + pdq(1, 1, 0)),
    arima111 = ARIMA(box_cox(GDP, lambda) ~ 1 + pdq(1, 1, 1)),
    arima112 = ARIMA(box_cox(GDP, lambda) ~ 1 + pdq(1, 1, 2)),
    arima113 = ARIMA(box_cox(GDP, lambda) ~ 1 + pdq(1, 1, 3)),
    arima210 = ARIMA(box_cox(GDP, lambda) ~ 1 + pdq(2, 1, 0)),
    arima211 = ARIMA(box_cox(GDP, lambda) ~ 1 + pdq(2, 1, 1)),
    arima212 = ARIMA(box_cox(GDP, lambda) ~ 1 + pdq(2, 1, 2)),
    arima213 = ARIMA(box_cox(GDP, lambda) ~ 1 + pdq(2, 1, 3)),
    arima310 = ARIMA(box_cox(GDP, lambda) ~ 1 + pdq(3, 1, 0)),
    arima311 = ARIMA(box_cox(GDP, lambda) ~ 1 + pdq(3, 1, 1)),
    arima312 = ARIMA(box_cox(GDP, lambda) ~ 1 + pdq(3, 1, 2)),
    arima313 = ARIMA(box_cox(GDP, lambda) ~ 1 + pdq(3, 1, 3))
  )
  1. choose what you think is the best model and check the residual diagnostics;
fit |>
  glance() |>
  arrange(AICc) |>
  select(.model, AICc)
# A tibble: 16 × 2
   .model    AICc
   <chr>    <dbl>
 1 arima110  657.
 2 arima011  659.
 3 arima111  659.
 4 arima210  659.
 5 arima012  660.
 6 arima112  661.
 7 arima211  661.
 8 arima310  662.
 9 arima013  662.
10 arima312  663.
11 arima311  664.
12 arima113  664.
13 arima212  664.
14 arima313  665.
15 arima213  666.
16 arima010  668.

The best according to the AICc values is the ARIMA(1,1,0) w/ drift model, the same model that was chosen by ARIMA().

best_fit <- us_economy |>
  model(ARIMA(box_cox(GDP, lambda) ~ 1 + pdq(1, 1, 0)))
best_fit |> report()
Series: GDP 
Model: ARIMA(1,1,0) w/ drift 
Transformation: box_cox(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
best_fit |> gg_tsresiduals()

augment(best_fit) |> features(.innov, ljung_box, dof = 1, lag = 10)
# A tibble: 1 × 4
  Country       .model                                         lb_stat lb_pvalue
  <fct>         <chr>                                            <dbl>     <dbl>
1 United States ARIMA(box_cox(GDP, lambda) ~ 1 + pdq(1, 1, 0))    3.81     0.923

The residuals pass the Ljung-Box test, but the histogram looks negatively skewed.

  1. produce forecasts of your fitted model. Do the forecasts look reasonable?
best_fit |>
  forecast(h = 10) |>
  autoplot(us_economy)

These look reasonable. Let’s compare a model with no transformation.

fit1 <- us_economy |> model(ARIMA(GDP))
fit1 |>
  forecast(h = 10) |>
  autoplot(us_economy)

Without a transformation, ARIMA() selects an ARIMA(0,2,2) model. Notice the effect of the transformation on the forecasts. Increase the forecast horizon to see what happens. Notice also the width of the prediction intervals.

us_economy |>
  model(
    ARIMA(GDP),
    ARIMA(box_cox(GDP, lambda))
  ) |>
  forecast(h = 20) |>
  autoplot(us_economy)

The forecasts from the transformed data grow faster than those from the untransformed data, and their prediction intervals are wider and asymmetric.

  1. compare the results with what you would obtain using ETS() (with no transformation).
us_economy |>
  model(ETS(GDP)) |>
  forecast(h = 10) |>
  autoplot(us_economy)

ETS() selects an ETS(M,A,N) model. Its point forecasts are similar to those from the ARIMA model for the untransformed data (26.1 and 25.6 trillion dollars in 2027), and lower than those from the ARIMA model for the transformed data (28.2 trillion dollars in 2027). However, the ETS forecast intervals are much wider.

fpp3 9.11, Ex9

Consider aus_arrivals, the quarterly number of international visitors to Australia from several countries for the period 1981 Q1 – 2012 Q3.

  1. Select one country and describe the time plot.
aus_arrivals |>
  filter(Origin == "Japan") |>
  autoplot(Arrivals)

  • There is an increasing trend to about 1996, and slowly decreasing thereafter.
  • The seasonal shape has changed considerably over time.
  1. Use differencing to obtain stationary data.
aus_arrivals |>
  filter(Origin == "Japan") |>
  gg_tsdisplay(
    Arrivals |> difference(lag = 4) |> difference(),
    plot_type = "partial"
  )

  1. What can you learn from the ACF graph of the differenced data?
  • The non-seasonal lags suggest an MA(1) component.
  • The seasonal lags suggest a seasonal MA(1) component.
  1. What can you learn from the PACF graph of the differenced data?
  • The non-seasonal lags might suggest an AR(1) component. But the lag 2 spike in the PACF is larger than the lag 2 spike in the ACF. So an MA(1) is probably better.
  • The seasonal lags show geometric decay which is consistent with a seasonal MA(1).
  1. What model do these graphs suggest?
  • The suggested model is an ARIMA(0,1,1)(0,1,1)_4.
  1. Does ARIMA() give the same model that you chose? If not, which model do you think is better?
fit <- aus_arrivals |>
  filter(Origin == "Japan") |>
  model(ARIMA(Arrivals))
report(fit)
Series: Arrivals 
Model: ARIMA(0,1,1)(1,1,1)[4] 

Coefficients:
          ma1     sar1     sma1
      -0.4228  -0.2305  -0.5267
s.e.   0.0944   0.1337   0.1246

sigma^2 estimated as 174801727:  log likelihood=-1330.66
AIC=2669.32   AICc=2669.66   BIC=2680.54
cf <- tidy(fit) |>
  select(term, estimate) |>
  tibble::deframe()

The resulting model has an additional seasonal AR(1) component compared to what I guessed. We can compare the two models based on the AICc statistic:

aus_arrivals |>
  filter(Origin == "Japan") |>
  model(
    guess = ARIMA(Arrivals ~ pdq(0, 1, 1) + PDQ(0, 1, 1)),
    auto = ARIMA(Arrivals)
  ) |>
  glance()
# A tibble: 2 × 9
  Origin .model     sigma2 log_lik   AIC  AICc   BIC ar_roots  ma_roots 
  <chr>  <chr>       <dbl>   <dbl> <dbl> <dbl> <dbl> <list>    <list>   
1 Japan  guess  177223035.  -1332. 2670. 2670. 2678. <cpl [0]> <cpl [5]>
2 Japan  auto   174801727.  -1331. 2669. 2670. 2681. <cpl [4]> <cpl [5]>

The automatic model is only slightly better than my guess based on the AICc statistic.

  1. Write the model in terms of the backshift operator, then without using the backshift operator.

(1 - B)(1 - B^4)(1 - \Phi B^4)y_t = (1 + \theta B)(1 + \Theta B^4) \varepsilon_t \left[ 1 - B - (1 + \Phi)B^4 + (1 + \Phi) B^5 + \Phi B^8 - \Phi B^9 \right]y_t = (1 + \theta B + \Theta B^4 + \theta\Theta B^5) \varepsilon_t

y_t - y_{t-1} - (1 + \Phi)y_{t-4} + (1 + \Phi) y_{t-5} + \Phi y_{t-8} - \Phi y_{t-9} = \varepsilon_t + \theta \varepsilon_{t-1} + \Theta \varepsilon_{t-4} + \theta\Theta \varepsilon_{t-5}.

y_t = y_{t-1} + (1 + \Phi)y_{t-4} - (1 + \Phi) y_{t-5} - \Phi y_{t-8} + \Phi y_{t-9} + \varepsilon_t + \theta \varepsilon_{t-1} + \Theta \varepsilon_{t-4} + \theta\Theta \varepsilon_{t-5}, where \theta = -0.42, \Phi = -0.23 and \Theta = -0.53.

fpp3 9.11 Ex 10

Choose a series from us_employment, the total employment in different industries in the United States.

  1. Produce an STL decomposition of the data and describe the trend and seasonality.
leisure <- us_employment |>
  filter(Title == "Leisure and Hospitality")
leisure |>
  autoplot(Employed)

  • The sudden change in the seasonal pattern is probably due to some change in the definition of who is counted in this group. So our STL decomposition will need to have a small seasonal window to handle that.
  • In addition, the variation changes a little as the level increases, so we will also use a square root transformation.
leisure |>
  model(STL(sqrt(Employed) ~ season(window = 7))) |>
  components() |>
  autoplot()

  • With such a long series, it is not surprising to see the seasonality change a lot over time.
  • The seasonal pattern changed in the 1990s to what it is now.
  • The period of change was rapid, and the seasonal component hasn’t fully captured the change, leading to some seasonality ending up in the remainder series.
  • The trend is increasing.
  1. Do the data need transforming? If so, find a suitable transformation.

Yes. A square root did ok – the remainder series is relatively homoscedastic. No transformation or log transformations led to the remainder series appearing to be heteroscedastic.

leisure |> features(Employed, guerrero)
# A tibble: 1 × 2
  Series_ID     lambda_guerrero
  <chr>                   <dbl>
1 CEU7000000001          -0.216
  • The automatically selected transformation (\lambda = -0.22) is a little stronger than logs.
  • My preference is for something weaker.
  • I think the automatic procedure is confusing the changing seasonality with the increasing variance.
  1. Are the data stationary? If not, find an appropriate differencing which yields stationary data.
leisure |>
  autoplot(sqrt(Employed) |> difference(lag = 12) |> difference())

The double differenced square-rooted data is close to stationary, although the variance is not constant (it is smaller in the 1970s and 1980s than before or after), and there are some large values around 1990 when the seasonal pattern changed.

  1. Identify a couple of ARIMA models that might be useful in describing the time series. Which of your models is the best according to their AICc values?
leisure |>
  gg_tsdisplay(
    sqrt(Employed) |> difference(lag = 12) |> difference(),
    plot_type = "partial"
  )

  • The ACF and PACF of the doubly differenced data suggest that an ARIMA(2,1,0)(0,1,1)_{12} would be a good start: the PACF has a significant spike at lag 2, and the ACF has a large negative spike at lag 12.
  • An alternative would be an ARIMA(0,1,2)(0,1,1)_{12}.
fit <- leisure |>
  model(
    arima210011 = ARIMA(sqrt(Employed) ~ pdq(2, 1, 0) + PDQ(0, 1, 1)),
    arima012011 = ARIMA(sqrt(Employed) ~ pdq(0, 1, 2) + PDQ(0, 1, 1))
  )
glance(fit)
# A tibble: 2 × 9
  Series_ID     .model      sigma2 log_lik   AIC  AICc   BIC ar_roots  ma_roots
  <chr>         <chr>        <dbl>   <dbl> <dbl> <dbl> <dbl> <list>    <list>  
1 CEU7000000001 arima210011 0.0380    207. -406. -406. -386. <cpl [2]> <cpl>   
2 CEU7000000001 arima012011 0.0381    206. -404. -404. -384. <cpl [0]> <cpl>

The ARIMA(2,1,0)(0,1,1)_{12} model is better.

  1. Estimate the parameters of your best model and do diagnostic testing on the residuals. Do the residuals resemble white noise? If not, try to find another ARIMA model which fits better.
fit |>
  select(arima210011) |>
  gg_tsresiduals()

The tails of the residual distribution are too long, and there is significant autocorrelation at lag 11, as well as some smaller significant spikes elsewhere.

fit <- leisure |>
  model(
    arima210011 = ARIMA(sqrt(Employed) ~ pdq(2, 1, 0) + PDQ(0, 1, 1)),
    arima012011 = ARIMA(sqrt(Employed) ~ pdq(0, 1, 2) + PDQ(0, 1, 1)),
    auto = ARIMA(sqrt(Employed))
  )
glance(fit)
# A tibble: 3 × 9
  Series_ID     .model      sigma2 log_lik   AIC  AICc   BIC ar_roots  ma_roots
  <chr>         <chr>        <dbl>   <dbl> <dbl> <dbl> <dbl> <list>    <list>  
1 CEU7000000001 arima210011 0.0380    207. -406. -406. -386. <cpl [2]> <cpl>   
2 CEU7000000001 arima012011 0.0381    206. -404. -404. -384. <cpl [0]> <cpl>   
3 CEU7000000001 auto        0.0365    226. -440. -440. -411. <cpl [2]> <cpl>   
fit |> select(auto) |> report()
Series: Employed 
Model: ARIMA(2,1,2)(0,1,1)[12] 
Transformation: sqrt(Employed) 

Coefficients:
         ar1      ar2      ma1     ma2     sma1
      1.6261  -0.9132  -1.4773  0.7937  -0.5443
s.e.  0.0400   0.0309   0.0535  0.0352   0.0340

sigma^2 estimated as 0.03655:  log likelihood=226.22
AIC=-440.44   AICc=-440.35   BIC=-411.26

The automatically selected ARIMA(2,1,2)(0,1,1)_{12} model is better than either of my selections.

fit |>
  select(auto) |>
  gg_tsresiduals()

The residuals look better, although there is still a significant spike at lag 11.

  1. Forecast the next 3 years of data. Get the latest figures from https://fred.stlouisfed.org/categories/11 to check the accuracy of your forecasts.
fc <- fit |>
  forecast(h = "3 years")
fc |>
  filter(.model == "auto") |>
  autoplot(us_employment |> filter(year(Month) > 2000))

Data downloaded from https://fred.stlouisfed.org/series/CEU7000000001.

update <- readr::read_csv("CEU7000000001.csv") |>
  mutate(
    Month = yearmonth(observation_date),
    Employed = CEU7000000001
  ) |>
  select(Month, Employed) |>
  as_tsibble(index = Month) |>
  filter(Month >= min(fc$Month))
fc |> accuracy(update)
# A tibble: 3 × 10
  .model      .type     ME  RMSE   MAE   MPE  MAPE  MASE RMSSE  ACF1
  <chr>       <chr>  <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 arima012011 Test  -2887. 3426. 2887. -22.3  22.3   NaN   NaN 0.706
2 arima210011 Test  -2885. 3425. 2885. -22.3  22.3   NaN   NaN 0.706
3 auto        Test  -2931. 3463. 2931. -22.7  22.7   NaN   NaN 0.705
fc |>
  filter(.model == "auto") |>
  autoplot(us_employment |> filter(year(Month) > 2000)) +
  geom_line(data = update, aes(x = Month, y = Employed), col = 'red')

The forecasts from the automatically selected model are shown in blue, and the actual values in red. The initial forecasts look great, but then the pandemic led to a huge impact on the employment in this industry. All three models have similar (and large) errors over the three years, with the automatically selected model slightly worse than the other two.

  1. Eventually, the prediction intervals are so wide that the forecasts are not particularly useful. How many years of forecasts do you think are sufficiently accurate to be usable?
# Width of 95% prediction intervals relative to the point forecasts
fit |>
  select(auto) |>
  forecast(h = "10 years") |>
  mutate(h = row_number()) |>
  filter(h %in% c(1, 6, 12, 24, 36, 60, 120)) |>
  hilo(level = 95) |>
  unpack_hilo(`95%`) |>
  mutate(rel_width = (`95%_upper` - `95%_lower`) / .mean) |>
  as_tibble() |>
  select(Month, h, .mean, rel_width)
# A tibble: 7 × 4
     Month     h  .mean rel_width
     <mth> <int>  <dbl>     <dbl>
1 2019 Oct     1 16722.    0.0116
2 2020 Mar     6 16592.    0.0350
3 2020 Sep    12 17241.    0.0443
4 2021 Sep    24 17616.    0.0772
5 2022 Sep    36 17997.    0.112 
6 2024 Sep    60 18771.    0.190 
7 2029 Sep   120 20803.    0.418

The 95% prediction intervals are narrow: their width is about 4% of the point forecast after one year, 11% after three years, 19% after five years and 42% after ten years. On that basis, the forecasts would be usable for about five years, after which the intervals become too wide to be of much use. However, the intervals assume that the future is like the past, and they do not allow for shocks such as the pandemic. The actual values fell outside the 95% intervals from March 2020 onwards, so in practice only the first 5 months of forecasts were usable.