How does the forecast package compute prediction intervals?

Question

I'm currently working with random walks with drift in R, I use the rwf formula from the forecast package and I wonder how the prediction intervals are computed. As I understand it, for the random walk with drift $$\gamma_t=\gamma_{t-1}+a+\epsilon_t$$ where $\epsilon\sim \mathcal{N}(0,\sigma^2)$,

the rwf function takes both uncertainty from the parameter estimate $\hat{a}$ and from the error term into account, but how exactly are they calculated?

I would say that the variance of the $m$th ahead forecast is equal to $$\mathbb{V}\gamma_{n+m}=m^2\frac{\hat{\sigma}^2}{n-1}+m\hat{\sigma}^2$$Then to compute the lower and uppper limits of the prediction band we would write $$\gamma_n+m\hat{a}\pm 1.96 \sqrt{m^2\frac{\hat{\sigma}^2}{n-1}+m\hat{\sigma}^2}$$ but this does not correspond to what the forecast package gives.

Answer 1

You can inspect the code of rwf by simply typing rwf.

In the case of a random walk with drift, an intercept-only OLS model is fitted to first differences:

$$\gamma_t-\gamma_{t-1}\sim a+\epsilon,$$

where $\epsilon\sim N(0,\sigma^2)$. Let $\hat{\sigma}^2$ denote the estimate of the residual variance $\sigma^2$ and $\hat{\sigma}_a$ the estimated standard error of the drift term $a$.

Per the code, the forecasted variance of $m$-step ahead future realizations is then estimated as

$$ \hat{s}^2 := (m\hat{\sigma}_a)^2+m\hat{\sigma}^2,$$

and this is then used with a standard $z$ score to yield prediction intervals:

$$ \gamma_n+m\hat{a}\pm z\sqrt{\hat{s}^2}.$$