How to calculate pinball loss for quantiles and for point forecasts?

Question

I have a few general questions about pinball loss:

1. Is a pinball loss typically calculated for each point in the forecast horizon or is it calculated across all points in the forecast horizon?

2. How is pinball loss typically calculated for a point forecast, which does not provide prediction quantiles?

3. Can someone provide a simple example or link to code for pinball loss. I have some sample code below, but it calculates loss at each point in the forecast horizon, maybe this is incorrect.

Formula for Pinball Loss

Pinball loss functional form I am using is.

$S_\tau(y, q) = 1_{y \geq q}\tau(y-q) + 1_{q > y}(1 - \tau)(q - y) = (y-q)(\tau - 1_{y - q < 0})$

Where $q$ is the quantile forecast, $y$ is the observed actual, $\tau$ is the quantile level ($\tau \in Q$ where $Q$ is the quantiles $[0.1, 0.2, \dots, 0.9]$), and $1$ is the indicator function. (Quantile regression for the qualifying match of GEFCom2017 probabilistic load forecasting)

pb_score <- function(y, q, tau){
  # y = actual value; q = quantile forecasted value; tau = quantile level
  indicator <- ifelse(y - q < 0, 1, 0)
  score <- (y - q) * (tau - indicator)
  return(score)
}

Sample of pinball loss calculation

# calculate the pinball loss for 1 series in the syph data set
# pinball loss is calculated at each point in the forecast
# horizon. 
require(tidyverse)   # for data manipulation
require(ZIM)         # for syph data
require(forecast)    # forcast package
require(janitor)     # clean names
data(syph)

# prep data
syph_ts <- ts(syph[,-1:-2], frequency = 52)
set.seed(42)

# ets
ets_fc <- function(ts_data, h){
  f = janitor::clean_names(as_tibble(forecast(train, h = h, level = seq(60,90,10))))
  horizon <- seq(1:h)
  cbind(horizon, f)
}

train <- subset(syph_ts[, 1],end=length(syph_ts[, 1])-9)
test <- subset(syph_ts[, 1],start=length(syph_ts[, 1])-8)

h <- length(test)
ets_fcast <- ets_fc(train, h = h)

# pinball loss
pb_score <- function(y, q, tau){
  # y = actual value; q = quantile forecasted value; tau = quantile level
  indicator <- ifelse(y - q < 0, 1, 0)
  score <- (y - q) * (tau - indicator)
  return(score)
}

y <- test

# ets pinball score
for(i in 1:h){
  ets_score <- rowMeans(rbind(
    pb_score(y = y[i], ets_fcast[i, ]$lo_90, tau = .1),
pb_score(y = y[i], ets_fcast[i, ]$lo_80, tau = .2),
    pb_score(y = y[i], ets_fcast[i, ]$lo_70, tau = .3),
pb_score(y = y[i], ets_fcast[i, ]$lo_60, tau = .4),
    pb_score(y = y[i], ets_fcast[i, ]$point_forecast, tau = .5),
pb_score(y = y[i], ets_fcast[i, ]$hi_60, tau = .6),
    pb_score(y = y[i], ets_fcast[i, ]$hi_70, tau = .7),
pb_score(y = y[i], ets_fcast[i, ]$hi_80, tau = .8),
    pb_score(y = y[i], ets_fcast[i, ]$hi_90, tau = .9)
  ))
}

Answer 1

1. Typically, loss functions for forecasts are computed for each horizon separately. This is partly because performance usually degrades the further out you go, and any improvements in the short term loss might be difficult to pick up when the long term loss is numerically much larger. Each horizon is effectively a different target and it makes sense to consider the model's performance on each one separately. There's nothing stopping you from averaging loss across horizons if that is meaningful for your problem domain, though.

   However, either way, you should repeat your training/testing split on a recursive or rolling basis, obtaining multiple loss values from different origins for each horizon, and average those. The single train/test split you have gives you effectively just one observation of performance for each horizon. This is true for any loss measure (not just pinball loss).

2. You can interpret a "pure" point forecast (with no associated prediction interval) as predicting a degenerate distribution with a mass at that point, i.e. $P[X_{t+h} = x] = 1$, where $x$ is your point forecast. In that case, for any $0 < \alpha < 1$, the predicted $\alpha$-quantile is simply $x$. This will probably not be a very good predicted quantile, but it does allow you to compare a point forecast method with other more complicated methods. It would not be unusual to find that a point forecast method gets a bad pinball loss in the tails but an actually better loss in the middle of the distribution than a model that tries to accurately reproduce all quantiles.

3. What you have there looks fine, aside from the fact that you will get better loss estimates from doing this from multiple origins. You should try plotting the result across quantile levels to see which parts of the distribution each model captures best.

   For example, here's your syph_ts[,1] and the performance for a normal white noise model, a normal AR(1) model, and the degenerate point forecast from the AR(1) model, where each model was first trained on the first 100 points, and then rolled forward 1 period at a time, and then the loss is averaged across folds for each horizon separately. This is the result one step ahead (here QS stands for quantile score, which is just twice the pinball loss):

And 10 steps ahead:

We can see that the AR(1) model dominates the white noise process in the short term but not the long term, and that the degenerate point forecast is worse in the tails but not too bad as an estimate in the 50-85% range (this is because the point forecast is biased high...).

If you are not really interested in the performance of the forecast at individual quantile levels, you might want to look at the continuous ranked probablity score (CRPS), which is the integral of the quantile score across quantile levels: