estimation of a yearly prediction interval for monthly data

Question

I have a monthly time series and my objective is to provide my client with the next 12 point forecasts along with a yearly forecast. To obtain the yearly forecast, I simply summed up the 12 points forecasts. And, to associate this yearly forecast with a prediction interval, I came up with the following solution which does not make an assumption of the distribution of the residuals:

• Predict the 12 point forecasts using a model (ex. ARIMA or Prophet).
• Check that the residuals are uncorrelated. If it is the case then:
  • To each of the 12 point forecasts yhat_i, add an epsilon e_i drawn from the train residuals (yhat’_i = yhat_i + e_i) .
  • Sum the noisy yhat_i to obtain a noisy yearly forecast.
  • Repeat steps 2 and 3 a large number of times (ex. 1000).
  • Select the quantiles (alpha/2, 1-alpha/2) which constitute the yearly prediction interval with alpha being the risk level.

Let me know what you think 😊

Answer 1

That won't work for all models.

Here's an example where it won't work. Suppose the data comes from a random walk without drift:

$$Y_t = Y_{t-1} + \varepsilon_t$$

The $h$-ahead point forecasts are:

$$\mathbb{E}(Y_{t+h} | Y_1,...,Y_t) = Y_t$$

You want to simulate the yearly sum $\sum_{h=1}^{12} Y_{t+h}$ conditional on $Y_t$, which you can write as:

$$\sum_{h=1}^{12} Y_{t+h} = (Y_t + \varepsilon_{t+1}) + \ldots + (Y_t + \varepsilon_{t+1} + \ldots \varepsilon_{t+12}) = 12Y_t + \sum_{h=1}^{12} (13-h)\varepsilon_{t+h}$$

The residuals are just the $\varepsilon_t$, which in your suggestion you would resample. Let $\left(\tilde{\varepsilon}^{(1)},\ldots,\tilde{\varepsilon}^{(12)}\right)$ be such a sample. You would suggest to simulate from the yearly sum like this:

$$(Y_t + \tilde{\varepsilon}^{(1)}) + \ldots + (Y_t + \tilde{\varepsilon}^{(12)}) = 12Y_t + \sum_{h=1}^{12} \tilde{\varepsilon}^{(h)}$$

This is not the same distribution. Here's an illustration:

# Generate some data
set.seed(123)
x <- cumsum(rnorm(300))

residuals <- diff(x)

nsim <- 5000
h <- 12

# Sample of residuals
eps <- matrix(sample(residuals, nsim * h, replace = TRUE), nrow = nsim, ncol = h)

# Sample the correct paths, passing the shocks through the dynamics
sims_correct <- matrix(NA_real_, nsim, h)
for(i in 1:h) {
  if (i == 1) {
    sims_correct[,i] <- x[length(x)] + eps[,i]
  } else {
    sims_correct[,i] <- sims_correct[,i-1] + eps[,i]  
  }
}

# Sample paths incorrectly, adding the shock on top of the point forecast
sims_wrong <- matrix(NA_real_, nsim, h)
for (i in 1:h) {
  sims_wrong[,i] <- x[length(x)] + eps[,i]
}

# Add up the simulated months
yearly_sim_correct <- rowSums(sims_correct)
yearly_sim_wrong <- rowSums(sims_wrong)

# Compare density
plot(density(yearly_sim_correct), main = "Yearly predicted density", ylim =     c(0,0.15))
lines(density(yearly_sim_wrong),col='red')
legend("topleft", legend = c("Correct", "Wrong"), col = 1:2, lty=1)

The issue is that you're just adding noise "on top" of the point forecast, but in many time series models, shocks are defined in a way that they are fed forward through the dynamics. That's how autocorrelation is constructed from independent shocks. In the random walk process, each new shock has a permanent effect on the future level and appears in the distribution of each value from then on. You need to do the same with resampled residuals so that you obtain the correct joint pathwise distribution and hence the correct summed yearly distribution.

So, this won't work for most ARIMA models. There are models where your idea will work as-is, namely models that have a deterministic trend function $f$ and iid noise, like this:

$$Y_t = f(t) + \varepsilon_t$$