How do the forecast intervals from an AR model behave when the time series is inherently stationary?

Question

I'm trying to wrap my head around two contradictory intuitions behind how forecast intervals should behave when we use an AR process to model a stationary time series:

• (a) On one hand, since the time series is stationary, the variance is constant and therefore the forecast intervals should stay the same throughout the forecast horizon, since they are proportional to the standard deviation.

• (b) One the other hand forecast intervals should widen with each forecasting step, since intuitively, the uncertainty of our forecast should increase the further away it moves from the present.

To test this, I did the following:

n <- 150
eps <- rnorm(n)
x0 <- rep(0, n)
for (i in seq.int(2, n))
x0[i] <- 0.9 * x0[i-1] + eps[i]
plot(ts(x0))

library(forecast)
fit <- auto.arima(ts(x0))
fit
    Series: ts(x0) 
    ARIMA(1,0,0) with non-zero mean 

    Coefficients:
             ar1     mean
           0.8566  -1.0876
     s.e.  0.0412   0.5513

Ok, so far so good: The model fitted by auto.arima() is pretty close to the process I specified.

fcst <- forecast(fit, h=50)
autoplot(fust, ylab = 'AR(1) Process + intervals') 

The result I get is confusing. It supports neither (a) nor (b). The forecast intervals get larger and larger for the first few steps, then they stabilize and stay constant for the rest of the forecast horizon.

Hyndman and Athanasopoulos mention this in fpp:

In general, prediction intervals from ARIMA models increase as the forecast horizon increases. For stationary models (i.e., with $d = 0$) they will converge, so that prediction intervals for long horizons are all essentially the same.

But they don't elaborate on the reasons.

• What is the intuition behind the fact that the fact the intervals grow for a while then converge?
• And what is the math behind it?

Answer 1

You intuition is somewhat flawed here. Contrary to what you have asserted in (a), stationarity means that the marginal variance of the time-series values is constant, but this does not imply that the standard error of the prediction is constant. For future values that are close in time to your last observation, the standard error for the prediction interval starts off small (since there are less error terms between the two times) and then grows monotonically up to a fixed limit. To see why this is the case, I will show you some of the mathematics of prediction for the model you are using.

---

Prediction for stationary Gaussian $\text{AR}(1)$ model: For this model, each future value you are predicting can be written in terms of the observed sample values and unobservable errors as:

$$X_{n+k} = (1-\kappa^k) \mu + \kappa^k x_n + \sum_{i=0}^{k-1} \kappa^i \varepsilon_{n+k-i}.$$

where $\mu$ is the mean, $\kappa$ is the auto-regressive parameter, and $\sigma$ is the error standard deviation. Hence, treating the sample $\mathbf{x} = (x_1,...,x_n)$ as fixed values you have:

$$\mathbb{E}(X_{n+k}) = (1-\kappa^k) \mu + \kappa^k x_n \quad \quad \quad \quad \quad \mathbb{V}(X_{n+k}) = \sigma^2 \cdot \frac{1-\kappa^{2k}}{1-\kappa^2}.$$

To construct a prediction interval, you replace the parameters in the above equation with their estimators, and then use this to form a pivotal quantity for the prediction interval. This gives you the following predictions and corresponding standard errors:

$$\hat{X}_{n+k} = (1-\hat{\kappa}^k) \mu + \hat{\kappa}^k x_n \quad \quad \quad \quad \quad \hat{\text{se}}_{n+k} = \hat{\sigma} \cdot \sqrt{\frac{1-\hat{\kappa}^{2k}}{1-\hat{\kappa}^2}}.$$

As $k \rightarrow \infty$ you have:

$$\hat{X}_{n+k} \rightarrow \mu \quad \quad \quad \quad \quad \hat{\text{se}}_{n+k} \rightarrow \frac{\hat{\sigma}}{\sqrt{1-\hat{\kappa}^2}}.$$

In your analysis, the estimated auto-regressive parameter is near one, so the decay of powers of this parameter is quite slow. As you can see from the above equations, the prediction starts off close to the observed sample value, but converges to the true mean. Similarly, the standard error grows in size and converges to a constant value determined by the auto-regressive parameter and the error standard deviation. This is exactly what you see in your plotted prediction intervals.