[cite,parencite]nameyeardelim,  textcite

#  D-Vine GAM Copula based Quantile Regression with Application to Ensemble
Postprocessing

David Jobst Corresponding author, University of Hildesheim, Institute of
Mathematics and Applied Informatics, Samelsonplatz 1, 31141 Hildesheim,
Germany<EMAIL_ADDRESS>0000-0002-2014-3530 Annette Möller Bielefeld
University, Faculty of Business Administration and Economics,
Universitätsstraße 25, 33615 Bielefeld, Germany, annette.moeller@uni-
bielefeld.deHelmholtz Centre for Infection Research (HZI), Inhoffenstraße 7,
38124 Braunschweig, Germany<EMAIL_ADDRESS>0000-0001-9386-1691 Jürgen Groß University of Hildesheim, Institute of
Mathematics and Applied Informatics, Samelsonplatz 1, 31141 Hildesheim,
Germany<EMAIL_ADDRESS>0000-0002-3861-4708

###### Abstract

Temporal, spatial or spatio-temporal probabilistic models are frequently used
for weather forecasting. The D-vine (drawable vine) copula quantile regression
(DVQR) is a powerful tool for this application field, as it can automatically
select important predictor variables from a large set and is able to model
complex nonlinear relationships among them. However, the current DVQR does not
always explicitly and economically allow to account for additional covariate
effects, e.g. temporal or spatio-temporal information. Consequently, we
propose an extension of the current DVQR, where we parametrize the bivariate
copulas in the D-vine copula through Kendall’s $\tau$ which can be linked to
additional covariates. The parametrization of the correlation parameter allows
generalized additive models (GAMs) and spline smoothing to detect potentially
hidden covariate effects. The new method is called GAM-DVQR, and its
performance is illustrated in a case study for the postprocessing of 2
$\mathrm{m}$ surface temperature forecasts. We investigate a constant as well
as a time-dependent Kendall’s $\tau$. The GAM-DVQR models are compared to the
benchmark methods Ensemble Model Output Statistics (EMOS), its gradient-
boosted extension (EMOS-GB) and basic DVQR. The results indicate that the GAM-
DVQR models are able to identify time-dependent correlations as well as
relevant predictor variables and significantly outperform the state-of-the-art
methods EMOS and EMOS-GB. Furthermore, the introduced parameterization allows
using a static training period for GAM-DVQR, yielding a more sustainable model
estimation in comparison to DVQR using a sliding training window. Finally, we
give an outlook of further applications and extensions of the GAM-DVQR model.
To complement this article, our method is accompanied by an R-package called
gamvinereg on GitHub.

Keywords: conditional copula; vine copula; dependence modeling; quantile
regression; covariate effects; ensemble postprocessing; probabilistic
forecasting.

##  1 Introduction

Nowadays, weather forecasts are based on numerical weather prediction (NWP)
models which suffer from various uncertainties. In practice, ensemble
prediction systems (EPS) are commonly used to address these uncertainties.
Therefore, the NWP model is run multiple times with different model and/or
initial and boundary conditions [20, 34]. Afterwards a forecast ensemble is
generated, which can be seen as a probabilistic forecast allowing to quantify
forecast uncertainty [44].

However, the forecast ensemble usually suffers from biases and dispersion
errors and thus may benefit from statistical postprocessing using past data to
improve calibration and forecast skill. A popular postprocessing model is the
so called Ensemble Model Output Statistics (EMOS, [20]). This method is used
to obtain a full predictive distribution from the ensemble forecasts.
Originally, EMOS was developed for Gaussian distributed weather quantities,
e.g. temperature or air pressure. Later, machine learning methods such as
quantile regression forests (QRF, [55]) or gradient boosting EMOS (EMOS-GB,
[38]) have been investigated to extend the classical EMOS setting. [48]
compared distributional regression networks to QRF as well as EMOS-GB for the
postprocessing of 2 $\mathrm{m}$ surface temperature forecasts and found only
minor differences among these methods for longer training periods. Recently,
the D-vine copula based quantile regression (DVQR), which was developed by
[32] and further extended by Tepegjozova2022, Sahin2022, was used by
Moeller2018, Demaeyer2023 for the postprocessing of 2 $\mathrm{m}$ surface
temperature forecasts. [28] used the same method for the postprocessing of 10
$\mathrm{m}$ surface wind speed forecasts. In all three analyses, DVQR showed
comparable or sometimes even better results with respect to its competing
methods.

Reasons for the superior performance of DVQR are manifold. DVQR is a quantile
regression model that overcomes typical issues of quantile regression such as
quantile crossings, transformations, collinearity and the integration of
interactions of variables [32]. In addition, DVQR is able to model complex
nonlinear relationships between the explanatory variables and response while
imposing less restrictive model assumptions. Last but not least, it can
theoretically adapt any distribution shape.

One drawback in the current DVQR is the fact that it is not straightforward to
explicitly include covariate effects, such as temporal effects into the model.
This is one reason for estimating DVQR by sliding training windows, where the
complete model needs to be re-estimated for each prediction time point. This
can become computationally expensive, as the optimal sliding window size is
not known in advance and additionally needs to be determined. In the ensemble
postprocessing context the sliding window size depends on various factors,
such as considered variables, seasons, locations, etc., which should be
ideally taken into account. Therefore, [28] compared different types of
training periods in the DVQR model estimation, and detected that a reduction
in the computational complexity usually comes along with a worse predictive
performance.

In this work, we exactly tackle this problem and allow for arbitrary covariate
effects in the DVQR model, such as temporal, spatial or spatio-temporal ones.
To be more precise, the correlation among two variables according to Kendall’s
$\tau$ can depend on covariates and is subsequently used to calculate the
parameters for the bivariate copulas in the D-vine copula. For this, we
combine the work of Vatter2015, Vatter2018 who introduced parametric bivariate
copulas and later vine copulas depending on covariates with the DVQR proposed
by [32]. As the correlations linked to the copulas are parametrized in terms
of generalized additive models (GAMs, [23, 21]) and smoothing splines our
method will be called GAM-DVQR.

We apply the GAM-DVQR with covariates modeling temporal effect for the
postprocessing of 2 $\mathrm{m}$ surface temperature forecasts at 462
observation stations in Germany. The results show that our suggested method is
able to capture temporal covariate effects and can select important predictor
variables from a potentially large set. Furthermore, the correlation time-
dependent GAM-DVQR models show better results in comparison to the GAM-DVQR
model assuming constant correlations, and are able to significantly outperform
the benchmark methods EMOS and EMOS-GB. Last but not least, due to the use of
a static training period for GAM-DVQR, the model needs to be fitted only once
which is more efficient than estimating DVQR on a sliding window and therefore
makes it attractive for practical and operational use. To the best of the
authors knowledge, GAM-DVQR has not been suggested yet and further analyzed in
an application. Although the presented application of GAM-DVQR is concerned
with the meteorological field, our suggested method is broadly applicable to
various areas where any correlation dependent covariate effects are required
to be integrated.

The rest of the paper is organized as follows: Section 2 introduces the D-vine
copula based quantile regression methods including DVQR and GAM-DVQR. In
Section 3, the data set for our application is described. A brief overview of
the competing ensemble postprocessing methods is given in Section 4. Section 5
provides a short introduction to the commonly used verification measures in
the ensemble postprocessing field. In Section 6, we discuss the results of our
application. We close with a conclusion and outlook in Section 7.

##  2 D-vine copula based quantile regression methods

In this section, we outline the copula method requirements employed by our
postprocessing approaches.

###  2.1 Copulas

Multivariate standard distributions, e.g. the multivariate normal are often
restricted in their marginal behavior, as they assume that all marginals are
of the same type. The application of copulas allows to overcome this problem.
A $p$-dimensional copula $C$ is a multivariate distribution function on
$[0,1]^{p}$. According to Sklar’s theorem [51], for every multivariate
distribution function $F$ of $p$ continuous variables
$\bm{X}:=(X_{1},\ldots,X_{p})\in\mathbb{R}^{p}$ there exists a copula $C$,
such that

$\displaystyle F(x_{1},\ldots,x_{p})=C(F_{1}(x_{1}),\ldots,F_{p}(x_{p})),$
(2.1)

where $F_{j}$ for $j=1,\ldots,p$ are the marginal distribution functions and
$\bm{x}:=(x_{1},\ldots,x_{p})\in\mathbb{R}^{p}$ are the realizations of
$\bm{X}$. If all distribution functions are differentiable, the corresponding
$p$-dimensional joint density function $f$ can be expressed by

$\displaystyle f(x_{1},\ldots,x_{p})=c(F_{1}(x_{1}),\ldots,F_{p}(x_{p}))\cdot
f_{1}(x_{1})\cdots f_{p}(x_{p}),$ (2.2)

where $c$ denotes the copula density function of the copula $C$ and $f_{j}$
for $j=1,\ldots,p$ are the marginal density functions of the variables
$X_{1},\ldots,X_{p}$.

Therefore, every multivariate distribution function $F$ can be decomposed into
a copula $C$ modeling the dependence and its univariate marginal distributions
allowing to construct a wide range of distributions.

###  2.2 D-vine copula

Multivariate copulas such as e.g. the elliptical and the archimedean copulas
are often not adaptable enough, as they usually presume that the variables in
all the pairs have homogeneous dependence structures. [3, 4] extended the
theory about multivariate copulas by developing the so-called pair-copula
construction (PCC), where the joint dependence is build up by only bivariate
copulas using conditioning. As a PCC is not unique, [4] introduced a graphical
structure which is called regular vine. A regular vine consists of a set of
nested trees, where the edges in one tree become the nodes of the subsequent
one. In a regular vine consisting of $p$ variables, the nodes and edges in the
first tree represent the $p$ variables and unconditional dependence for $p-1$
variables, respectively. In the subsequent trees the conditional dependence of
a pair of variables conditioned on the variables they have in common is
modeled. A regular vine copula is obtained by specifying bivariate copulas, so
called pair-copulas, on each edge of the trees.

A D-vine is special class of a regular vine in which each tree is a path, i.e.
all nodes in the graph are connected to at most two others (see Figure 1).

Tree 1Tree 2Tree 30123$c_{0,1}$ $c_{1,2}$ $c_{2,3}$ 0,11,22,3$c_{0,2;1}$
$c_{1,3;2}$ 0,2;11,3;2$c_{0,3;1,2}$ Figure 1: 4-dimensional D-vine tree
structure and corresponding pair-copula densities.

Therefore, a D-vine copula is a regular vine copula, where the tree structure
is a D-vine. A node in a D-vine copula represents a certain variable, while an
edge between a pair of nodes corresponds to the dependence among the variables
associated with the respective nodes expressed by a pair-copula.

For a short overview about vine copulas see e.g. [7] and for a detailed
introduction see e.g. [6].

###  2.3 D-vine copula based quantile regression

D-vine copulas can be used in a univariate or multivariate regression context.
Our focus will be on the univariate setting, where it is possible to derive a
conditional D-vine copula density. For this, the leaf node in the first tree
of a D-vine needs to be the response variable. In the following, we denote $Y$
as response variable with marginal distribution function $F_{Y}$ and the $p$
predictor variables by $X_{1},\ldots,X_{p}$ with marginal distribution
functions $F_{1},\ldots,F_{p}$. The lower case letters of the response and
predictor variables represent the respective realizations. For a D-vine copula
with node order $(0,1,\ldots,p)$ corresponding to the variable order
$(Y,X_{1},\ldots,X_{p})$, $p\geq 2$, the conditional density of $Y$ given
$X_{1},\ldots,X_{p}$ can be obtained by

$\displaystyle f_{0|1,\ldots,p}(y|x_{1},\ldots,x_{p})$
$\displaystyle=\prod\limits_{j=2}^{p}c_{0,j;1,\ldots,j-1}(F_{0|1,\ldots,j-1}(y|x_{1},\ldots,x_{j-1}),F_{j|1,\ldots,j-1}(x_{j}|x_{1},\ldots,x_{j-1}))$
$\displaystyle\qquad\qquad\cdot c_{0,1}(F_{Y}(y),F_{1}(x_{1}))\cdot f_{Y}(y),$
(2.3)

where $F_{0|1,\ldots,j-1}$ and $F_{j|1,\ldots,j-1}$ denote the distribution
functions of the conditional random variables
$Y|X_{1}=x_{1},\ldots,X_{j-1}=x_{j-1}$ and
$X_{j}=x_{j}|X_{1}=x_{1},\ldots,X_{j-1}=x_{j-1}$, respectively, and can be
calculated recursively. Furthermore, $c_{0,j;1,\ldots,j-1}$ denotes the
bivariate copula (pair-copula) density of the bivariate distribution of
$(Y,X_{j})$ given $X_{1}=x_{1},\ldots,X_{j-1}=x_{j-1}$. To allow for easy
estimation, we make the simplifying assumption [53], i.e. we assume that the
pair-copulas of conditional distributions are independent of the values of
variables on which they are conditioned. Nevertheless, the pair-copula
densities of the higher tree levels depend on the conditioning values by its
arguments.

Based on the conditional D-vine copula distribution, [32] introduced a
quantile regression (DVQR). The conditional quantile function for a D-vine
copula with $p$ predictor variables $X_{1},\ldots,X_{p}$ at quantile level
$\alpha\in(0,1)$ can be calculated as

$\displaystyle
F^{-1}_{0|1,\dots,p}(\alpha|x_{1},\dots,x_{p}):=F_{Y}^{-1}\left(C_{0|1,\ldots,p}^{-1}(\alpha|F_{1}(x_{1}),\dots,F_{p}(x_{p}))\right),$
(2.4)

where $F_{Y}^{-1}$ is the inverse marginal distribution of the response
variable $Y$ and $C_{0|1,\ldots,p}^{-1}$ denotes the conditional D-vine copula
quantile function.

#### Estimation Procedure.

The estimation of the D-vine copula quantile regression follows a two-step
procedure called “inference for margins” [30]. Firstly, the marginal
distributions of all variables are estimated. This is necessary for
transforming the raw data of each variable to uniformly distributed data in
$[0,1]$ by the probability integral transformation (PIT). Consequently, we
obtain the realizations $v=F_{Y}(y)$ and $u_{i}=F_{i}(x_{i})$ of the random
variables $V=F_{Y}(Y)$ and $U_{i}=F_{i}(X_{i})$ for $i=1,\ldots,p$. The
marginal distributions can be estimated parametrically or non-parametrically.
Secondly, the conditional copula function $C_{0|1,\ldots,p}$ can be obtained
in a closed form by a composition of so-called $h$-functions associated with
the pair-copulas [29]. This two-step approach is very often preferred over
estimating the marginal distributions and copulas simultaneously, as the joint
estimation may be harder to implement, can be very time consuming or is
sometimes simply infeasible.

#### Order of Variables.

The only remaining question is in which order the PIT transformed variables
$V,U_{1},\ldots,U_{p}$ should appear in the D-vine. By construction, the
transformed response variable $V$ needs to be the leaf node in the first tree
of the D-vine (node 0 in Figure 1). As the order of the predictors
$U_{1},\ldots,U_{p}$ is usually not obviously predetermined, one can select
the most informative predictors and order them according to their predictive
strength. To do so, [32] propose a sequential forward selection approach to
select the most important predictors by improving an (AIC/BIC)-corrected
conditional log-likelihood.

For the demonstration of the DVQR algorithm, we assume, that $k-1$ predictors
have already been selected and the current D-vine has the ordering
$(V,U_{l_{1}},\ldots,U_{l_{k-1}})$, where
$\left\\{l_{1},\ldots,l_{k-1}\right\\}\subset\left\\{1,\ldots,p\right\\}$.
Using each remaining predictor $U_{j}$ with
$j\in\left\\{1,\ldots,p\right\\}\backslash\left\\{l_{1},\ldots,l_{k-1}\right\\}$,
the current D-vine is estimated for $(V,U_{l_{1}},\ldots,U_{l_{k-1}},U_{j})$.
In each step of the DVQR method, the optimal pair-copulas according to the
minimum AIC/BIC-conditional log-likelihood or maximum conditional log-
likelihood are chosen. Having estimated the necessary pair-copulas to extend
the current D-vine for each of the $U_{j}$, we update the model by adding the
variable which yields to the lowest AIC/BIC- or highest conditional log-
likelihood of the model.

###  2.4 D-vine GAM copula based quantile regression

In the original formulation of DVQR, [32] assume parametric pair-copulas,
where the copula parameters are constant. A natural extension of the
parametric pair-copulas includes additional effects of covariates, e.g. in
time and/or space into the copula parameters by modelling them as functions of
such covariates. This statistical tool is called conditional copula, which has
already been discussed earlier by [46] using a fully parametric approach for
the copula parameter estimation. Later, [15] proposed a non-parametric version
and [1] a semi-parametric conditional copula model. [58] were the first to
suggested an alternative approach based on generalized additive models (GAMs,
[23, 21]) and spline smoothing for the copula parameter. While the previous
mentioned approaches are restricted to bivariate copulas only, [59] extend the
idea of conditional copulas to higher dimensions. More precisely, they used
the GAM based bivariate copula framework as suggested by [58] for the
construction of vine copulas.

In our proposed D-vine GAM copula based quantile regression (GAM-DVQR), we use
the GAM based bivariate copulas as suggested by [58] for the D-vine copula
quantile regression to include effects of covariates. The sequential forward
variable selection algorithm of DVQR is used in the same way for GAM-DVQR. The
difference between these two methods is mainly in the estimation of the
bivariate copulas, which will be briefly illustrated in the following.

For a vector of $q$ covariates $\bm{Z}:=(Z_{1},\ldots,Z_{q})\in\mathbb{R}^{q}$
with realizations $\bm{z}:=(z_{1},\ldots,z_{q})\in\mathbb{R}^{q}$, a
parametric form is assumed for the conditional copula densities $c(\cdot,\
\cdot;\eta(\bm{z}))$, where the copula parameter $\eta(\bm{z})$ depends on the
covariates $\bm{z}$. For frequently used copula families, bijective
transformations between the copula parameter $\eta(\bm{z})$ and Kendall’s
$\tau(\bm{z})$ can be derived (see Table 1). Note, that additional exogenous
variables which do not belong to the set of predictor variables in the D-vine
such as, e.g. temporal, spatial or other variables can be chosen as
covariates. Nonetheless, it is theoretically possible to incorporate predictor
variables from the D-vine as covariates as well.

Copulas | $\eta(\bm{z})$ | $\tau(\bm{z})$
---|---|---
Gaussian, Student-$t$ | $\sin\left(\frac{\pi}{2}\tau(\bm{z})\right)$ | $\frac{2}{\pi}\arcsin\left(\eta(\bm{z})\right)$
Gumbel, Gumbel-180∘ | $\frac{1}{1-\tau(\bm{z})}$ | $1-\frac{1}{\eta(\bm{z})}$
Gumbel-90∘, Gumbel-270∘ | $-\frac{1}{1+\tau(\bm{z})}$ | $-1-\frac{1}{\eta(\bm{z})}$
Clayton, Clayton-180∘ | $\frac{2\tau(\bm{z})}{1-\tau(\bm{z})}$ | $\frac{\eta(\bm{z})}{\eta(\bm{z})+2}$
Clayton-90∘, Clayton-270∘ | $\frac{2\tau(\bm{z})}{1+\tau(\bm{z})}$ | $\frac{\eta(\bm{z})}{2-\eta(\bm{z})}$
Table 1: Mappings between the copula parameter and Kendall’s $\tau$. The
degrees represent the amount of rotation of the respective copula, e.g. a
Gumbel-180∘ copula is a Gumbel copula rotated by 180∘ counterclockwise.

Due to the one-to-one mappings between the copula parameters and Kendall’s
$\tau$, we re-parameterize all conditional copulas in the D-vine copula as
functions of the corresponding Kendall’s $\tau$ and write $c(\cdot,\
\cdot;\tau(\bm{z}))$. The modeling of the Kendall’s $\tau$ instead of the
actual copula parameter might seem unnecessary at first sight. However, two
useful properties arise from this approach [58]. Firstly, a dependence measure
such as Kendall’s $\tau$ is easier to interpret than a copula parameter.
Secondly, this approach makes it simpler to compare different types of copula
families, as there exists a natural relationship between the copula parameter
and Kendall’s $\tau$. If the actual copula parameter is modeled, it often
becomes necessary to specify a link function to ensure the predefined range of
the copula parameter. Depending on the copula family, different link functions
need to be selected, resulting in possible miss-specifications of the link
function (see, e.g., [35]) and in a comparison which is not standardized.

Therefore, [58] suggest to model the change in the correlation with respect to
the covariates as

$\displaystyle
g^{-1}(\tau(\bm{u},\bm{v};\bm{\alpha},\bm{s})):=\bm{\alpha}\bm{u}^{T}+\sum\limits_{k=1}^{K}s_{k}(\bm{v}_{k}),$
(2.5)

where $g^{-1}(\tau):=2\cdot\mathrm{artanh}(\tau)$ is the inverse link function
between the GAM and Kendall’s $\tau$ to ensure the parameter range,
$\bm{u}\in\mathbb{R}^{J}$ and $\bm{v}\in\mathbb{R}^{K}$ are subsets of the
covariate $\bm{z}$ or products thereof to consider interactions, and
$\bm{\alpha}\in\mathbb{R}^{J}$ is a vector of parameters for the linear
component. The mappings $s_{k}:\mathbb{S}_{k}\to\mathbb{R}$ are smooth
functions supported on closed intervals $\mathbb{S}_{k}\subset\mathbb{R}$ for
$k=1,\ldots,K$, i.e. $s_{k}\in C^{2}(\mathbb{S}_{k})$ admits a finite-
dimensional basis-quadratic penalty representation such as natural cubic
splines, cyclic cubic splines or tensor product splines. Moreover,
$\bm{s}:=(\bm{s}_{1},\ldots,\bm{s}_{K})\in\mathbb{R}^{M}$ denotes the
parameter vector for the $K$ smooth functions $s_{k}$ with a total of
$M:=\sum_{k=1}^{K}m_{k}$ parameters.

Models as in Equation (2.5) are called partially linear models [22], as they
consist of a linear component $\bm{\alpha}\bm{u}^{T}$ and a non-linear
component $\sum_{k=1}^{K}s_{k}(\bm{v}_{k})$. The maximum penalized log-
likelihood estimates of the parameters $\bm{\alpha},\bm{s}$ are obtained by
iteratively reweighted generalized ridge regression. For a more technical
description of the copula parameter estimation as well as for extensive
simulation studies for the suggested conditional copulas and vine copulas, see
Vatter2015, Vatter2018. Moreover, it should be mentioned that Kendall’s $\tau$
will be modeled by only one single model specified in Equation (2.5) for all
bivariate copulas in the D-vine copula. As stated above, we will use the term
“predictor variables” for the variables in the D-vine and the term
“covariates” for the variables included in Equation (2.5) for modeling the
conditional bivariate copula.

To complement this work, [27] developed an R-package called gamvinereg for the
GAM-DVQR, which is based on the R-package gamCopula by [59] and on the code of
the R-package vinereg by [41].

##  3 Data

To illustrate the capabilities of GAM-DVQR, we present a case study for the
postprocessing of 2 $\mathrm{m}$ surface temperature forecasts initialized at
1200 UTC for a forecast lead time of 24 $\mathrm{h}$. The 2 $\mathrm{m}$
surface temperature observations are provided by [12] with maximal 5% missing
observations at each synoptic observation station between January 2, 2015 to
December 31, 2020, which leads to 462 observation stations (see Figure 2).

Figure 2: Observation stations for 2 $\mathrm{m}$ surface temperature.

The ensemble forecasts are provided by the [13], consisting of $m=50$
perturbed ensemble members. These forecasts are initialized at 1200 UTC on a
grid of 0.25∘ $\times$ 0.25∘. The gridded data is bilinearly interpolated to
the observation stations.

Additionally to the target variable (2 $\mathrm{m}$ surface temperature) we
add several auxiliary ensemble predictor variables, for an overview see Table
2. We calculate the 10 $\mathrm{m}$ surface wind speed by
$\text{ws10m}:=\sqrt{\text{u10m}+\text{v10m}}$ and the 2 $\mathrm{m}$ surface
relative humidity is approximated by
$\mathrm{r2m}:=\exp\left(\frac{17.625\cdot\mathrm{d2m}}{243.04+\mathrm{d2m}}\right)/\exp\left(\frac{17.625\cdot\mathrm{t2m}}{243.04+\mathrm{t2m}}\right)$
according to [2]. Furthermore, we calculate the variable sine- and cosine-
transformed day of the year (doy) abbreviated as $\sin$ and $\cos$ via
$\sin\left(\frac{2\pi\cdot\mathrm{doy}}{365.25}\right)$ and
$\cos\left(\frac{2\pi\cdot\mathrm{doy}}{365.25}\right)$ for
$\mathrm{doy}\in\left\\{1,2,\ldots,366\right\\}$, respectively.

In the following, ensemble forecasts are comprised to their mean and standard
deviation, where for a weather variable $v$,

$\overline{X}_{v}:=\frac{1}{m}\sum\limits_{i=1}^{m}X_{v,i}\quad\text{and}\quad
S_{v}:=\sqrt{\frac{1}{m-1}\sum\limits_{i=1}^{m}(\overline{X}_{v}-X_{v,i})^{2}},$
(3.1)

will denote the ensemble mean and standard deviation, respectively, from an
$m$ member ensemble $X_{v,1},\ldots,X_{v,m}$. The further variables
$w\in\\{\text{doy},\sin,\cos\\}$ will be denoted by $X_{w}$. The corresponding
realizations are indicated by lowercase letters. The response variable 2
$\mathrm{m}$ surface temperature is represented by $Y_{\mathrm{t2m}}$ with
realization $y_{\mathrm{t2m}}$. In the following, we will use the term reduced
variable set for $\overline{X}_{\mathrm{t2m}},S_{\mathrm{t2m}}$ including the
response variable $Y_{\mathrm{t2m}}$. Furthermore, we designate the set of the
mean and standard deviation of the first 10 weather variables in Table 2 incl.
the response variable as extended variable set.

Variable | Description
---|---
t2m | 2 $\mathrm{m}$ surface temperature
d2m | 2 $\mathrm{m}$ surface dewpoint temperature
pr | surface pressure
sr | surface solar radiation
u10m | 10 $\mathrm{m}$ surface $u$-wind speed component
v10m | 10 $\mathrm{m}$ surface $v$-wind speed component
r2m | 2 $\mathrm{m}$ surface relative humidity
tcc | total cloud cover
ws10m | 10 $\mathrm{m}$ surface wind speed
wg10m | 10 $\mathrm{m}$ surface wind gust
doy | day of the year
sin | sine-transformed day of the year
cos | cosine-transformed day of the year
Table 2: Potential predictor variables.

Eventually, we will use the period of 2015-2019 as training set and the
complete year 2020 as independent validation set. The implementation of some
of the methods requires the tuning of specific hyperparameters and selection
of marginal distributions for which the final specifications can be found on
GitHub. To avoid overfitting in the model selection process, we further split
the training set into the period of 2015-2018 which is used for pure training
of the models, while the year 2019 is used for testing in the model selection
process. After finalizing the choice of the most suitable model variant based
on the testing period, the entire training period from 2015-2019 is used to
fit that model for the final evaluation on the validation set. All
computations on the data set [26] will be carried out using the statistical
software R running version 3.6.3 by [47].

##  4 Ensemble postprocessing methods

In this section, we briefly describe the compared ensemble postprocessing
techniques. All methods will be applied locally, i.e. for each station a
separate model is estimated.

###  4.1 Ensemble model output statistics

Ensemble model output statistics (EMOS), also known as non-homogeneous
regression, is a parametric postprocessing method proposed by [20]. This
approach is based on the idea of distributional regression, assuming a
predictive distribution family $\mathcal{D}(\mu,\sigma,\nu,\varphi)$, where
the parameters $\mu$, $\sigma$, $\nu$, and $\varphi$ indicate location, scale,
shape, and degrees of freedom, respectively. There are link functions to
connect parameters with corresponding predictors
$\bm{x}_{\mu},\bm{x}_{\sigma},\bm{x}_{\nu},\bm{x}_{\varphi}$ via
$\mu:=h_{\mu}(\bm{x}_{\mu}),\sigma:=h_{\sigma}(\bm{x}_{\sigma}),\nu:=h_{\nu}(\bm{x}_{\nu}),\varphi:=h_{\varphi}(\bm{x}_{\varphi})$
in order to retain parameter ranges. In this context the predictors are
typically ensemble members and summary statistics thereof. The predictive
distribution is selected with regard to the type of weather quantity, e.g. a
Gaussian normal distribution for the variable 2 $\mathrm{m}$ surface
temperature as suggested by [20]. As the logistic and skewed logistic
distribution as well as the skew normal distribution show only minor
differences with respect to the performance of the Gaussian normal
distribution Gebetsberger2019, Taillardat2021 we assume the latter in the
following, i.e. $Y_{\mathrm{t2m}}\sim\mathcal{N}(\mu,\sigma)$ with location
parameter $\mu\in\mathbb{R}$, scale parameter $\sigma>0$ as well as inverse
link functions $h_{\mu}^{-1}:=\text{id}$, $h_{\sigma}^{-1}:=\log$.

In its basic formulation EMOS uses the reduced variable set with predictors
$\overline{X}_{\mathrm{t2m}}$ and $S_{\mathrm{t2m}}$ and connects the Gaussian
(transformed) distribution parameters to the predictors via the linear
relationships

$\displaystyle\mu:=a_{0}+a_{1}\overline{x}_{\text{t2m}},\quad\log(\sigma):=b_{0}+b_{1}\log(s_{\mathrm{t2m}}).$
(4.1)

The coefficients $a_{0},a_{1},b_{0},b_{1}\in\mathbb{R}$ are estimated e.g. by
a sliding training window. However, to take the strong seasonal periodic
patterns of $Y_{\mathrm{t2m}}$ (e.g. higher values in the summer period, lower
values in the winter period) into account and to facilitate a fair comparison
with the other methods, we add the sine- and cosine-transformed day of the
year $X_{\sin},X_{\cos}$ as further predictors to the equation of both
parameters. Therefore, we assume the conditional predictive distribution

$\displaystyle f(y_{\mathrm{t2m}}|x_{\text{t2m},1},\ldots,$ $\displaystyle
x_{\text{t2m},m},x_{\sin},x_{\cos})\sim\mathcal{N}(\mu,\sigma),$ (4.2)
$\displaystyle\mu:=a_{0}+a_{1}x_{\sin}+a_{2}x_{\cos}+a_{3}\overline{x}_{\text{t2m}},$
$\displaystyle\quad\log(\sigma):=b_{0}+b_{1}x_{\sin}+b_{2}x_{\cos}+b_{3}s_{\mathrm{t2m}},$
(4.3)

with coefficients $a_{i},b_{i}\in\mathbb{R}$ for $i=0,\ldots,3$, as e.g. [24,
8]. Consequently, we incorporate the seasonality by seasonal varying
intercepts as introduced in Equation (4.3) in comparison to the basic
formulation in Equation (4.1).

The coefficients of the parameters in Equation (4.3) are estimated by
optimizing the sum of a proper verification score over the training period
between 2015 and 2018 with the Broyden-Fletcher-Goldfarb-Shannon (BFGS)
algorithm. We investigate the optimization with respect to two different
scores, namely the CRPS (continuous ranked probability score) and the LogS
(logarithmic score). Then, we choose the best performing specification (CRPS
or LogS) according to the mean CRPS over all testing days in 2019 and all
stations. The implementation is based on the R-package crch by [37].

###  4.2 Gradient boosted ensemble model output statistics

[38] suggested an extension of EMOS which allows to select the most relevant
predictor variables $X_{1},\ldots,X_{p}$ for the model by a gradient-boosting
approach (EMOS-GB). This approach is especially useful if the amount of
potential predictor variables is large. Similar to the EMOS model the
conditional normal distribution
$f(y|x_{1},...,x_{p})\sim\mathcal{N}(\mu,\sigma)$ with

$\displaystyle\mu$ $\displaystyle:=a_{0}+a_{1}x_{1}+\ldots+a_{p}x_{p},\quad
a_{0},a_{1},\ldots,a_{p}\in\mathbb{R},$ (4.4) $\displaystyle\log(\sigma)$
$\displaystyle:=b_{0}+b_{1}x_{1}+\ldots+b_{p}x_{p},\quad
b_{0},b_{1},\ldots,b_{p}\in\mathbb{R},$ (4.5)

is assumed. We include the sine- and cosine transformed day of the year
$X_{\sin},X_{\cos}$ into the extended variable set to account for seasonality.
This results in $p=22$ predictor variables (see Table 2) for each distribution
parameter.

The boosting procedure initializes all coefficients for $\mu$ and $\sigma$ at
zero and to iteratively updates only the one which corresponds to the
predictor improving predictive performance the most. Using the gradient of the
loss function, the predictor with the highest correlation to the gradient is
selected and then the corresponding coefficient is updated by taking a step in
the direction of steepest descent of the gradient. This procedure is carried
out until a stopping criterion is reached to avoid overfitting.

The implementation is based on the R-package crch by [37]. We tune the
gradient-boosted EMOS (EMOS-GB) model by grid search with respect to the loss
function (LogS or CRPS), maximum number of boosting iterations (100, 500,
1000, 2000), stopping criterion (AIC, BIC), and step size (0.05, 0.1, 0.2) on
the training data set between 2015 and 2018. Then, we choose the best
performing model version according to the mean CRPS over all testing days in
2019 and stations for the validation period.

###  4.3 D-vine copula based quantile regression

We apply the D-vine copula based quantile regression (DVQR) as proposed by
[32] and further described in Section 2.4, where we use the reduced variable
set by minimizing the BIC-corrected conditional log-likelihood.

To take account of the seasonality of our variables, DVQR is first estimated
on a refined rolling training period [40] with window size
$n\in\left\\{10,15,20,\ldots,100\right\\}$, for which we use the days
$\\{t-n,\ldots,t-2,t-1\\}$ in the year where the forecast day $t$ lives and
the days $\\{t-n,\ldots,t-2,t-1,t,t+1,t+2,\ldots,t+n\\}$ in the previous $k=4$
years. The optimal window size $n$ for the validation period is determined
based on the minimal mean CRPS over all testing days in 2019 and all stations.
Afterwards, DVQR is estimated using $k=5$ years with the estimated optimal
length $n$. Due to very high computational costs of approximately 7 hours to
estimate the models in 2019 for one station and our limited resources of one
CPU with 40 cores and 62.5 GB RAM to determine the optimal window length $n$
we can not investigate DVQR on the extended variable set. This high
computational burden underlines again the need of an alternative approach such
as GAM-DVQR in a higher-dimensional context.

#### Marginal distributions.

The marginal distributions are fitted via kernel density estimates using the
Gaussian kernel.

#### Bivariate copulas.

We allow all bivariate copulas in the R-package vinereg, that is, elliptical
copulas (Gaussian and Student-$t$), archimedean copulas as well as rotated
versions thereof (Clayton, Gumbel, Frank, Joe, BB1, BB6, BB7 and BB8), and the
nonparametric Independence and Transformation Kernel copula (TLL). Elliptical
copulas have an elliptical shape in the contour plot. In Figure 3 for example
we can detect a Gaussian copula between $S_{\mathrm{r2m}}$ and
$S_{\mathrm{wg10m}}$ and a Student-$t$ copula for the variable pair
$S_{\mathrm{ws10m}}$ and $S_{\mathrm{wg10m}}$. The Gaussian copula has no tail
dependence, while the Student-$t$ copula only captures symmetric tail
dependence. Any departures from elliptical shapes may indicate to include non-
Gaussian dependencies. Therefore, archimedean copulas are provided which
exhibit a pear or bone shape in the contour plot allowing to detect lower
and/or upper tail dependence (except of the Frank copula). In Figure 3 we see
a contour shape indicating a Frank copula between
$\overline{X}_{\mathrm{t2m}}$ and $\overline{X}_{\mathrm{sr}}$, a Gumbel
copula between $\overline{X}_{\mathrm{r2m}}$ and $\overline{X}_{\mathrm{tcc}}$
and a Clayton copula between $\overline{X}_{\mathrm{ws10m}}$ and
$S_{\mathrm{wg10m}}$. The nonparametric Independence copula has a circular
shape in the contour plot (see e.g. between $\overline{X}_{\mathrm{r2m}}$ and
$S_{\mathrm{ws10m}}$ in Figure 3) and the Transformation Kernel copula can
approximate any dependence shape (see e.g. between
$\overline{X}_{\mathrm{u10m}}$ and $\overline{X}_{\mathrm{ws10m}}$ in Figure
3). Consequently, we cover lots of possible dependence patterns by this copula
set.
The implementation of the DVQR is based on the R-package vinereg by [41].

Figure 3: Empirical normalized copula contour plots (lower triangle), PIT
histograms (diagonal) and scatterplots including Kendall’s $\tau$ correlation
(upper triangle) for station Munich in the training data.

###  4.4 D-vine GAM copula based quantile regression

For the GAM-DVQR as explained in Section 2.4 we consider two cases: Estimation
of the GAM-DVQR on the reduced variable set and on the extended variable set.

#### Marginal distributions.

In both variable sets, we determine the marginal distributions for each
variable using distributional regression via generalized additive models for
location ($\mu$), scale ($\sigma$) and shape ($\nu$) (GAMLSS, [49]).
Therefore, we assume for all considered weather variables a distribution
$\mathcal{D}(\mu,\sigma,\nu,\varphi)$. As all weather variables show a
seasonal periodic behavior we model the distribution parameters by using the
notation of Section 4.1 via

$\displaystyle h_{\mu}^{-1}(\mu)$
$\displaystyle=a_{0}+a_{1}x_{\sin}+a_{2}x_{\cos},\quad$ $\displaystyle
h_{\sigma}^{-1}(\sigma)=b_{0}+b_{1}x_{\sin}+b_{2}x_{\cos},$ (4.6)
$\displaystyle h_{\nu}^{-1}(\nu)$ $\displaystyle=c_{0},\quad$ $\displaystyle
h_{\varphi}^{-1}(\varphi)=d_{0},$ (4.7)

with real valued coefficients using the sine- and cosine transformed day of
the year $X_{\sin},X_{\cos}$ as linear covariates. To keep the marginal models
simple and to have comparable settings, the parameters $\nu$ and $\varphi$ are
assumed to be constant. For each variable we allow a set of potential
distribution families, from which the best is chosen with respect to the mean
BIC over all stations in the whole training period. Although gamvinereg allows
to select the best performing marginal distribution for each station, we
decided against this procedure for a fairer comparison of the methods and a
more standardized verification. Besides of the actually tested distribution
families, e.g. censored or truncated versions could be additionally analyzed.
However, initial tests showed only small differences to the considered
distribution families. Additionally, [31] outlined in a simulation study that
misspecified margins are only problematic for copulas if they are severely
misspecified. They studied e.g. fitting normal margins when the true margins
are exponential.

#### Bivariate copulas.

The GAM-copula family set for modeling the dependencies consists of the
Gaussian, Student-$t$, double Clayton type I-IV and double Gumbel type I-IV
copula as implemented in the R-package gamCopula by [42]. The double Clayton
and Gumbel copula type consist of additional rotated versions of the Clayton
and Gumbel copula, respectively to cover negative dependence as well. All
predictor variables and therefore copulas will be selected by minimizing the
BIC-corrected conditional log-likelihood as for DVQR. In both, the reduced and
extended variable set, each Kendall’s $\tau$ linked to a pair-copula is
modeled by two different approaches: We assume either a constant or a time-
dependent correlation in Equation (2.5) and link it to the covariates by one
of the following linear models without a non-linear component

$\displaystyle
g^{-1}(\tau(\bm{u},\bm{v};\bm{\alpha},\bm{s}))=\begin{cases}\alpha_{0},&\text{constant
correlation (C)},\\\
\alpha_{0}+\alpha_{1}u_{\sin}+\alpha_{2}u_{\cos},&\text{time-dependent
correlation (T1)},\end{cases}$ (4.8)

where $\alpha_{0},\alpha_{1},\alpha_{2}\in\mathbb{R}$, and covariates
$u_{\sin},u_{\cos}$ denote the sine- and cosine transformed day of the year.
The need for a time-adaptive correlation between the predictor variables is
illustrated in Figure 4, where the empirical Kendall’s $\tau$ correlation as
well as its predictions clearly change over the day of the year. This aspect
will be further outlined in Section 6.

Figure 4: Empirical Kendall’s $\tau$ (purple) using the refined rolling
training period of window size $n=25$, Kendall’s $\tau$ prediction (darkgreen)
and 95% confidence band (lightgreen) using the time-dependent correlation
model (T1) for station Fürstenzell in the training data.

As the amount of predictor variables can become large in the extended variable
set it might not be clear anymore how the temporal correlation can be
appropriately modelled by a linear component. Thus, we additionally
investigate a time-dependent correlation model using a cyclic cubic spline $s$
depending on the covariate day of the year (doy) via

$\displaystyle
g^{-1}(\tau(\bm{u},\bm{v};\bm{\alpha},\bm{s})))=s(v_{\text{doy}}),\quad\text{time-
dependent correlation (T2)}.$ (4.9)

With this time-dependent non-linear correlation model we can take account of
even more flexible (unknown) changes in Kendall’s $\tau$ than by a linear
model. This approach could be beneficial in higher trees of the D-vine as
well.

Consequently, with GAM-DVQR we introduce an ensemble postprocessing model,
which is able to select the most important predictor variables from a large
set, taking account of the temporal correlation changes among the predictor
variables at the same time. As this model is estimated only once on a static
training period, it benefits from a longer consistent training period, while
being more efficient than DVQR using a sliding training window.

Model | Marginal specifications | Correlation specifications | Variable set
---|---|---|---
GAM-DVQR-C | GAMLSS | constant correlation (C) | reduced & extended
GAM-DVQR-T1 | GAMLSS | time-dependent correlation (T1) | reduced & extended
GAM-DVQR-T2 | GAMLSS | time-dependent correlation (T2) | extended

Table 3: Overview of marginal and correlation parameter specifications.

An overview of the selected marginal and correlation parameter specifications
can be found in Table 3. The estimation of GAM-DVQR is carried out by using
the R-package gamvinereg of [27].

##  5 Verification methods

Gneiting2005, Gneiting2007 claim, that the general aim of probabilistic
forecasting is to maximize the sharpness of the predictive distribution
subject to calibration. Calibration refers to the statistical consistency
between the predictive cumulative distribution function (CDF) $F$ and the
associated observation $Y$. Sharpness concerns the spread of the predictive
distribution $F$. The more concentrated the forecast, the sharper the
forecast, and the sharper the better, subject to calibration. In the
following, we present methods to measure calibration and sharpness which are
be used in the subsequent application.

Figure 5: Verification rank histogram of the raw ensemble aggregated over all
stations and time points in the validation period.

#### Visual assessment of calibration.

Dawid1984, Gneiting2007a call a continuous predictive probabilistic forecast
$F$ calibrated if $F(Y)$ is uniformly distributed. A so called probability
integral transform (PIT) histogram can be used as visual tool for the
evaluation of the calibration, where the PIT values are received by evaluating
the predictive CDF $F$ at the validating observations. Any departures from
uniformity of the PIT histogram can indicate that the predictive distribution
$F$ is miscalibrated in some way. A discrete counterpart of the PIT histogram
is the so called verification rank histogram displaying the histogram of ranks
of observations with respect to the corresponding ordered ensemble forecasts
[56]. In the case of a calibrated $m$-member ensemble, the ranks should be
uniformly distributed on the set $\\{1,\ldots,m+1\\}$.

#### Uncertainty quantification.

A further tool for assessing the calibration of a predictive distribution is
the coverage of a $(1-\alpha)\cdot 100\%$ central prediction interval,
$\alpha\in(0,1)$, which is the proportion of validating observations between
the lower and upper $\frac{\alpha}{2}$-quantiles of the predictive
distribution [18]. Assuming a calibrated predictive distribution, then
$(1-\alpha)\cdot 100\%$ of observations should fall within the range of the
central prediction interval. Sharpness of a predictive distribution can be
validated using the width of a $(1-\alpha)\cdot 100\%$ central prediction
interval [18]. Sharper distributions correspond to narrower prediction
intervals. Having a $m$-member forecast ensemble, we use a
$\frac{m-1}{m+1}\cdot 100\%$ central prediction interval corresponding to the
nominal coverage of the raw forecast ensemble and consequently allowing for a
direct comparison of all probabilistic forecasts. The target coverage rate for
an $m=50$ member ensemble is approximately 96.08%.

#### Scoring rules.

Proper scoring rules rate calibration and sharpness properties simultaneously
and thus play important roles in the comparative evaluation and ranking of
competing forecasts [17]. An attractive proper scoring rule in weather
forecasting is the continuous ranked probability score (CRPS, [36]), which is
defined as

$\displaystyle\text{CRPS}(F,y):=\int\limits_{-\infty}^{\infty}(F(z)-\mathds{1}\\{z\geq
y\\})^{2}\,\mathrm{d}z,$ (5.1)

where $F$ is the predictive cumulative distribution function, $y$ is the
true/observed value and $\mathds{1}$ denotes the indicator function. [19] show
that, if $F$ has a finite first moment the CRPS can be approximated by

$\displaystyle\text{CRPS}(F,y)\approx\frac{1}{K}\sum\limits_{k=1}^{K}|z_{k}-y|-\frac{1}{2K^{2}}\sum\limits_{k=1}^{K}\sum\limits_{k^{\prime}=1}^{K}|z_{k}-z_{k^{\prime}}|,$
(5.2)

where $z_{k}:=F^{-1}\left(\frac{k}{K+1}\right)$,
$z_{k^{\prime}}:=F^{-1}\left(\frac{k^{\prime}}{K+1}\right)$ for
$k,k^{\prime}\in\\{1,\dots,K\\}$ and $F^{-1}$ denotes the quantile function of
$F$. The mean CRPS over a set of forecast cases is denoted by
$\overline{\text{CRPS}}$.

In practice, a probabilistic forecast is sometimes reduced to a point forecast
via a statistical summary function such as the mean or median. In this
situation, consistent scoring functions provide useful tools for forecast
evaluation and generate proper scoring rules [16]. For a set of $n$ forecasts
cases we employ

$\text{RMSE}:=\sqrt{\frac{1}{n}\sum\limits_{i=1}^{n}(\text{mean}(F_{i})-y_{i})^{2}}\quad\text{and}\quad\text{MAE}:=\frac{1}{n}\sum\limits_{i=1}^{n}|\text{median}(F_{i})-y_{i}|.$
(5.3)

The relative improvement of a forecast with respect to a given reference
forecast in terms of CRPS can be quantified by the continuous ranked
probability skill score (CRPSS) via

$\displaystyle\text{CRPSS}:=1-\frac{\overline{\text{CRPS}}}{\overline{\text{CRPS}}}_{\text{ref}},$
(5.4)

where $\overline{\text{CRPS}}_{\text{ref}}$ denotes the
$\overline{\text{CRPS}}$ of the reference forecast.

#### Statistical tests to compare predictive performance.

To evaluate the statistical significance of the differences in the forecasts
between two competing postprocessing models, we make use of the Diebold-
Mariano test [11] for the verification score time series of both models
separately at each station. Afterwards we use the Benjamini-Hochberg procedure
[5] suggested by [60] that allows to account for multiple testing regarding
different stations and to control the overall probability of type I error, for
which we choose $\alpha=0.05$ in the subsequent analysis.

For the verification of the methods we use the R-package eppverification by
[25].

#### Visual dependence assessment.

For visually assessing the dependence of two variables $Y$ and $X$ a so called
empirical normalized bivariate contour plot (see, e.g. Figure 3) can be used
by the R-package VineCopula of [43]. This plot is obtained by an approximation
of the copula density $c$ and to visualize the contours of the bivariate
density function

$\displaystyle
d(z_{Y},z_{X}):=c(\Phi(z_{Y}),\Phi(z_{X}))\phi(z_{Y})\phi(z_{X}),$ (5.5)

for the $\Phi^{-1}$ transformed copula data $Z_{Y}:=\Phi^{-1}(F_{Y}(Y))$ and
$Z_{X}:=\Phi^{-1}(F_{X}(X))$. For more details concerning these plots, see,
e.g. [6].

##  6 Results

In the following two subsections the results of the considered methods based
on the reduced end extended variable set are presented and discussed. Note,
that we additionally investigated time series models as marginal distributions
for the GAM-DVQR. However, the respective results turned out to be worse than
the ones presented here so we do not show them.

###  6.1 Reduced variable set

In this setting, we compare the raw ensemble, EMOS, DVQR and GAM-DVQR on the
reduced variable set. Figure 6 shows the PIT histograms for all postprocessing
methods. All methods show improved calibration properties in comparison to the
raw ensemble in Figure 5. However, they are all slightly skewed to the right
causing a small overdispersion more or less in the middle of the PIT
histograms and underdispersion at both ends. This impression is supported by
the values of the coverage score shown in Table 4.

Figure 6: PIT histograms of the considered methods, where the PIT values are
aggregated over all stations and time points in the validation period.

#### Verification scores.

When looking at Table 4, we observe that all methods improve the raw ensemble
with respect to CRPS by around 20%-29%, MAE around 12%-22% and RMSE around
13%-22%. While EMOS yields the lowest MAE and width, GAM-DVQR-T1 yields the
lowest CRPS and RMSE. Furthermore, using the procedure for testing the
significant differences in the performance between two methods (see Section
5), GAM-DVQR-T1 significantly outperforms EMOS at around 8% of all stations
with respect to CRPS. A reason for the better performance of GAM-DVQR-T1 over
EMOS might be that the GAM-DVQR models use non-Gaussian copulas (see Figure
7), while the basic assumption of EMOS is the Gaussian dependence. The
dependence between $Y_{\mathrm{t2m}}$ and $\overline{X}_{\mathrm{t2m}}$ in
Figure 7 can be described by a Student-$t$ copula, while for the relationship
between $\overline{X}_{\mathrm{t2m}}$ and $S_{\mathrm{t2m}}$ a Gumbel-270∘
(e.g. double Gumbel type II) or Clayton-90∘ (e.g. double Clayton type I)
copula could be estimated. A further reason for improved performance in
comparison to DVQR could be traced back to the longer and more consistent
training data which might lead to more stable estimations [33] for GAM-
DVQR-T1, while DVQR is estimated on a sliding window.

Method | CRPS | MAE | RMSE | Coverage | Width
---|---|---|---|---|---
Raw ensemble | 1.017 | 1.255 | 1.730 | 63.071 | 2.922
EMOS | 0.718 | 0.985 | 1.387 | 95.573 | 5.441
DVQR | 0.717 | 0.985 | 1.366 | 96.339 | 5.990
GAM-DVQR-C | 0.719 | 0.989 | 1.358 | 96.499 | 5.994
GAM-DVQR-T1 | 0.713 | 0.985 | 1.350 | 96.594 | 5.889
Table 4: Verification scores aggregated over all stations and time points in
the validation period. Bold values represent the best value for each score.
Figure 7: Empirical normalized contour plots (lower triangle), PIT histograms
(diagonal) and scatterplots including Kendall’s $\tau$ correlation (upper
triangle) for station Fürstenzell in the training data.

Model | $\alpha_{1,1,0}$ | $\alpha_{1,1,1}$ | $\alpha_{1,1,2}$ | $\alpha_{1,2,0}$ | $\alpha_{1,2,1}$ | $\alpha_{1,2,2}$ | $\alpha_{2,1,0}$ | $\alpha_{2,1,1}$ | $\alpha_{2,1,2}$
---|---|---|---|---|---|---|---|---|---
GAM-DVQR-T1 | 100 | 12 | 74 | 94 | 16 | 78 | 100 | 29 | 62

Table 5: Significance of the coefficients for the correlation time-dependent
GAM-DVQR models in % after applying the Benjamini-Hochberg procedure to the
$p$-value for each coefficient $\alpha_{i,j,k}$ over all stations.
$\alpha_{i,j,k}$ represents the $k$-th coefficient for the $j$-th bivariate
copula in the $i$-th tree of the D-vine with respect to Equation (4.8).

#### CRPS comparisons.

In the following, we compare the methods with respect to CRPS in more detail.
The time-dependent correlation GAM-DVQR models (T1) perform better than the
GAM-DVQR models with constant correlation (C) in terms of CRPS. These results
underline the need of a time-dependent correlation model within the GAM-DVQR
framework in this application, and also highlight that GAM-DVQR is able to
capture the temporal varying empirical correlation. This is further
illustrated in Figure 4 where the empirical and predicted Kendall’s $\tau$ for
a D-vine GAM copula with order
$Y_{\mathrm{t2m}}-\overline{X}_{\mathrm{t2m}}-S_{\mathrm{t2m}}$ are plotted.
Moreover, the percentage of the significant coefficients for the correlation
Equation (4.8) in Table 5 indicate that the suggested predictors for
identifying the changes in the correlation parameter seem appropriate.

Figure 8: Highest CRPSS of the considered methods over EMOS in % in the
validation period. CRPSS $>5\%$ are visualized in black for a better
representation. Numbers in brackets denote the count.

We further investigated the station-specific performance of our methods with
respect to CRPSS over the benchmark method EMOS. Figure 8 shows the method
with the highest CRPSS over EMOS, where the CRPSS values are encoded by
colours. The DVQR based methods yield at around 68% of all stations a positive
CRPSS (green colour scale), which implies a substantially better performance
of these methods over EMOS. Furthermore, we observe that GAM-DVQR-T1
outperforms the other models most frequently and yields the highest CRPSS at
219 stations, while DVQR is the preferred model at 159 stations.

###  6.2 Extended variable set

In this section we investigate the results for the raw ensemble, EMOS, DVQR
and GAM-DVQR on the extended variable set.

The PIT histograms in Figure 9 indicate that all methods are able to improve
the calibration, where the remaining deficiencies are less pronounced than in
Section 6.1, but still visible. Moreover, GAM-DVQR-T1 and GAM-DVQR-T2 seem to
yield the PIT histograms which are closest to a uniform distribution. This
impression is further supported by the nearly perfect coverage score of
96.079% for GAM-DVQR-T1 in Table 6. All in all, it appears that the GAM-DVQR
models provide a more pronounced calibration in terms of the PIT histograms as
well as better coverage values.

Figure 9: PIT histograms of the considered methods, where the PIT values are aggregated over all stations and time points in the validation period. Method | CRPS | MAE | RMSE | Coverage | Width
---|---|---|---|---|---
Raw ensemble | 1.017 | 1.255 | 1.730 | 63.071 | 2.922
EMOS-GB | 0.706 | 0.979 | 1.336 | 96.756 | 5.670
GAM-DVQR-C | 0.710 | 0.980 | 1.337 | 96.145 | 5.680
GAM-DVQR-T1 | 0.684 | 0.942 | 1.282 | 96.079 | 5.449
GAM-DVQR-T2 | 0.681 | 0.939 | 1.278 | 96.157 | 5.462
Table 6: Verification scores aggregated over all stations and time points in
the validation period. Bold values represent the best value for each score.

#### Verification scores.

All methods are able to improve upon the raw ensemble in terms of CRPS around
30%-33%, in terms of MAE around 22%-25% and in terms of RMSE around 23%-26%,
i.e. they all yield a more pronounced improvement in comparison to the methods
using only the reduced variable set in Section 6.1. This result outlines that
an appropriate selection of predictor variables can enhance model performance.
Furthermore, it should be pointed out that GAM-DVQR-T2 clearly outperforms all
other methods with respect to CRPS, MAE and RMSE followed by GAM-DVQR-T1 for
the coverage and width scores in Table 6.

Figure 10: Kendall’s $\tau$ predictions and its 95% confidence bands for
station Arkona with D-vine GAM copula cutout
$Y_{\mathrm{t2m}}-\overline{X}_{\mathrm{t2m}}-\overline{X}_{\mathrm{tcc}}-\overline{X}_{\mathrm{ws10m}}$
in the training data. The green colour represents the linear correlation model
(T1), the light blue colour the spline-based correlation model (T2) and the
purple color stands for the empirical Kendall’s $\tau$ using the refined
rolling training period with window size $n=25$.

We conclude that GAM-DVQR-T1 and GAM-DVQR-T2 perform better than EMOS-GB as
well as GAM-DVQR-C, as it can take better account of the temporal variation of
the predictor variables and the correlation among them. Based on the results
and on Figure 10 for GAM-DVQR-T2 in comparison to GAM-DVQR-T1, we deduce that
the spline based time-dependent correlation model can identify and describe
the behavior of the time-varying Kendall’s $\tau$ more accurately in a higher-
dimensional variable setting than the linear model for Kendall’s $\tau$.
Furthermore, if the correlations among the predictor variables are correctly
specified, this choice leads to a more reasonable variable selection for GAM-
DVQR, and facilitates more appropriate dependencies as well as interactions
between the variables than EMOS-GB.

Figure 11: Left: Boxplots of station-specific mean CRPS values of the
considered methods in the validation period. Right: Boxplots of the station-
specific CRPSS of the considered methods over EMOS-GB in the validation
period. Outliers $\pm 1.5\cdot\mathrm{IQR}$ are omitted for better visual
representation.

#### CRPS comparisons.

As in Section 6.1 we investigate the CRPS in more detail. We observe in the
boxplots in Figure 11 that GAM-DVQR-T1 and GAM-DVQR-T2 lead to overall lower
CRPS values than the other methods. Furthermore, both models have smaller
variance in the CRPS values. This point becomes clearer with respect to the
CRPSS of the GAM-DVQR methods over EMOS-GB. While GAM-DVQR-C has less
variation in the skill scores over EMOS-GB which might be led back to the
constant correlation parameter, GAM-DVQR-T1 and GAM-DVQR-T2 show clearly more
variance. GAM-DVQR-T2 yields the highest median CRPSS improvement of about
2.7% over EMOS-GB followed by GAM-DVQR-T1 with about 2.4% and GAM-DVQR-C with
about $-0.5\%$. The results with respect to CRPSS also show, that especially
the time-dependent GAM-DVQR models can appropriately capture the time-
dependent correlation between the variables.

Figure 12: Highest CRPSS of the considered methods over EMOS-GB in % in the
validation period. CRPSS $>12\%$ are visualized in black for a better
representation. Numbers in brackets denote the count.

Figure 12 shows the method with the highest CRPSS in colour over EMOS-GB. It
should be pointed out that the CRPSS of the considered methods over EMOS-GB
seems to depend on the elevation of the stations. The higher the stations are
located, i.e. the more we go into the south of Germany, the higher the CRPSS
of the GAM-DVQR methods over EMOS-GB become. At around 89% of the stations,
the GAM-DVQR methods perform better with respect to CRPSS than EMOS-GB.
Moreover, GAM-DVQR-T2 outperforms EMOS-GB at around 49% of all stations in
terms of CRPSS, followed by GAM-DVQR-T1 (34%) and GAM-DVQR-C (17%).

Figure 13: Percentage of pair-wise Diebold-Mariano (DM) tests for the 2
$\mathrm{m}$ surface temperature forecasts indicating statistically
significant CRPS differences after applying a Benjamini-Hochberg procedure to
account for multiple testing for a nominal level of 0.05 of the corresponding
one-sided tests.

#### Testing CRPS differences.

Finally, we have a look at the statistical significance of the differences in
the predictive performance with respect to CRPS between the methods. In Table
13 the $(i,j)$-entry in the $i$-th row and $j$-th column indicates the
percentage of tests where the null hypothesis of equal predictive performance
of the corresponding one-sided DM test is rejected in favor of the model in
the $i$-th row when compared to the model in the $j$-th column. The remainder
of the sum of $(i,j)$\- and $(j,i)$-entry to 100% is the percentage where the
score differences are not significant.

All ensemble postprocessing models significantly outperform the raw ensemble
for more than 90% of all stations with respect to CRPS. Furthermore, GAM-
DVQR-T2 yields significantly lower CRPS values than EMOS-GB for one third of
all stations, followed by GAM-DVQR-T1 (33.12%) and GAM-DVQR-C (1.52%). It
should also be highlighted that EMOS-GB (90.04%) performs at around 8% of the
stations significantly worse in comparison to the raw ensemble than EMOS
(98.70%), while GAM-DVQR-T1/T2 (95.24/95.45%) shows only at 1% of all stations
on the extended variable set significantly lower CRPS values in comparison to
its version on the reduced variable set (96.32%). We conclude that the GAM-
DVQR-T1/T2 models lead to more stable results over the raw ensemble compared
to EMOS and its gradient-boosted extension EMOS-GB, regardless of whether the
reduced or extended variable set is used.

##  7 Conclusion and outlook

D-vine GAM copula quantile regression (GAM-DVQR) is a powerful statistical
method which allows to select important predictor variables, to model
nonlinear relationships between the considered variables and to simultaneously
take account of covariate effects linked to the Kendall’s $\tau$ of a pair-
copula. We complement the presentation of this new method with the R-package
gamvinereg by [27].

In the application for ensemble postprocessing of 2 $\mathrm{m}$ surface
temperature forecasts, GAM-DVQR is able to capture temporal correlation and to
choose predictor variables accordingly. Furthermore, the main reasons for the
overall better performance of GAM-DVQR over the other methods can be traced
back to the modeling of the temporal behavior of the marginal distributions as
well as the correlations among the predictor variables. For the reduced and
extended variable set, the correlation time-dependent GAM-DVQR models
outperform the constant correlation GAM-DVQR models. This indicates the
presence of a non-constant correlation among the variables which needs to be
included into the model. Furthermore, the time varying GAM-DVQR yields
significant improvements over the benchmark methods EMOS and EMOS-GB. Due to
the static training period used for GAM-DVQR in comparison to the conventional
day-by-day sliding training window for DVQR, the estimation procedure is more
economical, can even result in better fits and makes it appealing for
practical and possibly operational use. In future research, we will
investigate an extension of the GAM-DVQR method allowing for spatial and
spatio-temporal effects for Kendall’s $\tau$. This can be specifically
relevant in the field of ensemble postprocessing, and the new approach can be
compared with other spatial or spatio-temporal postprocessing models, such as
e.g. Markovian EMOS by [39]. It might also be beneficial to test other
covariate effects besides of the mentioned ones, and to use GAM-DVQR for the
postprocessing of non-Gaussian weather quantities, such as wind speed or
precipitation.

In terms of the method itself, the extension of GAM-DVQR to very high-
dimensional settings is on the top of our agenda. The work of [50] can serve
as a starting point and we plan to compare our extension with suitable methods
in various fields.

Additionally, the method could be further refined to deal with discrete
variables, where the work of [45] can be considered. Last but not least, the
extension of our method to more general vine structures, e.g. C-vine
(canonical vine) or R-vine (regular vine) based on Tepegjozova2022, Zhu2021,
respectively, would allow for different dependence structures.

## Acknowledgements

We are grateful to the European Centre for Medium-Range Weather Forecasts
(ECMWF) and the German Weather Service (DWD) for providing forecasts and
observation data, respectively. Furthermore, the authors acknowledge support
of the research by Deutsche Forschungsgemeinschaft (DFG) Grant Number MO
3394/1-1, and by the Hungarian National Research, Development and Innovation
Office under Grant Number NN125679. Annette Möller acknowledges support by the
Helmholtz Association’s pilot project “Uncertainty Quantification”.

Appendix

##  Appendix A Hyperparameter specifications

Method | Hyperparameter | Value
---|---|---
EMOS | loss function | CRPS
| number of max. iterations | 5000
| stopping criterion | relative threshold $1\mathrm{e}^{-8}$
EMOS-GB | loss function | LogS
| number of max. iterations | 500
| step size | 0.05
| stopping criterion | AIC
DVQR | window size $n$ | 25
Table 7: Overview of the hyperparameter specifications.

##  Appendix B Marginal distributions

Distribution $\mathcal{D}$-sets:

  * •

$A:=\left\\{\mathcal{N}(\mu,\sigma),\mathcal{S}\mathcal{N}(\mu,\sigma,\nu),\mathcal{S}t(\mu,\sigma,\nu,\tau)\right\\}$,

  * •

$B:=\left\\{\text{logit}\mathcal{N}(\mu,\sigma),\text{logit}\mathcal{S}\mathcal{N}(\mu,\sigma,\nu),\text{logit}\mathcal{S}t(\mu,\sigma,\nu,\tau),\mathcal{B}(\mu,\sigma,\nu,\tau)\right\\}$,

  * •

$C:=\left\\{\log\mathcal{N}(\mu,\sigma),\log\mathcal{S}\mathcal{N}(\mu,\sigma,\nu),\log\mathcal{S}t(\mu,\sigma,\nu,\tau)\right\\}$,

where $\mathcal{N}$ denotes the Gaussian normal distribution, $t$ denotes the
Student-$t$ distribution and $\mathcal{B}$ represents the Beta distribution.
Furthermore $\mathcal{S}$ abbreviates the skewed version of a distribution,
e.g. $\mathcal{SN}$ denotes the skew Gaussian normal distribution and logit as
well as $\log$ denote the transformation of the response with the logit- or
log-function.

Variable | Tested $\mathcal{D}$-set | Selected $\mathcal{D}$ |  Remark
---|---|---|---
$Y_{\mathrm{t2m}}$ | $A$ | $\mathcal{N}(\mu,\sigma)$ |
$\overline{X}_{\mathrm{t2m}}$ | $A$ | $\mathcal{N}(\mu,\sigma)$ |
$\overline{X}_{\mathrm{d2m}}$ | $A$ | $\mathcal{S}t(\mu,\sigma,\nu,\tau)$ |
$\overline{X}_{\mathrm{pr}}$ | $A$ | $\mathcal{N}(\mu,\sigma)$ |
$\overline{X}_{\mathrm{sr}}$ | $A$ | $\mathcal{S}\mathcal{N}(\mu,\sigma,\nu)$ |
$\overline{X}_{\mathrm{u10m}}$ | $A$ | $\mathcal{N}(\mu,\sigma)$ |
$\overline{X}_{\mathrm{v10m}}$ | $A$ | $\mathcal{N}(\mu,\sigma)$ |
$\overline{X}_{\mathrm{r2m}}$ | $B$ | $\text{logit}\mathcal{S}t(\mu,\sigma,\nu,\tau)$ |
$\overline{X}_{\mathrm{tcc}}$ | $B$ | $\mathcal{B}(\mu,\sigma,\nu,\tau)$ |  raw data transformation based on [52]
$\overline{X}_{\mathrm{ws10m}}$ | $C$ | $\log\mathcal{S}t(\mu,\sigma,\nu,\tau)$ |
$\overline{X}_{\mathrm{wg10m}}$ | $C$ | $\log\mathcal{N}(\mu,\sigma)$ |
$S_{\mathrm{t2m}}$ | $C$ | $\log\mathcal{N}(\mu,\sigma)$ |
$S_{\mathrm{d2m}}$ | $C$ | $\log\mathcal{N}(\mu,\sigma)$ |
$S_{\mathrm{pr}}$ | $C$ | $\log\mathcal{S}t(\mu,\sigma,\nu,\tau)$ |
$S_{\mathrm{sr}}$ | $C$ | $\log\mathcal{S}t(\mu,\sigma,\nu,\tau)$ |
$S_{\mathrm{u10m}}$ | $C$ | $\log\mathcal{N}(\mu,\sigma)$ |
$S_{\mathrm{v10m}}$ | $C$ | $\log\mathcal{N}(\mu,\sigma)$ |
$S_{\mathrm{r2m}}$ | $B$ | $\text{logit}\mathcal{N}(\mu,\sigma)$ |
$S_{\mathrm{tcc}}$ | $B$ | $\text{logit}\mathcal{S}t(\mu,\sigma,\nu,\tau)$ |  raw data min-max-transformation & transformation based on [52]
$S_{\mathrm{ws10m}}$ | $C$ | $\log\mathcal{N}(\mu,\sigma)$ |
$S_{\mathrm{wg10m}}$ | $C$ | $\log\mathcal{N}(\mu,\sigma)$ |
Table 8: GAMLSS distribution selection.

## References

  * [1] Elif F. Acar, Radu V. Craiu and Fang Yao  “Dependence Calibration in Conditional Copulas: A Nonparametric Approach”  In _Biometrics_ 67.2 Wiley, 2010, pp. 445–453  DOI: 10.1111/j.1541-0420.2010.01472.x
  * [2] Oleg A. Alduchov and Robert E. Eskridge  “Improved Magnus Form Approximation of Saturation Vapor Pressure”  In _Journal of Applied Meteorology_ 35.4 American Meteorological Society, 1996, pp. 601–609  DOI: 10.1175/1520-0450(1996)035<0601:imfaos>2.0.co;2
  * [3] T. Bedford and R. M. Cooke  “Probability Density Decomposition for Conditionally Dependent Random Variables Modeled by Vines”  In _Annals of Mathematics and Artificial Intelligence_ 32.1/4 Springer ScienceBusiness Media LLC, 2001, pp. 245–268  DOI: 10.1023/a:1016725902970
  * [4] T. Bedford and R. M. Cooke  “Vines—A New Graphical Model for Dependent Random Variables”  In _The Annals of Statistics_ 30.4 Institute of Mathematical Statistics, 2002, pp. 1031–1068  DOI: 10.1214/aos/1031689016
  * [5] Yoav Benjamini and Yosef Hochberg  “Controlling the False Discovery Rate: A Practical and Powerful Approach to Multiple Testing”  In _Journal of the Royal Statistical Society: Series B (Methodological)_ 57.1 Wiley, 1995, pp. 289–300  DOI: 10.1111/j.2517-6161.1995.tb02031.x
  * [6] Claudia Czado  “Analyzing Dependent Data with Vine Copulas”  Springer International Publishing, 2019  DOI: 10.1007/978-3-030-13785-4
  * [7] Claudia Czado and Thomas Nagler  “Vine Copula Based Modeling”  In _Annual Review of Statistics and Its Application_ 9.1 Annual Reviews, 2022, pp. 453–477  DOI: 10.1146/annurev-statistics-040220-101153
  * [8] Markus Dabernig, Georg J. Mayr, Jakob W. Messner and Achim Zeileis  “Spatial ensemble post-processing with standardized anomalies”  In _Quarterly Journal of the Royal Meteorological Society_ 143.703 Wiley, 2017, pp. 909–916  DOI: 10.1002/qj.2975
  * [9] A. P. Dawid  “Present Position and Potential Developments: Some Personal Views: Statistical Theory: The Prequential Approach”  In _Journal of the Royal Statistical Society. Series A (General)_ 147.2 JSTOR, 1984, pp. 278  DOI: 10.2307/2981683
  * [10] Jonathan Demaeyer et al.  “The EUPPBench postprocessing benchmark dataset v1.0”  In _Earth System Science Data_ 15.6 Copernicus GmbH, 2023, pp. 2635–2653  DOI: 10.5194/essd-15-2635-2023
  * [11] Francis X. Diebold and Roberto S. Mariano  “Comparing Predictive Accuracy”  In _Journal of Business & Economic Statistics_ 13.3 Informa UK Limited, 1995, pp. 253–263  DOI: 10.1080/07350015.1995.10524599
  * [12] DWD Climate Data Center (CDC)  “Historische stündliche Stationsmessungen der Lufttemperatur und Luftfeuchte für Deutschland, Version v006”, 2018  DWD Climate Data Center (CDC)  URL: https://opendata.dwd.de/climate_environment/CDC/observations_germany/climate/hourly/air_temperature/historical/BESCHREIBUNG_obsgermany_climate_hourly_tu_historical_de.pdf
  * [13] European Centre for Medium-Range Weather Forecasts (ECMWF)  “Gridded forecast” Creative Commons Attribution 4.0 International (CC BY 4.0), 2021  URL: https://www.ecmwf.int
  * [14] Manuel Gebetsberger, Reto Stauffer, Georg J. Mayr and Achim Zeileis  “Skewed logistic distribution for statistical temperature post-processing in mountainous areas”  In _Advances in Statistical Climatology, Meteorology and Oceanography_ 5.1 Copernicus GmbH, 2019, pp. 87–100  DOI: 10.5194/ascmo-5-87-2019
  * [15] Irene Gijbels, Noel Veraverbeke and Marel Omelka  “Conditional copulas, association measures and their applications”  In _Computational Statistics & Data Analysis_ 55.5 Elsevier BV, 2011, pp. 1919–1932  DOI: 10.1016/j.csda.2010.11.010
  * [16] Tilmann Gneiting  “Making and Evaluating Point Forecasts”  In _Journal of the American Statistical Association_ 106.494 Informa UK Limited, 2011, pp. 746–762  DOI: 10.1198/jasa.2011.r10138
  * [17] Tilmann Gneiting, Fadoua Balabdaoui and Adrian E. Raftery  “Probabilistic forecasts, calibration and sharpness”  In _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ 69.2 Wiley, 2007, pp. 243–268  DOI: 10.1111/j.1467-9868.2007.00587.x
  * [18] Tilmann Gneiting and Adrian E Raftery  “Strictly Proper Scoring Rules, Prediction, and Estimation”  In _Journal of the American Statistical Association_ 102.477 Informa UK Limited, 2007, pp. 359–378  DOI: 10.1198/016214506000001437
  * [19] Tilmann Gneiting et al.  “Assessing probabilistic forecasts of multivariate quantities, with an application to ensemble predictions of surface winds”  In _TEST_ 17.2 Springer ScienceBusiness Media LLC, 2008, pp. 211–235  DOI: 10.1007/s11749-008-0114-x
  * [20] Tilmann Gneiting, Adrian E. Raftery, Anton H. Westveld and Tom Goldman  “Calibrated Probabilistic Forecasting Using Ensemble Model Output Statistics and Minimum CRPS Estimation”  In _Monthly Weather Review_ 133.5 American Meteorological Society, 2005, pp. 1098–1118  DOI: 10.1175/mwr2904.1
  * [21] P.J. Green and Bernard. W. Silverman  “Nonparametric Regression and Generalized Linear Models”  ChapmanHall/CRC, 1993  DOI: 10.1201/b15710
  * [22] Wolfgang Härdle, Hua Liang and Jiti Gao  “Partially Linear Models”  Physica-Verlag HD, 2000  DOI: 10.1007/978-3-642-57700-0
  * [23] T.J. Hastie and R.J. Tibshirani  “Generalized Additive Models”  Taylor & Francis, 1990
  * [24] S. Hemri et al.  “Trends in the predictive performance of raw ensemble weather forecasts”  In _Geophysical Research Letters_ 41.24 American Geophysical Union (AGU), 2014, pp. 9197–9205  DOI: 10.1002/2014gl062472
  * [25] David Jobst  “eppverification”, 2021  DOI: 10.5281/zenodo.5574572
  * [26] David Jobst, Annette Möller and Jürgen Groß  “Data set for the ensemble postprocessing of 2m surface temperature forecasts in Germany for 24 hours lead time”  Zenodo, 2023  DOI: 10.5281/zenodo.8127734
  * [27] David Jobst, Annette Möller and Jürgen Groß  “gamvinereg”, 2023  URL: https://github.com/jobstdavid/gamvinereg
  * [28] David Jobst, Annette Möller and JÃºrgen Groß  “D-vine-copula-based postprocessing of wind speed ensemble forecasts”  In _Quarterly Journal of the Royal Meteorological Society_ 149.755 Wiley, 2023, pp. 2575–2597  DOI: 10.1002/qj.4521
  * [29] Harry Joe  “Families of $m$-variate distributions with given margins and $m(m-1)/2$ bivariate dependence parameters”  In _Institute of Mathematical Statistics Lecture Notes - Monograph Series_ Institute of Mathematical Statistics, 1996, pp. 120–141  DOI: 10.1214/lnms/1215452614
  * [30] Harry Joe and James Jianmeng Xu  “The Estimation Method of Inference Functions for Margins for Multivariate Models”  The University of British Columbia, 1996  DOI: 10.14288/1.0225985
  * [31] Gunky Kim, Mervyn J. Silvapulle and Paramsothy Silvapulle  “Comparison of semiparametric and parametric methods for estimating copulas”  In _Computational Statistics & Data Analysis_ 51.6 Elsevier BV, 2007, pp. 2836–2850  DOI: 10.1016/j.csda.2006.10.009
  * [32] Daniel Kraus and Claudia Czado  “D-vine copula based quantile regression”  In _Computational Statistics & Data Analysis_ 110 Elsevier BV, 2017, pp. 1–18  DOI: 10.1016/j.csda.2016.12.009
  * [33] Moritz N. Lang et al.  “Remember the past: a comparison of time-adaptive training schemes for non-homogeneous regression”  In _Nonlinear Processes in Geophysics_ 27.1 Copernicus GmbH, 2020, pp. 23–34  DOI: 10.5194/npg-27-23-2020
  * [34] Martin Leutbecher and Tim N. Palmer  “Ensemble forecasting”  In _Journal of Computational Physics_ 227.7 Elsevier BV, 2008, pp. 3515–3539  DOI: 10.1016/j.jcp.2007.02.014
  * [35] Ker-Chau Li and Naihua Duan  “Regression Analysis Under Link Violation”  In _The Annals of Statistics_ 17.3 Institute of Mathematical Statistics, 1989  DOI: 10.1214/aos/1176347254
  * [36] James E. Matheson and Robert L. Winkler  “Scoring Rules for Continuous Probability Distributions”  In _Management Science_ 22.10 Institute for Operations Researchthe Management Sciences (INFORMS), 1976, pp. 1087–1096  DOI: 10.1287/mnsc.22.10.1087
  * [37] Jakob W. Messner, Georg J. Mayr and Achim Zeileis  “Heteroscedastic Censored and Truncated Regression with crch”  In _The R Journal_ 8.1, 2016, pp. 173–181  DOI: 10.32614/RJ-2016-012
  * [38] Jakob W. Messner, Georg J. Mayr and Achim Zeileis  “Nonhomogeneous Boosting for Predictor Selection in Ensemble Postprocessing”  In _Monthly Weather Review_ 145.1 American Meteorological Society, 2017, pp. 137–147  DOI: 10.1175/mwr-d-16-0088.1
  * [39] Annette Möller, Thordis L. Thorarinsdottir, Alex Lenkoski and Tilmann Gneiting  “Spatially adaptive, Bayesian estimation for probabilistic temperature forecasts”, 2015  arXiv:1507.05066 [stat.AP]
  * [40] Annette Möller et al.  “Vine copula based post-processing of ensemble forecasts for temperature”, 2018  arXiv:1811.02255 [stat.AP]
  * [41] Thomas Nagler  “vinereg: D-Vine Quantile Regression” R package version 0.8.2, 2022  URL: https://CRAN.R-project.org/package=vinereg
  * [42] Thomas Nagler and Thibault Vatter  “gamCopula: Generalized Additive Models for Bivariate Conditional Dependence Structures and Vine Copulas” R package version 0.0-7, 2020  URL: https://CRAN.R-project.org/package=gamCopula
  * [43] Thomas Nagler et al.  “VineCopula: Statistical Inference of Vine Copulas” R package version 2.4.1, 2020  URL: https://CRAN.R-project.org/package=VineCopula
  * [44] T. N. Palmer  “The economic value of ensemble forecasts as a tool for risk assessment: From days to decades”  In _Quarterly Journal of the Royal Meteorological Society_ 128.581 Wiley, 2002, pp. 747–774  DOI: 10.1256/0035900021643593
  * [45] Anastasios Panagiotelis, Claudia Czado and Harry Joe  “Pair Copula Constructions for Multivariate Discrete Data”  In _Journal of the American Statistical Association_ 107.499 Informa UK Limited, 2012, pp. 1063–1072  DOI: 10.1080/01621459.2012.682850
  * [46] A. J. Patton  “Applications of Copula Theory in Financial Econometrics”, 2002
  * [47] R Core Team  “R: A Language and Environment for Statistical Computing”, 2020  R Foundation for Statistical Computing  URL: https://www.R-project.org/
  * [48] Stephan Rasp and Sebastian Lerch  “Neural Networks for Postprocessing Ensemble Weather Forecasts”  In _Monthly Weather Review_ 146.11 American Meteorological Society, 2018, pp. 3885–3900  DOI: 10.1175/mwr-d-18-0187.1
  * [49] R. A. Rigby and D. M. Stasinopoulos  “Generalized additive models for location, scale and shape”  In _Applied Statistics_ 54, 2005, pp. 507–554
  * [50] Özge Sahin and Claudia Czado  “High-dimensional sparse vine copula regression with application to genomic prediction”, 2022  arXiv:2208.12383 [stat.ME]
  * [51] Abe Sklar  “Fonctions de Répartition à Dimensions et Leurs Marges”  In _Publications de L’Institut de Statistique de L’Université de Paris_ 8, 1959, pp. 229–231
  * [52] Michael Smithson and Jay Verkuilen  “A better lemon squeezer? Maximum-likelihood regression with beta-distributed dependent variables.”  In _Psychological Methods_ 11.1 American Psychological Association (APA), 2006, pp. 54–71  DOI: 10.1037/1082-989x.11.1.54
  * [53] Jakob Stöber, Harry Joe and Claudia Czado  “Simplified pair copula constructions—Limitations and extensions”  In _Journal of Multivariate Analysis_ 119 Elsevier BV, 2013, pp. 101–118  DOI: 10.1016/j.jmva.2013.04.014
  * [54] Maxime Taillardat  “Skewed and Mixture of Gaussian Distributions for Ensemble Postprocessing”  In _Atmosphere_ 12.8 MDPI AG, 2021, pp. 966  DOI: 10.3390/atmos12080966
  * [55] Maxime Taillardat, Olivier Mestre, MichaÃ«l Zamo and Philippe Naveau  “Calibrated Ensemble Forecasts Using Quantile Regression Forests and Ensemble Model Output Statistics”  In _Monthly Weather Review_ 144.6 American Meteorological Society, 2016, pp. 2375–2393  DOI: 10.1175/mwr-d-15-0260.1
  * [56] Olivier Talagrand, Robert Vautard and Bernard Strauss  “Evaluation of probabilistic prediction systems”  In _Proc. Workshop on Predictability_ , 1997, pp. 1–25  URL: https://www.ecmwf.int/en/elibrary/12555-evaluation-probabilistic-prediction-systems
  * [57] Marija Tepegjozova, Jing Zhou, Gerda Claeskens and Claudia Czado  “Nonparametric C- and D-vine-based quantile regression”  In _Dependence Modeling_ 10.1 Walter de Gruyter GmbH, 2022, pp. 1–21  DOI: 10.1515/demo-2022-0100
  * [58] Thibault Vatter and Valérie Chavez-Demoulin  “Generalized additive models for conditional dependence structures”  In _Journal of Multivariate Analysis_ 141 Elsevier BV, 2015, pp. 147–167  DOI: 10.1016/j.jmva.2015.07.003
  * [59] Thibault Vatter and Thomas Nagler  “Generalized Additive Models for Pair-Copula Constructions”  In _Journal of Computational and Graphical Statistics_ 27.4 Informa UK Limited, 2018, pp. 715–727  DOI: 10.1080/10618600.2018.1451338
  * [60] D. S. Wilks  “The Stippling Shows Statistically Significant Grid Points: How Research Results are Routinely Overstated and Overinterpreted, and What to Do about It”  In _Bulletin of the American Meteorological Society_ 97.12 American Meteorological Society, 2016, pp. 2263–2273  DOI: 10.1175/bams-d-15-00267.1
  * [61] Kailun Zhu, Dorota Kurowicka and Gabriela F. Nane  “Simplified R-vine based forward regression”  In _Computational Statistics Data Analysis_ 155 Elsevier BV, 2021, pp. 107091  DOI: 10.1016/j.csda.2020.107091