R Tutorial

Video tutorials are based on the old forecast package. Text tutorials have been updated to use fable and associated packages.

Video Tutorials

Evaluating Uncertainty Using Coverage

Evaluating How Forecast Accuracy Changes With Forecast Horizon

Written Tutorial

Setup

• Starting point from Evaluating Forecasts 1

library(tsibble)
library(fable)
library(feasts)
library(dplyr)

portal_data = read.csv("portal_timeseries.csv") |>
  mutate(month = yearmonth(date)) |>
  as_tsibble(index = month)

train <- filter(portal_data, month < yearmonth("2011 Dec"))
test <- filter(portal_data, month >= yearmonth("2011 Dec"))

ma2_model = model(train, ARIMA(NDVI ~ pdq(0,0,2) + PDQ(0,0,0)))
ma2_forecast = forecast(ma2_model, test)
models = model(train,
               ma2 = ARIMA(NDVI ~ pdq(0,0,2) + PDQ(0,0,0)),
               arima = ARIMA(NDVI))
forecasts = forecast(models, test)
autoplot(forecasts, train) + autolayer(test, NDVI)
accuracy(forecasts, test)

Introduction

• So far we’ve been dealing with point estimates
• Draw a time-series forecast with wide and narrow prediction intervals and a point outside the narrow, but inside the wide interval
• Point estimate comparisons say that these two forecasts are equivalent
• But the forecast with the wider interval is better
• This means we are just comparing the observed value to the mean predicted value
• Draw a set of axes with two parallel vertical lines
• Label 1 prediction and 1 observation

Coverage

• How do we think about uncertainty
• We have these blue prediction intervals, but how do we evalute them
• Prediction Interval: range of values in which a percentage of observations should occur
• 80% prediction interval is the range of values we expect 80% of observations to fall between

ma2_intervals <- hilo(ma2_forecast, level = 80) |>
  unpack_hilo("80%")
ma2_intervals$`80%_lower`
ma2_intervals$`80%_upper`

• We can find the observed points that occur in this range by checking for points that match both conditions
• So we want NDVI to be greater than lower and less than upper

in_interval <- test$NDVI > ma2_intervals$`80%_lower` &
  test$NDVI < ma2_intervals$`80%_upper`

• We can then determine what proportion of these values are TRUE, i.e., fall in the prediction interval

length(in_interval[in_interval == TRUE]) / length(in_interval)

• Called “coverage”
• Want coverage value to be as close to the value of the interval as possible
• So we want it to be close to 0.8

• Now it’s your turn.
• Write code to evaluate the coverage of the seasonal ARIMA model

• Let’s compare this result to the uncertainty of the seasonal model

arima_model = model(train, ARIMA(NDVI))
arima_forecast = forecast(arima_model, test)
arima_intervals <- hilo(arima_forecast, level = 80) |>
  unpack_hilo("80%")
in_interval = test$NDVI > arima_intervals$`80%_lower` & test$NDVI < arima_intervals$`80%_upper`
length(in_interval[in_interval == TRUE]) / length(in_interval)

• The full ARIMA is better because it is closer to the coverage interval of 0.8

Scores Incorporating Uncertainty

• Scores that incorporate uncertainty

• Reward prediction intervals that are just wide enough

• Instead of evaluating the mean prediction (point forecast) evaluate the prediction interval

• Add two sets of intervals on the axes

• Winkler Score

• Width of the prediction interval + a penalty for points outside the interval

• The width component rewards models with narrower prediction intervals

• The penalty rewards models without too many points outside the prediction intervals

• Penalties are calibrated to reward models with best coverage

• Include the Winkler score accuracy by adding two arguments

• A list, with a name for the score a function used to calculate it

• Any arguments that function requires (a prediction interval level in this case)

accuracy(forecasts, test, list(winkler = winkler_score), level = 80)

• Winkler requires choosing a single prediction interval

• Ideally we’d include information on lots of prediction intervals

• Instead of evaluting the mean check how likely an observation is given the full predicted distribution

• Add a distribution to the axes with the mode matching $\hat{y}$

• The best models most closely match the empirical distribution

• Doing this is technically complicated

• Continuous Rank Probability Score

• Scores each value relative to the predicted cumulative distribution function

• A value far from the mean is penalized less if the uncertainty is high

• We can add this by adding another model to our list

accuracy(forecasts, test, list(winkler = winkler_score, crps = CRPS), level = 80)

• And if we want to add back some of the other metrics we’d been seeing by default we can do that to

accuracy(forecasts, test, list(winkler = winkler_score, crps = CRPS, mae = MAE), level = 80)

Forecast horizon

• Forecasts generally get worse through time
• Can look at this by comparing the fit at different forecast horizons
• accuracy helps us do this with the by argument
• Evaluates the score separately for each unique value in a column

accuracies = accuracy(forecasts, test, list(winkler = winkler_score, crps = CRPS, mae = MAE), level = 80, by = c(".model", "month"))
accuracies

ggplot(data = accuracies, mapping = aes(x = month, y = crps, color = .model)) +
  geom_line()

• General trend towards increasing error with increasing horizon