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 😊

  • $\begingroup$ The basic problem with this solution is that it implicitly uses a model that is likely to differ (strongly) from the model you are using for forecasting. What would be the problem with applying the actual model you have adopted to vary the fit and the forecasts? That, at least, would be consistent. $\endgroup$
    – whuber
    May 31, 2022 at 22:16
  • $\begingroup$ Didn't your post at stats.stackexchange.com/questions/568828/…, which appears to have asked the same question (in a more constructive way), get an acceptable answer? $\endgroup$
    – whuber
    May 31, 2022 at 22:29
  • $\begingroup$ @whuber I am not sure I understand your comment. Would you please rephrase it? Also, I actually did ask the same question and got an acceptable answer but I could not use it as I code in python. I therefore wanted to develop the "simulate" function in python but needed feedback on it first. $\endgroup$
    – N-Light
    Jun 1, 2022 at 16:06

1 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
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))
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$$


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.