Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
2,375,935 characters · 27 sections · 47 citation commands
Multivariate Forecasting Evaluation: On Sensitive and Strictly Proper Scoring Rules
\lhead{\nouppercase{\leftmark}}
{\bf Keywords:} Scoring rules, Multivariate forecasting, Prediction evaluation, Ensemble forecasting, Energy score, Copula score, Strictly proper, Diebold-Mariano test
Forecasting evaluation is still highly discussed in different forecasting communities, like in meteorology, energy, economics, business or social and natural sciences. Also, in recent years, the term probabilistic forecasting (essentially density forecasting or quantile forecasting) became popular across all areas of application (gneiting2014probabilistic, reich2015probabilistic, hong2016probabilistic, lerch2017forecaster). It is clear that the evaluation of forecasts is crucial for assessing the quality of different forecasts. However, for 1-step ahead forecasts the discussion is much more comprehensive than for multiple step ahead predictions, say $H$-step ahead forecasts.
The evaluation of forecasted marginal distributions of $H$-step ahead forecasts have been treated in literature using so called univariate scoring rules. Essentially it breaks down to the 1-step ahead forecasting case for each marginal distribution. In contrast, the (fully) multivariate case where the complete $H$-dimensional distribution is usually not covered, in the forecasting itself such as in the evaluation. Obviously, this causes problems because the multivariate $H$-step ahead distribution is not fully explained by reporting $H$-step ahead marginal distributions or quantiles. A main reason for these circumstances is the lack of suitable evaluation measures that can judge the quality of forecasts adequately. We require approaches to assess not only the performance of a prediction with respect to the marginals but also to the dependency structures, we need suitable multivariate scoring rules.
The most well known multivariate scoring rule, the energy score, was introduced by gneiting2007strictly. For instance, in sari2016statistical, spath2015time, gneiting2008assessing, pinson2012adaptive, berrocal2008probabilistic, pinson2012evaluating, moller2013multivariate, baran2015joint, yang2015multi, moller2015spatially the energy score is used for evaluation. This is mainly in the area of meteorologic forecasting, primarily wind speed and wind power forecasting. In the wind power forecasting review of zhang2014review it is the only mentioned multivariate forecasting evaluation method.
The energy score is a strictly proper scoring, meaning that only the true model optimizes the corresponding score. However, in many empirical applications the energy score is still avoided. In the study of pinson2013discrimination, the discrimination ability of the energy score is checked in a simulation study while considering the bivariate normal distribution. They conclude that the energy score can not separate differences in the dependency structure well, which made the researchers and practitioners skeptical in using it. In this paper, we will discuss this study in detail and draw somewhat different conclusions.
In contrast, another multivariate scoring rule were introduced to overcome the reported problem in energy score. scheuerer2015variogram developed the variogram score, that is sensitive to changes in the correlation. However, it is only a proper scoring rule, but not a strictly proper one. Thus, it can not identify the true underlying model. Another reported plausible candidate is the multivariate log-score. However, it requires that we have multivariate density forecasts which is not the case in many applications. As discussed in lerch2017forecaster, this density may be approximated. But these approximation methods suffer efficiency in higher dimensions. The results depend crucially on the chosen approximation method. Additionally, we have the Dawid-Sebastiani score (see e.g. gneiting2007strictly), which evaluates the mean and covariance matrix of the distributions and corresponds to the log-score in the multivariate Gaussian settings. Moreover, a characterization only by the first two moments is not sufficient for many applications.
With all these considerations in mind, we introduce a new scoring rule that is sensitive to dependency changes in the distribution using copula theory. By Sklar's theorem, we are able to extract the dependency structure of our forecast (in the form of the copula) and construct a more sensitive measure for evaluating the dependencies. Afterwards, we combine the score of the copula with a score of the marginal distributions to obtain a proper and even strictly proper scoring rules. The approach is somehow flexible with respect to the copula and marginal scores that are used.
We start with the introduction of the aforementioned famous scoring rules in section (ref). Afterwards, we introduce our concept of marginal-copula scores. In section (ref) we give a short overview on how to report multivariate forecasts. Efficient estimation of scores is crucial in practice, we describe some approaches for different scoring rules in section (ref). Since the sensitivity of a score value alone does not indicate very much about the goodness-of-fit or the discrimination ability, we introduce the concept of the Diebold-Mariano test in section (ref). In section (ref), we apply the renowned scoring rules and the marginal-copula scores in three different synthetic case studies, compare and discuss the results. A real data example using airline passenger data is discussed in section (ref). We conclude the results in section (ref).
In this section we are going to introduce a concept of scoring rules to adequately assess the goodness-of-fit of multivariate probabilistic forecasts. In order to do that, we are going to introduce some well known scoring rules in a first step. Additionally, we present a new score which combines univariate scoring of the forecasted marginal distributions as well as multivariate scoring for the forecasted copula in a way that desirable preferences like (strict) propriety are preserved. We start with some notations and basics on scoring rules.
Let $\boldsymbol Y = (Y_1,\ldots, Y_H)'$ be an $H$-dimensional random variable with multivariate (cumulative) distribution function $\boldsymbol F_{\boldsymbol Y}$ which we want to forecast. The most standard example would be that $\boldsymbol Y$ is a $H$-step ahead forecast of a univariate time series, hence $H$ is the forecasting horizon. From $\boldsymbol Y$ we observe the realized respectively materialized vector $\boldsymbol y = (y_1,\ldots, y_H)'$. Let $\boldsymbol X = (X_1,\ldots, X_H)$ be the forecast random vector for $\boldsymbol Y$. In practice we have to think about the reporting of these forecasts (see Section (ref)). Here, we simply assume that we have the full multivariate (cumulative) distribution function $\boldsymbol F_{\boldsymbol X}$ as reported forecast available.
Given a forecast $\boldsymbol F_{\boldsymbol X}$ and a random variable $\boldsymbol Y$ we can define a scoring rule $\text{S}(\boldsymbol F_{\boldsymbol X}, \boldsymbol Y)$ that maps into ${\mathbb R}$ (or $\ov{{\mathbb R}}$). They are designed in such a way that a good forecast yields small values (positively orientated), that is if $\boldsymbol F_{\boldsymbol Y}$ is close (or closely related) to $\boldsymbol F_{\boldsymbol X}$. Such a scoring rule is proper if $\text{S}(\boldsymbol F_{\boldsymbol Y}, \boldsymbol Y) \leq \text{S}(\boldsymbol F_{\boldsymbol X}, \boldsymbol Y) $ holds for arbitrary random vectors $\boldsymbol X$, and strictly proper if equality holds only if $\boldsymbol F_{\boldsymbol X} = \boldsymbol F_{\boldsymbol Y}$. Hence the true distribution can be well separated.
In forecasting applications the random variable $\boldsymbol Y$ is usually observed (or realized / materialized), therefore we often write directly $\text{S}(\boldsymbol F_{\boldsymbol X}, \boldsymbol y)$. Note in some literature, see e.g. gneiting2007strictly, scoring rules are defined with inverse orientation. Obviously, all the theory holds as well. Here, we use the notation that is popular in applications, that means the smaller the score the better the forecast.
As pointed out in the introduction, at the moment there is basically only a limited amount of scoring rules available. Here we recall the the major important ones, as we will work with them in the latter part of this paper. We start with the most popular univariate scoring rule.
The continuous ranked probability score (CRPS) for a random variable $X$ with distribution function $F_{X}$ and given forecast $y$ is defined by
where $\widetilde{X}$ is an iid copy of $X$. It is a special case of the energy score, which we will introduce in the next section, for one dimension and $\beta=1$. It is a strictly proper scoring rule with respect to the distribution of $Y$. Another univariate score is for example the quantile loss or pinball score on a grid of quantiles ($0\leq\alpha_1<\alpha_2<\ldots <\alpha_N\leq 1$). Note that for an equidistant grid for $|\alpha_i - \alpha_{i+1}| \to 0$ the corresponding score converges to the CRPS, see e.g. nowotarski2017recent.
If $Y$ respectively $X$ has a density, the logarithmic score is also a candidate for univariate scoring. However, in practical applications we have the problem that for sophisticated models we do not have an explicit formula for the forecasted density available. Still, it may be approximated lerch2017forecaster.
For more details on univariate scoring rules we recommend gneiting2007strictly. In the subsequent sections, we will introduce some multivariate scores.
For the $H$-dimensional random variable $\boldsymbol X$ and the observation vector $\boldsymbol y$ of the target distribution $\boldsymbol Y$ the (Euclidean) energy score (see gneiting2007strictly) is given by
where $\widetilde{\boldsymbol X}$ is an i.i.d. copy of $\boldsymbol X$, so it is drawn independently from the same distribution $\boldsymbol F_{\boldsymbol X}$ as $\boldsymbol X$. Moreover, $\beta \in (0,2)$ and $\|\cdot\|_2$ is the the Euclidean norm. As described in gneiting2007strictly, $\text{ES}_\beta ( \boldsymbol F_{\boldsymbol X}, \boldsymbol y)$ is a strictly proper evaluation score for all $\beta$, but $\beta=1$ seems to be the standard choice in application. szekely2013energy points out that for heavy tailed data the choice of small $\beta$ values should be favored to guarantee that the corresponding moments exists.
The second multivariate score of interest is the variogram score (see scheuerer2015variogram) which is motivated by the variogram that is a popular tool in geostatistics, see cressie1985fitting. The variogram score is defined by
with $p>0$ and weight matrix $\boldsymbol W = (w_{i,j})_{i,j}$. In applications standard cases for $p$ are $p=0.5$ and $p=1$, the weight matrix is usually chosen with $w_{i,j}=1$. It is designed to capture differences in the correlation structure.
Another alternative is motivated by the multivariate logarithmic score (or log-score) which is available if $\boldsymbol X$ has a density $\boldsymbol f_{\boldsymbol X}$. It is given by $$ \text{LogS}( \boldsymbol F_{\boldsymbol X} , \boldsymbol y) = \log(\boldsymbol f_{\boldsymbol X}( \boldsymbol y )) .$$ If $\boldsymbol Y$ has a density, this score is strictly proper as well. Moreover, it has locality properties which might be an advantage, depending on the application. As mentioned in the introduction, in many applications a density forecast is not available. Therefore the multivariate log-score has a limited range of applications, even though the density might be approximated, see lerch2017forecaster.
Under the normality assumption for $\boldsymbol X$ (up to some dropped constants) the log-score yields the Dawid-Sebastiani score, which is given by
with $\boldsymbol \mu_{\boldsymbol X}$ and $\boldsymbol \Sigma_{\boldsymbol X}$ as mean and covariance matrix of $\boldsymbol X$ and $\det$ as determinant. It is clear this score is only proper but not strictly proper as it only matches the mean and the covariance matrix of the forecast with the observations.
After we introduced the most well known and adapted scoring rules, we present a new scoring approach using copulas in the subsequent part. The idea is to describe a given multivariate distribution by its marginals and the copula for the dependency structure. Hence, we have to combine a scoring rule for the marginals with one for the copula. We start with the description of the marginal score.
Let $\text{MS}_h$ be a scoring rule for the random variable $Y_h$, the $h$-th marginal score. $\text{MS}_h$ is a function depending on the forecast $X_h$ of $Y_h$ and the observation $y_h$, so precisely $\text{MS}_h = \text{MS}_h(F_{X_h}, y_h)$ with $F_{X_h}$ as forecasted (cumulative) distribution function of $X_h$, which is the $h$-th marginal distribution of $\boldsymbol X$. Following gneiting2007strictly, we could use every plausible scoring rule for univariate distributions, as the aforementioned CRPS. We construct a marginal score as follows. Let $\text{\textbf{MS}} = (\text{MS}_1, \ldots, \text{MS}_H)'$ denote the vector of marginal scores for each dimension and $\boldsymbol a = (a_1,\ldots, a_H)'$ be a coefficient vector with $a_h>0$, then we define the (joint) marginal score by $$\text{MS}(\boldsymbol a) = \boldsymbol a' \text{\textbf{MS}} = \sum_{h=1}^H a_h \text{MS}_h .$$ It is obvious, that if $\text{MS}_h$ is a (strictly) proper scoring rule for $Y_h$ then the resulting marginal score $\text{MS}(\boldsymbol a)$ is (strictly) proper for $\boldsymbol Y$ for any choice of $\boldsymbol a$. However, the most intuitive choice is clearly $\boldsymbol a = \boldsymbol 1_H$ which equally weights the marginals. This is usually a plausible choice as all the $i$-th step ahead forecasts live on the same scale. In many applications $\boldsymbol a = \boldsymbol 1_H /H$ or $\boldsymbol a = \boldsymbol 1_H$ is used by default, e.g. in the Global Energy Forecasting Competition (GEFCom2014) (hong2016probabilistic). Still, other weighting schemes for evaluation the marginals are plausible as well. Any marginal score can correctly identify the true marginal distribution. However, the major drawback of those scores is that they are not feasible for identification of the correlation structure of the distribution. In order to account for that, we combine the marginal score with another score for the dependency structure. Here the copula comes into play.
First, we briefly recall some copula theory. The core of all copula theory is Sklar's theorem. It states that any random vector $\boldsymbol Y = (Y_1,\ldots,Y_H)'$ can be represented by its marginal distributions $F_{Y_1}, \ldots, F_{Y_H}$ and the copula function $\boldsymbol C_{\boldsymbol Y}$ defined on a $H$-dimensional unit cube. Vice versa, given marginal distribution and a copula function, they describe the joint distribution function. If the marginal distributions $F_{Y_1}, \ldots, F_{Y_H}$ are continuous this representation is unique.
As $F_{Y_1}, \ldots, F_{Y_H}$ characterize the marginal distributions, the copula $\boldsymbol C_{\boldsymbol Y}$ represents the full dependency structure of $\boldsymbol Y$. We want make use of this characteristic, as e.g. the energy score shows weaknesses in discriminating well between dependency structures. Further, remember that the copula $\boldsymbol C_{\boldsymbol Y}$ is a (cumulative) distribution function - so we can apply all characteristics to them.
Let $\text{CS}$ denote a copula score, i.e. an evaluation score for the copula $\boldsymbol C_{\boldsymbol Y}$, which is a multivariate (cumulative) distribution function. If the $\text{CS}$ is a strictly proper scoring rule it can identify the true dependency structure of $\boldsymbol Y$. In the forecasting evaluation literature there is no specific focus on the copulas, even though their are often used to link marginal distribution forecast with a copula forecast to represent the full multivariate distribution or adjust ensemble forecasts (e.g. moller2013multivariate, madadgar2014towards). However, as every copula is just a special multivariate distribution we can use the forecasting evaluation and scoring rule literature on multivariate forecasting evaluation to construct copula scores. Essentially we have three possible options, the energy score, the variogram score and the Dawid-Sebastiani score. Note that only the former one is a strictly proper scoring rule and thus ad-hoc preferable. Still, we will discuss all options within the paper.
Considering that we have the marginal score $\text{MS}$ and the copula score $\text{CS}$ we have to link both to create a proper scoring rule. Therefore, we need a transformation $g:{\mathbb R}\times{\mathbb R} \to {\mathbb R}$ that preserves the properties of $\text{MS}$ and $\text{CS}$. It is easy to observe that the functions $g$ should satisfy some monotonicity conditions, e.g. $g$ is strictly monotone in both components. In more detail, we need that $g$ is strictly isotonic (also known as $2$-monotonic) on the supports. Thus, we require $g$ such that $g(x_1,y_1) - g(x_1,y_2) - g(x_2,y_1) + g(x_2,y_2) > 0 $ holds for $x_1,x_2\in\text{supp}(\text{MS})$ and $y_1,y_2\in\text{supp}(\text{CS})$ with $x_1<x_2$ and $y_1<y_2$. The choice $g(x,y) = x+y$ is strictly isotonic on ${\mathbb R}\times {\mathbb R}$, while $g(x,y)=xy$ works on $(0,\infty)\times(0,\infty)$. Other options would be feasible as well, but for practical application a simple transformation $g$ is favored.
The construction of a combined score is done as follows: Let $\text{MS}(\boldsymbol a)$ be the marginal score and $\text{CS}$ be the copula score, then we propose the marginal-copula score by
Hence, formally we choose $g(\text{MS}(\boldsymbol a),\text{CS}) = \text{MS}(\boldsymbol a) \cdot \text{CS} $ which is a multiplicative structure.
Note that the additive option $g(\text{MS}, \text{CS}) = \lambda \text{MS} + (1-\lambda)\text{CS}$ with $\lambda\in (0,1)$ looks appealing at the beginning as well. However, it turns out to be impractical in application as $\text{MS}$ and $\text{CS}$ live on different scales. Consider for example a forecast for $\boldsymbol Y$ and $c \boldsymbol Y$ with $c>0$. For most marginal scores it follows that $\text{MS}(F_{\boldsymbol Y}, \boldsymbol Y ; \boldsymbol a) \neq \text{MS}(F_{c\boldsymbol Y}, c\boldsymbol Y ; \boldsymbol a)$. For the popular CRPS we even have $\text{MS}(F_{\boldsymbol Y}, \boldsymbol Y ; \boldsymbol a) = \frac{1}{c} \text{MS}(F_{c\boldsymbol Y}, c\boldsymbol Y ; \boldsymbol a)$. On the other hand, for the copula it holds $\text{CS}(\boldsymbol C_{\boldsymbol Y}, \boldsymbol U_{\boldsymbol Y} ) = \text{CS}(\boldsymbol C_{c\boldsymbol Y}, \boldsymbol U_{c\boldsymbol Y} )$ with $\boldsymbol U_{\boldsymbol Y} \sim \boldsymbol C_{\boldsymbol Y}$ and $\boldsymbol U_{c\boldsymbol Y} \sim \boldsymbol C_{c\boldsymbol Y}$. So the marginal score changes with its scale but the copula does not change at all. Thus the additive approach does not seem to be practical for empirical application. One could solve this problem with an adequate scaling for the marginal score but since in general practice the real distribution is unknown it is difficult to find proper bounds.
Moreover, the multiplicative structure of (ref) can be justified by looking at the 2-dimensional normal distribution as in the energy score study in pinson2013discrimination. For the bivariate normal distribution of $\boldsymbol Y = (Y_1,Y_2)'$ only the correlation $\rho$ determines the copula. If $Y_1$ and $Y_2$ have a certain unit (e.g. $\$$, $kWh$ or $m^2/s$), then $\mu_1$, $\mu_2$, $\sigma_1$ and $\sigma_2$ have the same unit, but $\rho$ has no unit at all. Usually the marginal scores inherits the units (e.g. $\$$, $kWh$ or $m^2/s$) but the copula score can never have a unit. Therefore the definition $c MS + (1-c)CS$ provides interpretation troubles.
Here we want to summarize everything discussed above using the corollary that follows directly from the theorem of Sklar:
Unfortunately, this holds only for continuous random variables. In case of non continuous marginals the theorem holds for normal (but not strict) propriety at least. Still, many applications are covered.
Another feature that we want to point out is that the copula is invariant against strictly monotonic transformations. So if we have strictly monotonic transformations $g_h$ then the copula of $(g_1(Y_1), \ldots,g_H(Y_H))'$ has the same copula as $\boldsymbol Y = (Y_1,\ldots, Y_H)'$. Thus, only the marginal score components $\text{MS}_h$ get affected by strictly monotonic transformations; especially the scoring rule remains strictly proper after monotonic transformations of the random variables $\boldsymbol Y$ of interest.
In the next part, we discuss in more details possible choices for the copula score $\text{CS}$. We discuss the energy score, the variogram score and the Dawid-Sebastiani score applied to copulas. Remember that latter score evaluates only the first two moments. As all copulas have uniform marginals, they have all the same mean and variance structure, thus only the correlation structure would be evaluated. The log-score is not suitable for application as we would require a copula density forecast which is rarely available in practice.
We are interested in evaluating the fit of the copula $\boldsymbol C_{\boldsymbol X}$ belonging to the forecast $\boldsymbol F_{\boldsymbol X}$ to the copula $\boldsymbol C_{\boldsymbol Y}$ of $\boldsymbol Y$. Usually we require for forecasting evaluation the observed value $\boldsymbol y=(y_1,\ldots, y_H)'$ of $\boldsymbol Y$. In the copula case we apply the probability integral transformation of each component to receive $$\boldsymbol U_{\boldsymbol Y} = (U_{\boldsymbol Y,1}, \ldots, U_{\boldsymbol Y,H})' = (F_{Y_1}(Y_1), \ldots, F_{Y_H}(Y_H))' = \boldsymbol F_{\boldsymbol Y} (\boldsymbol Y) $$ which has uniform marginals and define the (pseudo) copula observations as $$\boldsymbol u_{\boldsymbol Y} = (u_{\boldsymbol Y,1}, \ldots, u_{\boldsymbol Y,H})' = (F_{Y_1}(y_1), \ldots, F_{Y_H}(y_H))' = \boldsymbol F_{\boldsymbol Y} (\boldsymbol y) .$$ Note that in application we can never observe $\boldsymbol u_{\boldsymbol Y}$ as the required marginal distributions $F_{Y_h}$ are unknown. We will further address this issue in section (ref).
Now, we define the copula energy score ($\text{CES}$) as the energy score of the copula scaled by $H^{-\frac12}$
where $\widetilde{\boldsymbol U}_{\boldsymbol X}$ is an iid copy of $\boldsymbol U_{\boldsymbol X}$, so $\boldsymbol U_{\boldsymbol X}, \widetilde{\boldsymbol U}_{\boldsymbol X} \stackrel{\text{iid}}{\sim} \boldsymbol C_{\boldsymbol X}$. The additional $ \text{lb}_{\text{CES}}$ term is inserted because it corresponds to the lower bound where we choose $ \text{lb}_{\text{CES}} = \frac{1}{4}- \frac{1}{2 \sqrt{6}}$ and discussed below. Note that in gneiting2007strictly a more general version of the energy score is introduced which involves an additional parameter $\beta$ which is chosen here (but also in other studies like pinson2013discrimination) to be $1$. Furthermore we consider the energy-score with negative orientation, thus its optimum is a minimum (for more details see gneiting2007strictly). This is in line with many applications.
In (ref) we defined the copula energy score scaled by $H^{-\frac12}$ to adjust for the dependency in the dimension $H$. This allows a better comparison of score values for different forecast horizons, for instance. Clearly, it would be nice to have a copula score which is scaled in such a way that it allows easy interpretation, like the correlation which is a bounded measure and takes values between -1 and 1 and allows for interpretation of the linear dependency structure. However, such a scaling is not feasible for the copula energy score as its outcome will always depend on the true distribution $\boldsymbol F_{\boldsymbol Y}$, especially the lower and upper bounds will change. Still, both the lower an upper bound grow in $H$ always with a rate of $\sqrt{H}$. Hence, we suggest to scale the outcome with $H^{-\frac12}$.
Since we use negatively orientated scoring rules, the lower bound is the score of the optimal forecast. We know that for a given copula $\boldsymbol C_{\boldsymbol X}$ the energy score $\text{CES}( \boldsymbol C_{\boldsymbol X}, \boldsymbol u_{\boldsymbol Y} )$ is minimized if the forecasted distribution $\boldsymbol X$ equals the true distribution $\boldsymbol Y$ since the energy score is strictly proper.
We studied the lower bound of the energy score of the copula. However, the authors were not able to proof a strict lower bound. Additional research is ongoing. Instead, we consider the lower bound $\text{lb}_{\text{CES}} = \frac{1}{4}- \frac{1}{2}\frac{1}{ \sqrt{6}}$ where the $\frac{1}{4}$ corresponds to a lower bound of the first term ${\mathbb E}\left( \| \boldsymbol U_{\boldsymbol X} - \boldsymbol u_{\boldsymbol Y} \|_2 \right)$ in (ref), and the $\frac{1}{ \sqrt{6}}$ to an upper bound of second term ${\mathbb E}\left( \| \boldsymbol U_{\boldsymbol X} - \widetilde{\boldsymbol U}_{\boldsymbol X} \|_2 \right)$. In the Appendix we show upper and lower bounds for both terms in Lemma (ref) and Lemma (ref).
As mentioned, we can apply the variogram score as a copula score $\boldsymbol C_X$ as well. We define the variogram copula score as a scaled variogram score of the copula:
with $\boldsymbol U_{\boldsymbol X} = (U_{\boldsymbol X,1},\ldots,U_{\boldsymbol X,H})' \sim \boldsymbol C_{\boldsymbol X}$, $p>0$ and weight matrix $\boldsymbol W = (w_{i,j})_{i,j}$. The standard cases are $p=0.5$, $1$ and $2$. For the weight weight matrix we assume $w_{i,j}=1$ which gives a scaling coefficient of $\boldsymbol 1' \boldsymbol W \boldsymbol 1= H^2$. Similarly to the energy score the scaling is motivated by the fact to adjust the dependency on the dimension of the resulting score. As for the energy score case we can derive an upper scaling bound
which justifies the scaling constant. Obviously, the lower bound of $\text{CVS}_{p}$ is zero, as for is holds $ \text{VS}_{\boldsymbol W,p}( \boldsymbol M_{H}, \boldsymbol U_{\boldsymbol Y}) =0$ if $\boldsymbol U_{\boldsymbol Y} \sim \boldsymbol M_H $.
As the variogram score is only a proper scoring rule, the resulting MS-$\text{CVS}_{p}$ score can not be strictly proper. However, from practical perspective the use of MS-$\text{CVS}_{p}$ can be more appealing as the disadvantages occur mainly in the marginal distributions. There is no very simple example where a MS-$\text{CVS}_{p}$ score is not strictly proper if MS is strictly proper.
Finally, we can also consider the Dawid-Sebastiani score (DSS). As mentioned above, the DSS can only evalutate the correlation structure in the forecast, so it is only appealing if we are interested in evaluating the linearly dependency structure in the copula of the forecasts, but not recommended in general.
with $\boldsymbol \mu_{\boldsymbol U_{\boldsymbol X}}$ and $\boldsymbol \Sigma_{\boldsymbol U_{\boldsymbol X}}$ as mean and covariance matrix of $\boldsymbol U_{\boldsymbol X}\sim \boldsymbol C_{\boldsymbol X}$. From properties of the uniform distribution it follows that $\boldsymbol \mu_{\boldsymbol U_{\boldsymbol X}} = \frac{1}{2}\boldsymbol 1$ and $\boldsymbol \Sigma_{\boldsymbol U_{\boldsymbol X}} = \boldsymbol S \boldsymbol R_{\boldsymbol U_{\boldsymbol X}} \boldsymbol S$ where $\boldsymbol S = \frac{1}{\sqrt{12}} \boldsymbol I$ and correlation matrix $\boldsymbol R_{\boldsymbol U_{\boldsymbol X}}$. So we receive
as $ \det( \boldsymbol \Sigma_{\boldsymbol U_{\boldsymbol X}}) = 12^{-H} \det( \boldsymbol R_{\boldsymbol U_{\boldsymbol X}} ) $. Hence, due to (ref) a scaling by $\frac{1}{H}$ might be considered for the Copula Dawid-Sebastiani score to adjust for impact of the dimension.
In contrast to the copula energy score and copula variogram score we do not apply any scaling constant. The reason is that DSS is an unbounded score. For instance, when considering a copula with the the constant correlation matrix $\boldsymbol R_{\boldsymbol U_{\boldsymbol X}}(\delta) = (1-\delta) \boldsymbol I + \delta \boldsymbol 1\boldsymbol 1'$ then it holds for the limit $ \lim_{\delta \to 1} \det( \boldsymbol R_{\boldsymbol U_{\boldsymbol X}}(\delta) ) = 0$. Additionally, the second term is a quadratic form and minimal if $\boldsymbol u_{\boldsymbol y} - \frac{1}{2}\boldsymbol 1$. Thus, the copula Dawid-Sebastiani score has no lower bound and can take negative values.
Remember, we want to forecast $\boldsymbol Y$ by $\boldsymbol X$. In statistical communities, often $\widehat{\boldsymbol Y}$ is usually used instead of $\boldsymbol X$ to emphasize relationship to the target $\boldsymbol Y$.
In practice, there are two options for describing the forecasting distribution $\boldsymbol F_{\boldsymbol X}$ of $\boldsymbol X$ or to report this forecast:
If we are going for option (ref), then we have to report the forecasted distribution $\boldsymbol F_{\boldsymbol Y}$ explicitly or an equivalent representation. A popular way of providing the equivalent information is by reporting an estimate $\boldsymbol f_{\boldsymbol X}$ of the joint density $\boldsymbol f_{\boldsymbol Y}$ which exists for (absolutely) continuous distribution. Alternatively, $\boldsymbol F_{\boldsymbol X}$ can be described using copulas, in detail we have to report estimates $f_{X_{1}},\ldots, f_{X_{H}}$ of the marginal distributions $f_{Y_1},\ldots, f_{Y_H}$ of $\boldsymbol f_{\boldsymbol Y}$ and a copula estimate $\boldsymbol C_{\boldsymbol X}$ for the copula $\boldsymbol C_{\boldsymbol Y}$. For our copula-score based evaluation approach this is a very appealing way of reporting, unfortunately not a very common one. The last option we want to mention is to report an estimator $\boldsymbol \varphi_{\boldsymbol X}$ of the characteristic function defined through $\boldsymbol \varphi_{\boldsymbol Y}(\boldsymbol z) = {\mathbb E}(\exp(i\boldsymbol z'\boldsymbol Y))$ which determines uniquely $\boldsymbol F_{\boldsymbol Y}$.
From the evaluation point of view we must be able to draw from the reported distribution on a computer; or solve some characteristics (esp. some expected values) of the reported distribution which is only feasible in limited cases. Unfortunately for sophisticated forecasting models, an explicit description of $\boldsymbol F_{\boldsymbol X}$ (or an equivalent counterpart) is either very challenging or not feasible. Still, option (ref) can be applied in such cases. Therefore option (ref) turned into the standard alternative in application.
In the second reporting option (ref) we provide $M$ independent simulations resp. draws from the underlying model of $\boldsymbol X$ which essentially describes $\boldsymbol F_{\boldsymbol X}$. We denote the paths/trajectories of the ensemble by ${\mathcal{X}} = \boldsymbol X^{(1)},\ldots,\boldsymbol X^{(M)}$, so $\boldsymbol X^{(j)}$ is an independent and identically distributed (iid) copy of $\boldsymbol X$. This methodology is basically a Monte-Carlo simulation, but it also known as ensemble forecasting, ensemble simulation, path simulation, trajectory simulation, scenario simulation or scenario generation (but depending on the community the meaning might be differ slightly as well). Within this paper we assume that the sample size denoted by $M$ is large, so the true underlying distribution function is described well by the given simulated sample. This approximation goes back to the multivariate Glivenko-Cantelli theorem which characterize that the multivariate empirical (cumulative) distribution function is converging almost surely to the drawn distribution.
For comparing forecasts in a scientific way it is not sufficient to do forecasting evaluation based on a single forecast $\boldsymbol Y$ of $\boldsymbol X$. We require multiple forecasts in a (pseudo-)out-of-sample study design to make significance statements later on. Therefore we assume to have $N$ forecasts of $\boldsymbol Y$ available. How these forecasted values are generated is theoretically not relevant as long only information is taken into account that was available at the time point of generating the forecast.
Now, we assume to have a univariate time series $(Y_t)_{t\in {\mathbb Z}}$ that we are interested to forecast. We denote these multiple targets within the forecasting window by $\boldsymbol Y_{ i} = (Y_{T+s_i+1},\ldots, Y_{T+s_i+H})'$ with $s_i< s_{i+1}$ ($1\leq i<N$). So $\boldsymbol Y_{1},\ldots, \boldsymbol Y_{N}$ are ordered in time. The forecast of $\boldsymbol Y_{1},\ldots, \boldsymbol Y_{N}$ are denoted by $\boldsymbol X_{1},\ldots, \boldsymbol X_{N}$ in accordance with the notation of the previous section. Note, that the forecasting periods of $\boldsymbol Y_{i}$ and $\boldsymbol Y_{i+1}$ my overlap, so if e.g. $s_1=0$, $s_2=1$ with $H=24$ we have forecasts $\boldsymbol X_{1} = (X_{T+1},\ldots, X_{T+24}) $ and $\boldsymbol X_{2} = (X_{T+2},\ldots, X_{T+25} )$ of $\boldsymbol Y_{1} = (Y_{T+1},\ldots, Y_{T+24}) $ and $\boldsymbol Y_{2} = (Y_{T+2},\ldots, Y_{T+25} )$. Obviously, this design should be tailored in accordance with the corresponding application. From the statistical point of view is useful to have nice properties of both the forecasts $\boldsymbol X_{i}$ and the true $\boldsymbol Y_{i}$ in the evaluation window like stationarity, ergodicity or finite variance. However, we can only influence the forecasting model, so $\boldsymbol X_{i}$.
For forecasting models which evaluate historic data we recommend the usage of a rolling window forecasting study using $N$ equidistant windows. It is the most suitable way to evaluate forecasts. In general there are several options for the design of a rolling window forecasting study. But they all share together that the in-sample data length is fixed and the in-sample window moves equally distant across the time range. So the window shift are $s_i = K(i-1)$. If $K<H$ then we have overlapping forecasting horizons. If $K\geq H$ the forecasting horizons in the forecasting study do not overlap which is favorable from the statistical point of view - but does not necessarily meet practical applications. Here the special case $K=H$ is usually chosen. A slightly different design of an rolling window forecasting study with $s_i = 0.25 H(i-1)$ and $H=12$ is visualized in Figure (ref) of the real data application section along with forecasting reporting option reporting option (ref). In many practical applications such a design is clearly favorable. Note that for forecasting evaluation expanding windows forecasting studies are not suitable. More details on the design of forecasting studies, can be found in diebold2015comparing.
Finally, we require for the evaluation that we have observations of the full forecasting window of $\boldsymbol Y_{1},\ldots, \boldsymbol Y_{N}$ available. We denote these observations by $\boldsymbol y_{1}, \ldots, \boldsymbol y_{N}$. This is important for the evaluation method as well. So we do not want to compare the forecasted distributions of $\boldsymbol X_{1},\ldots, \boldsymbol X_{N}$ with the true underlying $\boldsymbol Y_{1},\ldots, \boldsymbol Y_{N}$ but with the observations $\boldsymbol y_{1}, \ldots, \boldsymbol y_{N}$. So from the statistical point of view we require a so called one sample testing framework for comparing distributions.
So far we had only a look at the theoretical properties score of a forecasting method. However, as always in statistics, the theoretical scores involve objects that are unknown in practice. Therefore we assume that we are in a forecasting study framework, where we have the observations $\boldsymbol y_1,\ldots,\boldsymbol y_N$ of $\boldsymbol Y_1,\ldots, \boldsymbol Y_N$ such as the ensemble forecasts ${\mathcal{X}}_i = (\boldsymbol X_{i}^{(1)},\ldots, \boldsymbol X_{i}^{(M)})'$ of the forecasting distribution $\boldsymbol X_i$ for $\boldsymbol Y_i$ available.
Subsequently, we give an overview on estimators for the standard evaluation methods, namely the energy score, the variogram score and the Dawid-Sebastiani score. Afterwards we deal with the problem of estimating copula observations and the marginal copula score.
When estimating (ref), we have essentially two options, matching the two possible representations. The most standard form is to utilise the first one. It solves the corresponding integral numerically by replacing the unknown cdf by an ecdf. $\widehat{\text{CRPS}}_{i,h} = \int_{-\infty}^{\infty} (\widehat{F}_{Y_{i,h}}(z) - \mathbbm{1}\{y_h < z\})^2 \,d z $ with standard estimator for $F_{Y_{i,h}}$ is the empirical distribution function (ecdf) of $h$-th coordinate of the ensemble ${\mathcal{X}}_i$. This is given by $\widehat{F}_{Y_{i,h}}(z) = \widehat{F}_{Y_{i,h}}(z; {\mathcal{X}}_i) = \frac{1}{M} \sum_{j=1}^M \mathbbm{1}\{X_{i,h}^{(j)} \leq z\}$.
Alternatively, we can utilize the second term in (ref) and estimate the two expected values. This is e.g. done in taieb2016forecasting in a load forecasting context. It additionally helps for interpretation. We discuss this resulting estimator as a special case of the estimation of the energy score in the subsequent subsection.
Given the scores $\widehat{\text{\textbf{CRPS}}}_{i} = (\widehat{\text{CRPS}}_{i,1},\ldots, \widehat{\text{CRPS}}_{i,H})$ for each marginal distribution we can compute the corresponding marginal score by $ \widehat{\text{CRPS}}_i(\boldsymbol a) = \boldsymbol a'\widehat{\text{\textbf{CRPS}}}_{i}$ with weight vector $\boldsymbol a = (a_1,\ldots, a_H)'$.
First, we can rewrite the energy score definition (ref) as
If reporting option (ref) is used, we might be able to compute the two terms in (ref) explicitly, e.g. under the normality assumption see pinson2013discrimination. If we are able to do use explicit or numeric approximation formulas for the two terms in (ref), we should prefer them. They are usually more precise and much faster than the simulation approaches.
However, in the standard case of reporting option (ref) we have to estimate the first term in $\text{ED}_{\beta,i}$ and the second term $\text{EI}_{\beta,i}$ in (ref) given the $M$ samples $\boldsymbol X_{i}^{(1)},\ldots,\boldsymbol X_{i}^{(M)}$ of the distribution $\boldsymbol F_{\boldsymbol X}$. The estimation of the first term $\text{ED}_{\beta,i}$ is straight forward, we simply estimate the expectation by its sample mean. So we have
The second term $\text{EI}_{\beta,i}$ has multiple plausible options for the estimation. The reason is that the definition in (ref) of $\text{EI}_{\beta,i}$ suggests we require these independent copies $\widetilde{\boldsymbol X}_{i}$ of our forecasting distribution. As $\boldsymbol X_{i}^{(1)},\ldots,\boldsymbol X_{i}^{(M)}$ are independently drawn, we can simply use the first have of the data set as draws from $\boldsymbol X$ in (ref) and the second half as draws from the iid copy $\widetilde{\boldsymbol X}$. Hence we can simply use the estimator
where we assume without loss of generality that $M$ is even. Note that the sum in $\widehat{\text{EI}}^{\text{iid}}_{i,\beta}$ contains only $M/2$ summands. But these summands have nice statistical properties, as they are iid. Also moller2013multivariate and junk2014comparison mention the resulting estimator for the energy score as a plausible option.
Still, we can use an alternative estimator that uses more summands and provides in average higher accuracy. Such an estimator is given by
for an integer $K$ with $1\leq K\leq M$ where we set for notational purpose $\boldsymbol X_{i}^{(M+k)} = \boldsymbol X_{i}^{(k)}$.
The expression $\widehat{\text{EI}}^{K\text{band}}_{i,\beta}$ has the advantage by using a larger amount of summands for approximating the sum. This should increase the precision in general. However, as the elements of the sum become pairwise dependent, this improvement is weaker as if they would be all independent.
For the important case $K=1$ we have $M$ summands across the summation diagonal, for $K=M$ have $M(M-1)/2$ different pairs $\boldsymbol X_{i}^{(j)} - \boldsymbol X_{i}^{(l)}$ that are summed up. The latter cases uses the maximal amount of information for estimating $\text{EI}_{\beta,i}$. In the R packages Rpackage_scoringRules and Rpackage_energy, only this estimator is provided. It has the mentioned advantage of high accuracy but a high computation demand for large values of $M$ as the amount of summands is increasing quadratically in $M$. Thus, if $M$ is large (which is recommended to choose), we suggest to use a small $K$ to keep the computational complexity manageable. Then the energy score is estimated by $$ \widehat{\text{ES}}^{K\text{band}}_{i,\beta} = \widehat{\text{ED}}_{i,\beta} - 0.5 \widehat{\text{EI}}^{K\text{band}}_{i,\beta} .$$
The computation is similarly to the energy score. If we apply the reporting option (ref), we might be able to solve some special cases explicitly, like the case of multivariate normality.
If we consider the simulation reporting option (ref), we have can estimate the expected values in (ref) by the sample means. Thus, we have
for a given $\boldsymbol W$ and $p$. Due to symmetries $y_{i,j}-y_{i,k} = y_{i,k}-y_{i,j}$ and $X_{i,j}-X_{i,k} = X_{i,k}-X_{i,j}$ the computational cost can be halved in the implementation. Thus it holds
The Dawid-Sebastiani score depends only on the first two moments. Under reporting option (ref), we are usually able to compute the corresponding moments explicitly. Under the reporting option (ref), we estimate the moments by their sample counterparts:
Then we have for the mean and covariance matrix estimates
where the latter one is the unbiased sample covariance matrix. The resulting plug-in estimator for the Dawid-Sebastiani score is
Note that we require $M>H$ to receive a positive semidefinite sample covariance matrix.
Remember that for all copula scores we are interested in evaluating $\boldsymbol U_{\boldsymbol Y_i} = (U_{\boldsymbol Y_i,1}, \ldots, U_{\boldsymbol Y_i,H})' = (F_{Y_{i,1}}(Y_{i,1}), \ldots, F_{Y_{i,H}}(Y_{i,H}))' $. There we defined the copula observations $\boldsymbol u_{\boldsymbol Y_i} = (u_{\boldsymbol y_i,1}, \ldots, u_{\boldsymbol y_i,H})' = (F_{Y_{i,1}}(y_{i,1}), \ldots, F_{Y_{i,H}}(y_{i,H}))'$. However, as mentioned already above a crucial problem here is that the definition of $\boldsymbol u_{\boldsymbol Y_i}$ includes the marginals $F_{Y_{i,h}}$ of the true underlying distribution which is unknown in practice. Thus, the only way around this to estimate $F_{Y_{i,h}}$ to get an estimator for $\boldsymbol u_{\boldsymbol Y_i}$.
Of course the standard estimator for $F_{Y_{i,h}}$ is the empirical distribution function (ecdf) of the ensemble ${\mathcal{X}}_i$. However, we can use the standard estimator $\widehat{F}_{Y_{i,h}}(z) = \widehat{F}_{Y_{i,h}}(z; {\mathcal{X}}_i) = \frac{1}{M} \sum_{j=1}^M \mathbbm{1}\{\boldsymbol x_i^{(j)} \leq z\}$.\footnote{ The mid-point rule $\widehat{F}_{Y_{i,h}}^{\text{mid}}(z) = \frac{1}{2M} \sum_{j=1}^M \mathbbm{1}\{\boldsymbol x_i^{(j)} \leq z\} + \mathbbm{1}\{\boldsymbol x_i^{(j)} < z\}$ might be even a better choice.} Then we receive estimated copula observations $$\widehat{\boldsymbol u}_{\boldsymbol Y_i} = \widehat{\boldsymbol u}_{\boldsymbol Y_i}({\mathcal{X}}_i) = (\widehat{F}_{Y_{i,1}}(y_{i,1}; {\mathcal{X}}_i), \ldots, \widehat{F}_{Y_{i,H}}(y_{i,H}; {\mathcal{X}}_i))'.$$
Intuitively, the use of $\widehat{\boldsymbol u}_{\boldsymbol Y_i}$ instead of $\boldsymbol u_{\boldsymbol Y_i}$ should not be a problem if $\widehat{F}_{Y_{i,h}}$ is a good estimator for $F_{Y_{i,1}}$; and here intuition is right. Unfortunately it is not a good choice if this is not the case. This quickly leads to problems concerning the application of $\widehat{\boldsymbol u}_{\boldsymbol Y_i}$ for the evaluation of copula based scores. There are situations with a misspecified forecasting marginal distribution and the use of $\widehat{\boldsymbol u}_{\boldsymbol Y_i}$ as copula observation estimator where the corresponding copula score of the misspecified forecast is smaller than the true forecasting distribution. If for example $X_{i,h}$ of the forecast $\boldsymbol X_i$ has a variance that is too large then $\widehat{u}_{\boldsymbol Y_i,h}({\mathcal{X}}_i)$ takes too many values around the center (usually around 0.5) and too less values in the tails around $0$ and $1$. This is a problem for practical application and restricts the potential applications of $\widehat{\boldsymbol u}_{\boldsymbol Y_i}$ to the comparison of forecasts with the same marginals which is far away from practical needs.
Fortunately, there is a relatively simple way to avoid the problem with the misspecified marginals. We can force the misspecified marginal distribution to be a (plausible) draw from a copula, so that the general monotonic ordering across is preserved, but such that they have uniform marginals. To perform such an adjustment we have to learn about the distribution of the estimated copula observations $\widehat{\boldsymbol u}_{\boldsymbol Y_i}$. Given a single $\widehat{\boldsymbol u}_{\boldsymbol Y_i}$ this is impossible. Therefore, we explore the distributional structure across the full out-of-sample period $\widehat{\boldsymbol u}_{\boldsymbol Y_1},\ldots, \widehat{\boldsymbol u}_{\boldsymbol Y_N}$. The most obvious adjustment is to adjust the copula observations on each marginal $\widehat{u}_{\boldsymbol Y_1,h},\ldots, \widehat{u}_{\boldsymbol Y_N,h}$ so that they follow a perfect uniform distribution on $(0,1)$. This can be easily done by ranking the estimated observations $\widehat{u}_{\boldsymbol Y_1,h},\ldots, \widehat{u}_{\boldsymbol Y_N,h}$ for each $h$. Hence, denote $R_{i,h}$ the rank of $\widehat{u}_{\boldsymbol Y_i,h}$ within $\widehat{u}_{\boldsymbol Y_1,h},\ldots, \widehat{u}_{\boldsymbol Y_N,h}$. Then we define the adjusted estimated copula observations by
Note that for the ranks $R_{i,h}$ there can appear ties due to the stepwise structure of the ecdf. To preserve the optimal distribution, we highly suggest break the ties in the ranks at random (uniform sampling). This ranking procedure guarantees that $\widehat{u}^{*}_{\boldsymbol Y_1,h},\ldots, \widehat{u}^{*}_{\boldsymbol Y_N,h}$ have perfect uniform marginals, i.e. they take the values $\frac{1}{2N},\ldots, \frac{2N-1}{2N}$. Perfect means here that for these values the corresponding ecdf minimizes the Kolmogorov-Smirnov distance (but also the L\'evy distance or L\'evy metric which characterizes convergence in distribution) with respect to the uniform distribution. This procedure assumes that across the sample the comonotonicity within the implied copula structure is preserved.
For the forecasted copula $\boldsymbol C_{\boldsymbol X_i}$ of $\boldsymbol X_i$ we simply consider the empirical copula of the ensemble ${\mathcal{X}}_i$ with elements $\boldsymbol x_i^{(j)} = (x_{i,1}^{(j)},\ldots, x_{i,H}^{(j)})'$ as a corresponding estimator. Here we suggest to consider
with the ranks $$ \widetilde{R}_{i,j,h} = \frac{1}{2} \sum_{k=1}^M \mathbbm{1}\{ x_{i,h}^{(k)} \leq x_{i,h}^{(j)} \} + \mathbbm{1}\{ x_{i,h}^{(k)} < x_{i,h}^{(j)} \}$$ where we consider the mid-point rule for the ecdf. Similarly as above potential ties should be broken at random (uniform sampling).
Whenever, we want to estimate a marginal copula score (ref) we apply the plug-in principle. Still, remember that (ref) has a multiplicative structure and it holds ${\mathbb E}(X){\mathbb E}(Y) ={\mathbb E}( XY) - \operatorname{{\mathbb C} ov}(X,Y)$. Thus, we get the plug-in estimator
where $\widehat{\text{MS}}_i(\boldsymbol a)$ and $\widehat{\text{CS}}_i$ are the sample means across the marginal and copula score. $\widehat{\sigma}_{\text{MS},\text{CS}} $ is the covariance between $\text{MS}$ and $\text{CS}$ across $i$. Here we assume implicitly some covariance stationarity across the out-of-sample window study between the marginal and copula score.
For applications the estimated marginal scores $\widehat{\text{MS}}_i$ can be estimated by the CRPS and the copula score can be either the estimated energy score $\widehat{\text{ES}}_{i,\beta}$, the estimated variogram score $\widehat{\text{VS}}_{i,\boldsymbol W, p}$ or the estimated Dawid-Sebastiani score $\widehat{\text{DSS}}_{i}$ of the copula observations and the copula of the forecast. Here, we suggest the estimates $\widehat{u}^{*}_{\boldsymbol Y_i,h}$ (eqn. (ref)) and $\widehat{\boldsymbol C}_{\boldsymbol X_i}$ (eqn. (ref)) for the copula observations and the copula as derived in the previous section. As the Copula Dawid-Sebastiani score $\widehat{\text{DSS}}_{i}$ is not bonded from below, the resulting scoring rule is not strictly proper. In fact, it results in a useless rule for applications.
Let $\text{SC}_i$ be a score of the $i$-th forecasting experiment for $i=1,\ldots,N$ of the considered forecasting model. Further let $\text{SC}^*$ be the score of the true model, so $\boldsymbol Y \sim \boldsymbol F_{\boldsymbol X}$ with draws $\boldsymbol y$.
All covered (strictly) proper scoring rules are negatively oriented. Thus, it holds: the smaller the score the better the forecasting performance. Obviously, for practical application the corresponding sample mean
will be the most relevant criterion when comparing the predictive performance of two forecasts.
In this section we shortly introduce two further evaluation methods, based on the scores. First, we consider the relative change in scores as used for deriving conclusions in pinson2013discrimination. Secondly, we consider the Diebold-Mariano test, which allows for significance statements.
In pinson2013discrimination the relative change in the score with respect to the best forecast is considered:
where $\ov{\text{SC}}^* = 1/N \sum_{i=1}^N \text{SC}^*_i$. Obviously, it holds $\text{RelCh}(\text{SC}^*) = 0$. The idea is to measure the sensitivity in the scores with respect to some biased non-optimal forecast in a relative manner.
The Diebold-Mariano (DM) test was designed for the point forecast evaluation in diebold1995comparing. However, as pointed out in e.g. diebold2015comparing the test design is very general and allows for several generalizations. For instance, in ziel2018day anduniejewski2018variance it was applied in multivariate settings. In moller2015spatially it is considered in an application for the energy score.
The DM-test requires a (pseudo-)out-of-sample forecasting study which aims to forecast $\boldsymbol Y_{1},\ldots,\boldsymbol Y_{N}$ by $\boldsymbol X_{1},\ldots,\boldsymbol X_{N}$ as described in section (ref). Its target is to compare the forecasting accuracy of two forecasts ${\mathbb A}$ and ${\mathbb B}$ on the same (pseudo-)out-of-sample environment based on a scoring rule $\text{SC}$.
The DM-test checks whether the scores $\ov{\text{SC}}^{{\mathbb A}}$ is significantly different from $\ov{\text{SC}}^{{\mathbb B}}$, see (ref). Therefore, score differences are required. Given the losses $\text{SC}_1^{{\mathbb A}}, \ldots, \text{SC}_N^{{\mathbb A}}$ and $\text{SC}_1^{{\mathbb B}}, \ldots, \text{SC}_N^{{\mathbb B}}$, we define the loss differences by $$\Delta_i^{{\mathbb A},{\mathbb B}} = \text{SC}_i^{{\mathbb A}} - \text{SC}_i^{{\mathbb B}}.$$ Now, the key idea of the DM-test is to check if the mean loss (or score) difference $$\ov{\Delta}^{{\mathbb A},{\mathbb B}} = \frac{1}{N} \sum_{i=1}^N \Delta^{{\mathbb A},{\mathbb B}}_i$$ is significantly different from zero or not.
Remember, scoring rules $\text{SC}_i$ have a certain direction, so that either a positive loss or a negative loss represents a better forecasting accuracy. As we consider only negatively oriented scores, a $\ov{\Delta}^{{\mathbb A},{\mathbb B}}$ which is significantly smaller than zero lead to the conclusion that ${\mathbb A}$ is significantly better than ${\mathbb B}$ with respect to the considered scoring rule.
In general, the distribution of $\boldsymbol X_{i}$ is different from $\boldsymbol X_{j}$ for $i\neq j$, in the consequence of the distribution of $\Delta_i^{{\mathbb A},{\mathbb B}}$ is usually different from $\Delta_j^{{\mathbb A},{\mathbb B}}$ as well. So $\ov{\Delta}^{{\mathbb A},{\mathbb B}}$ is the sum of $N$ different distributions. Moreover, $\Delta_i^{{\mathbb A},{\mathbb B}}$ and $\Delta_j^{{\mathbb A},{\mathbb B}}$ exhibit a specific dependency structure. Usually they are not even independent or uncorrelated. Hence, from the statistical point of view it is not possible to make further statements about $\ov{\Delta}^{{\mathbb A},{\mathbb B}}$ unless we make some assumptions for the sequence $(\Delta_i^{{\mathbb A},{\mathbb B}})_{i\in {\mathbb Z}}$ (or $\Delta_1^{{\mathbb A},{\mathbb B}}, \ldots, \Delta_N^{{\mathbb A},{\mathbb B}}$). The standard assumptions in the DM-test are:
Thus, $(\Delta_i^{{\mathbb A},{\mathbb B}})_{i\in {\mathbb Z}}$ is a weakly stationary (or covariance stationary) process. These assumptions are very restrictive, and can be relaxed substantially. Still, under these assumption given above it possible to show that $\ov{\Delta}^{{\mathbb A},{\mathbb B}}$ is asymptotic normal diebold2015comparing, i.e. it holds for $N\to \infty$ that
in distribution for where $\sigma(\ov{\Delta}^{{\mathbb A},{\mathbb B}}) = \sqrt{\gamma_{{\mathbb A},{\mathbb B}}(0) }$ is the standard deviation of $\ov{\Delta}^{{\mathbb A},{\mathbb B}}$. This allows the creation of asymptotic tests, namely the Diebold-Mariano test with the null hypothesis $\text{H}_0: \mu_{{\mathbb A},{\mathbb B}} = 0$. When replacing $\sigma(\ov{\Delta}^{{\mathbb A},{\mathbb B}})$ with a suitable estimator, e.g. the sample standard deviation $\widehat{\sigma}(\ov{\Delta}^{{\mathbb A},{\mathbb B}})$ the same distributional limit is attained.
Note that the DM-test is a asymptotic test, as a consequence the out-of-sample window length $N$ should be as large as possible. In the next section we also study in more detail the impact of the simulation sample size $M$ of the proposed evaluation procedure.
In the subsequent subsections we will carry out some simulation studies to compare the forecasting measure in various forecasting situations. These experiment are designed with $N$ replications (or instances) indexed by $i$ as described in the reporting section. For evaluation we will consider always 9 different multivariate measures:
The first three are the established multivariate forecasting criterion. The next two are marginal-copula scores with the CRPS as marginal score. The sixth score is the CRPS itself, which evaluated only the marginals. The latter three ones only evaluate the dependency structure using the copulas.
Here we replicate the study of pinson2013discrimination and enlarge the study design to some extend. This study was the reason for pinson2013discrimination to conclude that the energy score has bad discrimination abilities with respect to the dependency behavior. We consider a bivariate normal distribution with zero mean and variance of $1$ in each component. So the true model is $\boldsymbol Y \sim {\mathcal{N}}_2( \boldsymbol \mu, \boldsymbol \Sigma(\rho) )$ with $\boldsymbol \mu= (0,0)'$ and $\boldsymbol \Sigma(\rho) = \left(
\right)$ with correlation $\rho$. In the simulation setup we choose the true $\rho$ on a grid $\{-1,-0.8, \ldots, 0.8, 1\}$. For the forecasting model we consider $\boldsymbol Y \sim {\mathcal{N}}_2( \boldsymbol \mu, \boldsymbol \Sigma(\varrho) )$. Here we choose the forecasted correlation $\varrho$ out of a slightly denser grid $\{-1,-0.9, -0.8 \ldots, 1\}$.
In the experiment we draw $M=2^{14}=16384$ times for a single bivariate forecast. The window length is $N=2^9=512$. Then we compute the relative change in the score and the DM-test statistic with respect to the true model. Note that pinson2013discrimination evaluated only the relative change. The results are given in Figures (ref) and (ref).