RW + drift prediction interval differs from textbook

[david-cortes]

I am trying to get a prediction interval for a simple drift model. I observe some small differences between the numbers produced by fable and the numbers that I would get by following the procedure from the reference that the docs point to:
https://otexts.com/fpp3/prediction-intervals.html

Suppose I generate a random series like this:

set.seed(123)
s0 <- rgamma(1,15,1)
slope <- rgamma(1,2,1)
slen <- 20
line <- s0 + slope * seq(1, slen)
s <- line * runif(slen, min=.7, max=1.3)
print(s)

 [1] 11.45998 19.38955 27.64921 26.48457 28.24855 41.17304 34.62969 43.08775 44.09223
[10] 34.70767 60.59940 44.23600 40.24937 52.73816 79.06283 80.72730 76.69836 78.13231
[19] 97.72110 86.03193

According to the reference, I should be able to generate a forecast following a random walk with drift model as follows:

est_slope <- (s[slen] - s[1]) / (slen - 1)
seq_fc <- seq(1, 10)
forecasted <- s[slen] + est_slope * seq_fc
forecasted

 [1]  89.95677  93.88161  97.80645 101.73129 105.65613 109.58097 113.50581 117.43064
 [9] 121.35548 125.28032

So far, so good - I get the same result from fable:

tbl <- data.frame(
    series = s,
    day = as.Date("2020-01-01") + seq(0, slen-1)
)
tbl <- as_tsibble(tbl, index="day")
tbl %>%
    model(pred = RW(series ~ drift())) %>%
    forecast(h = 10) %$%
    .mean

 [1]  89.95677  93.88161  97.80645 101.73129 105.65613 109.58097 113.50581 117.43064
 [9] 121.35548 125.28032

---

Now, according to the book, I should be able to calculate an 80% prediction interval with the following formula:

fitted <- s[1] + est_slope * seq(0, slen-1)
res_std <- sqrt(sum((s - fitted)^2) / (slen - 1))
alpha <- qnorm(0.9)
sd_multiplier <- sqrt(seq_fc * (1 + seq_fc / slen))
upper_pi_textbook <- forecasted + alpha * res_std * sd_multiplier
upper_pi_textbook

 [1] 101.6944 110.8717 119.0827 126.8274 134.2930 141.5723 148.7186 155.7655 162.7355
[10] 169.6443

But this doesn't match with fable's:

upper_pi_fable <- tbl %>%
    model(pred = RW(series ~ drift())) %>%
    forecast(h = 10) %>%
    hilo(80) %>%
    mutate(upper = `80%`$upper) %$%
    upper
upper_pi_fable

 [1] 105.2335 116.0474 125.6243 134.6085 143.2401 151.6377 159.8690 167.9765 175.9887
[10] 183.9256

---

Full code:

library(magrittr)
library(dplyr)
library(fable)

#### random data
set.seed(123)
s0 <- rgamma(1,15,1)
slope <- rgamma(1,2,1)
slen <- 20
line <- s0 + slope * seq(1, slen)
s <- line * runif(slen, min=.7, max=1.3)

#### textbook's formula
est_slope <- (s[slen] - s[1]) / (slen - 1)
seq_fc <- seq(1, 10)
forecasted <- s[slen] + est_slope * seq_fc

fitted <- s[1] + est_slope * seq(0, slen-1)
res_std <- sqrt(sum((s - fitted)^2) / (slen - 1))
alpha <- qnorm(0.9)
sd_multiplier <- sqrt(seq_fc * (1 + seq_fc / slen))
upper_pi_textbook <- forecasted + alpha * res_std * sd_multiplier


### fable
tbl <- data.frame(
    series = s,
    day = as.Date("2020-01-01") + seq(0, slen-1)
)
tbl <- as_tsibble(tbl, index="day")

upper_pi_fable <- tbl %>%
    model(pred = RW(series ~ drift())) %>%
    forecast(h = 10) %>%
    hilo(80) %>%
    mutate(upper = `80%`$upper) %$%
    upper

[robjhyndman]

There are two issues here.

1. Your calculation of fitted values is incorrect. Fitted values are one-step forecasts, and for a random walk with drift model are given by $y_{t-1}+\hat{b}$ or c(NA, s[-slen] + est_slope) . Because there is a missing observation at the start, and one parameter to be estimated, the unbiased estimator of the residual variance is

\hat{\sigma)^2 = \frac{1}{T-2}\sum_{t=2}^T (y_t - \hat{y}_t)^2

or sum((s - fitted)^2, na.rm=TRUE) / (slen - 2). (For some reason the LaTeX won't render, so I've omitted the math mode).)

2. The current fable code actually approximates $\hat{\sigma}^2$ by the MSE in one place:
   $$\hat{\sigma)^2 = \frac{1}{T}\sum_{t=2}^T (y_t - \hat{y}_t)^2$$

Putting these together, the existing fable code computes the upper limit like this:

mse <- mean((s - fitted)^2, na.rm=TRUE)
bvar <- sum((s - fitted)^2, na.rm=TRUE) / (slen - 2) / (slen-1)
alpha <- qnorm(0.9)
forecasted + alpha * sqrt(mse * seq_fc + bvar*seq_fc^2)

We should probably fix the approximation, so it uses the unbiased estimate and not the MSE. Then the equation in the book would be correct, except that the $T$ should be $T-1$. I'll fix that too.

[robjhyndman]

fable now using unbiased sigma^2: 159f19e
The book has also been updated as described above.
Thanks for pointing out the inconsistency

[david-cortes]

Thanks, can verify that after the update the numbers from fable match with the textbook's formula.

[mitchelloharawild]

@robjhyndman, this change to the forecast variance now introduces a difference between fable and forecast intervals. Should {forecast} also be updated?

[robjhyndman]

Yes, we should update the forecast package to use the unbiased estimate of variance.