Interactive online version:

Quantile Loss (Quantile Score)#

The quantile loss (quantile score) is a consistent scoring function for an \(\alpha\)-quantile forecast, where \(\alpha \in (0, 1)\). The quantile loss function is also used in Machine Learning, in quantile regression, which focuses on estimating different quantiles of the conditional distributution of the response variable. In finance, quantiles are often referred to as value-at-risk.

The quantile loss function is defined as:

\[\begin{split} S(x, y) = \begin{cases} \alpha \cdot (y - x) & \text{if } y \geq x \\ (\alpha - 1) \cdot (y - x) & \text{if } y < x \end{cases}\end{split}\]

where \(S\) is the scoring function (here quantile loss), \(x\) and \(y\) are forecast and observation, respectively. Lower values of the quantile loss are better.

More information about the quantile loss function can be found in this article:

Gneiting, T. (2011). Making and evaluating point forecasts. Journal of the American Statistical Association, 106, 746-762. https://doi.org/10.1198/jasa.2011.r10138

Now let’s look at an example where we demonstrate how to use the quantile score:

[1]:

from scores.continuous import quantile_score
import numpy as np
import pandas as pd
import xarray as xr

In this example, we create observation data and two synthetic forecasts, all with dimensions of ('lat', 'lon', 'time'). One forecast system is better at targeting the median (FCST_50), while the other one is better at targeting the 90th percentile (FCST_90).

[2]:

lat = np.linspace(-90, 90, 10)
lon = np.linspace(-180, 180, 20)
times = pd.date_range('2023-11-19', periods=5)
forecast_90th = np.random.uniform(0.7, 1.0, size=(len(lat), len(lon), len(times)))
forecast_90th_da = xr.DataArray(
    forecast_90th,
    dims=('lat', 'lon', 'time'),
    coords={'lat': lat, 'lon': lon, 'time': times}
)
forecast_50th = np.random.uniform(0.3, 0.6, size=(len(lat), len(lon), len(times)))
forecast_50th_da = xr.DataArray(
    forecast_50th,
    dims=('lat', 'lon', 'time'),
    coords={'lat': lat, 'lon': lon, 'time': times}
)

obs = np.random.rand(len(lat), len(lon), len(times))
obs_da = xr.DataArray(
    obs,
    dims=('lat', 'lon', 'time'),
    coords={'lat': lat, 'lon': lon, 'time': times}
)

Here, we want to calculate the quantile score of these two forecast systems for the median (\(\alpha=0.5\)) and 90th percentile (\(\alpha=0.9\)). In the quantile_score function, users can either calculate the overall score (i.e., preserve_dims=reduce_dims=None), so it will return one value that is the overall mean generalised piecewise linear (GPL) score, or calculate the score over specific dimension(s) by using preserve_dims or reduce_dims arguments. First let’s calculate the overall score of these forecasts for their targeted quantile (you always need to ensure that verification is aligned with the forecast definition):

[3]:

overall_qs_med = quantile_score(fcst=forecast_50th_da, obs=obs_da, alpha=0.5)
overall_qs_90th = quantile_score(fcst=forecast_90th_da, obs=obs_da, alpha=0.9)
print(
    f"Overall GPL score for:"
    f"\n FCST_50 (targeted at 50th percentile):  {np.round(overall_qs_med.values, 3)}"
    f"\n FCST_90 (targeted at 90th percentile): {np.round(overall_qs_90th.values, 3)}"
)

Overall GPL score for:
 FCST_50 (targeted at 50th percentile):  0.127
 FCST_90 (targeted at 90th percentile): 0.049

Now, let’s use preserve_dims=['time'] in our calculation to get quantile score values for each of the time periods:

[4]:

overall_qs_med_time = quantile_score(fcst=forecast_50th_da, obs=obs_da, alpha=0.5, preserve_dims=['time'])
overall_qs_90th_time = quantile_score(fcst=forecast_90th_da, obs=obs_da, alpha=0.9, preserve_dims=['time'])
print(
    f"Quantile score values by preserving time dimension"
    f"\n FCST_50 (targeted at 50th percentile):  {np.round(overall_qs_med_time.values, 3)}"
    f"\n FCST_90 (targeted at 90th percentile): {np.round(overall_qs_90th_time.values, 3)}"
)

Quantile score values by preserving time dimension
 FCST_50 (targeted at 50th percentile):  [0.116 0.125 0.137 0.126 0.13 ]
 FCST_90 (targeted at 90th percentile): [0.051 0.046 0.053 0.048 0.048]

Users have the option to input an array of weights within quantile_score to compute a weighted average score. An example scenario could involve assigning weights to regions based on their population when computing the quantile score. The following is an example of using weights, where we give larger weights to a few locations (lat, lon).

[5]:

weights_data = np.ones_like(obs)
weights_data[:, -2:, -2:] = 2
weights = xr.DataArray(weights_data, dims=('lat', 'lon', 'time'))
weighted_time_loss_score_med = quantile_score(
    fcst=forecast_50th_da, obs=obs_da, alpha=0.5, preserve_dims=['time'], weights=weights
)
weighted_time_loss_score_90th = quantile_score(
    fcst=forecast_90th_da, obs=obs_da, alpha=0.9, preserve_dims=['time'], weights=weights
)
print(
    f"GPL scores by preserving 'time' dimension and using an array of weights are:"
    f"\n FCST_50 (targeted at 50th percentile): {np.round(weighted_time_loss_score_med.values, 3)},"
    f"\n FCST_90 (targeted at 90th percentile): {np.round(weighted_time_loss_score_90th.values, 3)}"
)

GPL scores by preserving 'time' dimension and using an array of weights are:
 FCST_50 (targeted at 50th percentile): [0.116 0.125 0.137 0.139 0.144],
 FCST_90 (targeted at 90th percentile): [0.051 0.046 0.053 0.052 0.052]