text:
\section{Introduction}
The European Union (EU) emissions trading system (ETS) allows polluters to trade emission allowances (European emission allowances, EUA). Their supply is fixed, while the demand is determined by the market participants.
Understanding the underlying dynamics of this market is essential for policymakers, commodity trading, and sustainability planning. This motivates exploring the determinants of the EUA price and its statistical properties.

The EU-ETS was set up in 2005, making it the first ETS worldwide. Not much later, research picked up the market. The EU-ETS has significant implications for energy companies and the industrial sector. As a side-effect, the introduction of the EU-ETS has strengthened the link between gas and power prices, leading to undesirable geopolitical uncertainties in the natural gas price. More than  fifteen years later, following the Russian invasion of Ukraine, the geopolitical uncertainties of natural gas supply are more pronounced than ever before. It is, therefore, plausible to assume that some dependence between energy assets and the EUA price exists. Early studies dealing with these relations support this view:
\citet{fezzi2009interaction} employ a cointegrated VAR model and find gas prices to induce quick reactions of EUA prices, while electricity prices show a delayed reaction to shocks in EUA prices.

The dependence of EU-ETS on other energy markets (including coal, oil, natural gas, and electricity) has  been analyzed more recently by \citet{chevallier2019conditional}, who find a negative link between carbon prices and oil and gas prices. \citet{meier2020commodities} also find that EUA prices are significantly influenced by fossil fuel prices, while \cite{zhou2020measurements} show evidence that carbon price itself drives the price of renewable energies.
Employing an artificial neural network,  ~\cite{garcia2020short} pronounce the effect of subjective economic and political decisions on the EUA price and reject the influence of the external variables considered (iron, electricity, and steel prices).
\citet{duan2021marginal} find a negative impact of energy prices on EUA prices. The effect is asymmetric, i.e., it is more pronounced in the lower quantiles of the EUA distribution.

Some studies focus rather on forecasting EUA prices, which is  particularly relevant for portfolio management. \citet{paolella2008econometric} look into forecastability and value-at-risk forecasting in the EU-ETS, employing an AR-GARCH model with generalized skewed student-t innovations. \citet{benz2009modeling} use a Markov-switching framework including GARCH effects to model EUA prices, while \citet{trabelsi2022co2} compare the GARCH approach to Generalized Autoregressive Score models.
Many of these studies find EUA prices to exhibit heavy tails, leverage effects, and asymmetries. Therefore student-t and skewed-student-t distributions are often used and have been shown to lead to superior performance in forecasting \citep{trabelsi2022co2, benz2009modeling, paolella2008econometric}.

Another strand of the literature deals with the volatility process of EUA prices and possible spillovers to other markets. \citet{dutta2019assessing} find structural breaks in the EUA volatility process. They show that the forecasting performance of GARCH models increases (while persistence decreases) when accounting for structural breaks. \citet{hanif2021nonlinear} investigate volatility spillovers between EUA prices and energy indices using copulas. Finally, \citet{reboredo2014volatility} analyzes volatility spillovers between oil and EUA markets. Their results suggest leverage effects and volatility dynamics but do not confirm significant spillovers.

The brief literature review above highlights several key characteristics of EUA prices which have implications for modeling and forecasting. First, quite a few studies point towards a relationship to the prices of related commodities like oil, gas, and coal. Second, a normality assumption is likely to fail due to the excess kurtosis and possible skewness of the data. Third, there is evidence for volatility dynamics and even structural breaks.

In this paper, we focus on probabilistic multi-step ahead forecasting of EUA prices. Building upon the aforementioned findings of earlier studies, we propose a VECM-Copula-GARCH model to jointly model EUA, gas, oil and coal prices as a four-variate system. As a combination of three different components the model is able to account for possible long-run equilibrium relations between carbon and fuel prices, autoregressive effects, conditional heteroskedasticity in the individual series as well as for dependence structures between them in a simultaneous yet flexible way. In contrast to earlier combination models which are estimated in stages, we opt for a one-step estimation approach to improve efficiency.

The contributions of this paper are manifold. First and foremost, we jointly model EUA, oil, coal, and natural gas prices using a probabilistic approach. The obtained multivariate predictive distribution can be used for sampling future trajectories and therefore enables probabilistic multi-step forecasts which exploit the dependencies of EUA on other commodity markets. We evaluate the proposed model in an extensive rolling window forecasting study. Several strictly proper scoring rules are employed to evaluate the multivariate forecasts as well as the marginal forecasts. Moreover, we test for significant differences between the forecasts of competing benchmarks. Finally, we compare our results to two standard benchmarks.

The remainder of this paper is organized as follows. The data are presented in Section~\ref{sec_data}. Section~\ref{sec_model} introduces the proposed VECM-Copula-GARCH model and discusses its one-step estimation by maximum likelihood (ML). The set-up for the forecasting study and forecast evaluation measures are outlined in Section~\ref{sec_evaluation}, while Section~\ref{sec_results} discusses the results. Finally, Section~\ref{sec_conclusion} concludes.

\section{Data description and preliminary analysis}\label{sec_data}

The data, plotted in Figure \ref{Fig_dataplot}, comprises 1957 daily observations of short-term futures prices covering the period from December 30, 2014, until May 12, 2022, on four commodities, traded on the Intercontinental Exchange (ICE): European Allowances (carbon emissions), natural gas (Dutch TTF natural gas futures), coal (Rotterdam coal futures), and brent oil. Oil and coal prices, originally given in USD, have been converted from USD into EUR using the daily exchange rate.

\begin{figure}[!h]
    \centering
    \resizebox{\textwidth}{!}{%
        \fbox{\input{anc/log_ts.tex}}
    }
    \caption{Time series plot of the data in levels. The dashed line marks the Russian invasion of Ukraine.}
    \label{Fig_dataplot}
\end{figure}

Raw data have been transformed by taking natural logarithms following the literature convention  to stabilize their variance. Unit root testing by the augmented Dickey-Fuller (ADF) test points to non-stationarity in the levels and stationarity in the first differences of all four series\footnote{Results are omitted for brevity but are available upon request.}. As a next step, their cointegration properties as a system are examined by Johansen's likelihood ratio (LR) trace test. The results are summarized in Table~\ref{Tcoint}.
Despite the data plot suggesting co-movements between the series and in contrast to previous studies \citep[e.g.,][]{Lahiani2017}, no cointegrating relations emerge based on the formal test results. Possible explanation could be found in the non-linear transformation applied to the data, heteroskedasticity and/or structural breaks in our sample.
Guided by earlier results in the literature, we shall nevertheless allow for the possibility to include long-term equilibrium relationships in our proposed model. This will allow us to evaluate whether they contribute for better forecasting performance.

\begin{table}
    \begin{center}
        \begin{threeparttable}
            \caption{Johansen's LR trace cointegrating rank test, 4-variate system with EUA, oil, coal and gas prices in levels  after logarithmic transformation}\label{Tcoint}
            \begin{tabular}{l|llrrr}
                \toprule
                             &                     & \multicolumn{3}{c}{Critical values}                         \\
                $	H_0	$      & LR trace log-levels & $	10 \%	$                           & $	5\%	$   & $	1\%	$   \\ \midrule
                $	r=0	$      & $31.30 $            & $	49.65	$                           & $	53.12	$ & $	60.16	$ \\
                $	r \leq 1	$ & $15.85 $            & $	32.00	$                           & $	34.91	$ & $	41.07	$ \\
                $	r \leq 2	$ & $5.63 $             & $	17.85	$                           & $	19.96	$ & $	24.60	$ \\
                $	r \leq 3	$ & $0.58 $             & $	7.52	$                            & $	9.24	$  & $	12.97	$ \\
                \bottomrule
            \end{tabular}
            \begin{tablenotes}
                \begin{scriptsize}

                    \item The lag order is 2. Results remain qualitatively the same for lag orders 3, 4 and 5.
                \end{scriptsize}
            \end{tablenotes}
        \end{threeparttable}
    \end{center}
\end{table}

\section{The VECM-Copula-GARCH model}\label{sec_model}
\subsection{Model design}

The preliminary analysis in the previous section reveals several key characteristics of the data which we explicitly account for in our proposed model. First, the variables are non-stationary.  While cointegration is not formally detected in the transformed data, we wish to explore whether allowing for it may nevertheless improve forecasting performance, so we employ a VECM for the mean process, and experiment with different cointegrating ranks, including rank zero, i.e. no cointegration. Second, consistent with the literature, we consider univariate GARCH models for the conditional heteroscedasticity of each series. Third, we model the conditional cross-sectional dependence structure of the individual time series by a copula, thereby accounting for both linear and non-linear dependence in a consistent manner. For more flexibility, the dependence is allowed to vary over time, which is implemented by modeling the copula parameters by GARCH-type processes.
Hence our proposed model can be regarded as a VECM-Copula-GARCH-type model with time-varying dependence parameters. As such it adds sophisticated dependency modeling to an otherwise standard in the financial econometrics multivariate conditional mean and variance model.

Denoting the $K$-dimensional vector of observable variables at time $t$ by $\boldsymbol X_t$, our general model can be formulated in terms of the conditional joint distribution $\boldsymbol F_{\boldsymbol X_t|\mathcal{F}_{t-1}}$, where $\mathcal{F}_{t}$ is the sigma field generated by all information available up to and including time $t$. Using Sklar's Theorem, we specify the joint conditional distribution $\boldsymbol F_{\boldsymbol X_t|\mathcal{F}_{t-1}}$ of the model variables in terms of their marginal distributions
$F_{k,X_{k,t}|\mathcal{F}_{t-1}}$ for $k=1,\ldots, K$, and the copula $\boldsymbol C_{\boldsymbol U_{k,t}|\mathcal{F}_{t}}$, as is popular in many financial applications \citep{jondeau2006copula, hu2006dependence}. Thus, for the realized values $\boldsymbol x_t= (x_{1t},\ldots, x_{Kt})^\intercal$ it holds  that
\begin{align}
    \boldsymbol F_{\boldsymbol X_t|\mathcal{F}_{t-1}}(\boldsymbol x_t)
    = \boldsymbol C_{\boldsymbol U_{k,t}|\mathcal{F}_{t}}(\boldsymbol u_t),
    \label{eq_F_marg_cop_decomp}
\end{align}
where $\boldsymbol u_t=(u_{1t},\ldots, u_{Kt})^\intercal$ and $u_{kt} = F_{k,X_{k,t}|\mathcal{F}_{t-1}}(x_{kt})$, $k=1, \ldots, K$.
That is, we specify $K$ independent marginal model components along with a copula which captures their dependency structure. In order to simplify the notation we drop the conditioning on $\mathcal{F}_{t-1}$ and $\mathcal{F}_{t}$ in equation~\eqref{eq_F_marg_cop_decomp} and write henceforth $\boldsymbol F(\boldsymbol x_t) = \boldsymbol C(\boldsymbol u_t)$.

The model can thus be specified as follows:
\begin{equation}
    \boldsymbol F(\boldsymbol x_t) = \boldsymbol C\left[\mathbf{F}(\boldsymbol x_t;\mbox{\boldmath$\mu$}_t, \mbox{\boldmath$\sigma$}_{t}^2, \mbox{\boldmath$\nu$}, \mbox{\boldmath$\lambda$});\Xi_t, \Theta\right], \label{eq_copula}
\end{equation}
where the $m$-dimensional vector $\Xi_{t}$ comprises the time-varying dependence parameters, while the remaining time-invariant copula parameters are gathered in the vector $\Theta$. Specifically, we take $\boldsymbol C$ to be the $t$-copula to allow for heavy tails, in which case $\Xi_{t}$ stands for the time-varying correlation matrix, while $\Theta$ denotes the degrees of freedom parameter.
The individual marginal distributions, succinctly represented by the vector $\mathbf{F} = (F_1, \ldots, F_K)^{\intercal}$ in Eq. \eqref{eq_copula}, are specified in terms of the $(K\times 1)$ vectors of time-varying marginal means $\mbox{\boldmath$\mu$}_t$=$(\mu_{1,t}, \ldots, \mu_{K,t})^{\intercal}$ and variances $\mbox{\boldmath$\sigma$}^2_t$ $=(\sigma^2_{1,t}, \ldots, \sigma^2_{K,t})^{\intercal}$, respectively. In order to allow for both skewness and heavy tails in the marginals, we consider the generalized non-central $t$-distribution for each $F_i$, $i=1,\ldots, K$, with individual degrees and freedom and non-centrality parameters gathered in \mbox{\boldmath$\nu$} $= (\nu_1, \ldots, \nu_K)^{\intercal}$ and  \mbox{\boldmath$\lambda$} $= (\lambda_1, \ldots, \lambda_K)^{\intercal}$. For estimation each $F_i$ has been parameterized in such a way that its expectation and the variance are given by the parameters to be estimated.

We now turn our attention to the mean and variance processes $\mbox{\boldmath$\mu$}_t$ and $\mbox{\boldmath$\sigma$}^2_t$. The temporal evolution of the means vector is modeled by a VECM model as
\begin{align}
    \Delta \mbox{\boldmath$\mu$}_t := \Delta \boldsymbol x_{t} = & \Pi \boldsymbol x_{t-1} + \Gamma \Delta \boldsymbol x_{t-1} + \epsilon_t, \label{eq_vecm}
\end{align}
where $\Pi = \alpha \beta^{\intercal}$, $\alpha, \beta \in {\mathbb R}^{K\times r}$, is the cointegrating matrix of rank $r$, $0 \leq r\leq K$.
\\
The dynamics of the individual variances are modeled by univariate GARCH models:
\begin{align*}
    \sigma_{i,t}^2 = & \omega_i + \alpha^+_{i} (\epsilon_{i,t-1}^+)^2 + \alpha^-_{i} (\epsilon_{i,t-1}^-)^2 + \beta_i \sigma_{i,t-1}^2,
\end{align*}
where $\epsilon_{i,t-1}^+ = \max\{\epsilon_{i,t-1}, 0\}$ and $\epsilon_{i,t-1}^- = \min\{\epsilon_{i,t-1}, 0\}$, $i=1, \ldots, K$, respectively. Separating the coefficients for positive and negative innovations allows the model to capture leverage effects.

Finally, the temporal evolution of the $\frac{K \times (K-1)}{2}$ primary dependence parameters of the copula, gathered in the (scaled) correlation matrix $\Xi_t$, is modeled similarly to a GARCH process:
\begin{align*}
    \Xi_{t} =    & \Lambda\left(\mbox{\boldmath$\xi$}_{t}\right),
    \\
    \xi_{ij,t} = & \eta_{0,ij} + \eta_{1,ij} \xi_{ij,t-1} + \eta_{2,ij} z_{i,t-1} z_{j,t-1},
\end{align*}
where $\xi_{ij,t}$ is a latent process, $z_{i,t}$ denotes the $i$-th element of the standardized residual from  Eq. \eqref{eq_vecm}, $i=1,\ldots, K$,  and  $\Lambda(\cdot)$ is a link function. It ensures that $\Xi_{t}$ does not exceed its support space and remains semi-positive definite by replacing negative values in the singular value decomposition of $\Xi_{t}$ by $0.00001$.

\subsection{Relevant nested models}

The proposed model \eqref{eq_copula} is the most general model incorporating cointegration, conditional heteroskedasticity and time-varying dependence parameters. This model is denoted by \textbf{$\text{VECM}_{\text{lev}, \text{ncp}}^{r, \sigma_t, \rho_t}$}. Here the subscripts $\textbf{lev}$ and $\textbf{ncp}$ denote joint estimation along with the other model parameters of the non-centrality and leverage parameters, respectively, while the superscripts signify
the assumed cointegrating rank ($r$), the time-varying conditional variance ($\sigma_t^2$) and the time-varying dependence parameters $(\rho_t)$.
The performance of this model will be compared to that of several nested models, following the same notational convention.  \textbf{$\text{RW}$}  will stand for a model with no cointegration, so that the dynamics of the mean process can be thought of as a $K$-dimensional random walk. Assuming constant conditional variance and dependence parameters will be denoted by superscripts $\sigma^2$ and $\rho$, respectively, where the dependence on $t$ is suppressed. Hence, e.g., \textbf{$\text{VECM}^{r2, \sigma, \rho}$} will denote the VECM-Copula-GARCH model of cointegrating rank $r=2$, where the dependence parameters and variances are constant in time, and no leverage or non-centrality parameters are estimated.

\subsection{Estimation}
Despite having three different model components, all parameters can be estimated jointly by maximum likelihood (ML). Using conditional independence, the joint likelihood can be factorized as  \begin{align*}
    L = f_{X_1} \prod_{i=2}^T f_{X_i|\mathcal{F}_{i-1}},
\end{align*}
with the multivariate conditional density at any time point $t$ given by
\begin{align*}
    f_{X_t}(\mathbf{x_t} | \mathcal{F}_{t-1}) = c\left[\mathbf{F}(\mathbf{x}_t;\mbox{\boldmath$\mu$}_t, \mbox{\boldmath$\sigma$}_{i,t}, \mbox{\boldmath$\nu$}, \mbox{\boldmath$\lambda$});\Xi_t, \Theta\right] \cdot \prod_{i=1}^K f_{X_{i,t}}(\mathbf{x}_t;\mbox{\boldmath$\mu$}_t, \mbox{\boldmath$\sigma$}_{i,t}, \mbox{\boldmath$\nu$}, \mbox{\boldmath$\lambda$}).
\end{align*}
Here the copula density $c$ can be derived analytically from the copula as $c(\mathbf{u}) = \partial_{\mathbf{u}} C(\mathbf{u})$. The probability integral transform, $F_i(x_{i,t}|\mathcal{F}_{t-1};\mu_{i,t},\sigma_{i,t}, \nu_i, \lambda_i)$, maps the observations for the $i$-th variable on to $[0,1]$.

\section{Forecasting study}\label{sec_evaluation}

\subsection{Study design}

To perform $h$-step ahead probabilistic forecasting from time point $t$ the following procedure is employed. First the model is fitted and the values of all time varying parameters, $\Xi_t$, $\mbox{\boldmath$\mu$}_t$ and $\mbox{\boldmath$\sigma$}_{i,t}$ are determined using the information contained in $\mathcal{F}_{t}$. Then $n$ samples are drawn from the $K$-dimenional copula, $C\left(\cdot;\Xi_t, \Theta\right)$. The $i$-th marginal sample from the copula, $i=1,\hdots,K$, is quantile-transformed using the corresponding $i$-th quantile function,  $F_i^{-1}(\cdot;\mbox{\boldmath$\mu$}_t, \mbox{\boldmath$\sigma$}_{i,t}, \mbox{\boldmath$\nu$}, \mbox{\boldmath$\lambda$})$. The point forecast is estimated by the mean. Using the point forecast, all time-varying parameters are updated and the procedure is repeated $h-1$ times. Thus the generated samples will have a cross-sectional dependence structure imposed by the copula and marginal distributions according to the univariate densities. Expected values, variances and other (mixed) moments can be approximated with arbitrary accuracy, which is limited only by computational power.

\subsection{Forecast evaluation}
Evaluation of the probabilistic forecast is conducted using different evaluation measures. The $K$-variate forecast will be evaluated by the energy score (ES). The energy score of a multivariate probabilistic forecast with distribution $F$ is given by
\begin{align*}
    \text{ES}_t(F, y) = \mathbb{E}_{F} \left(||\tilde{\mathbf{x}}_t - \mathbf{x}_t||_2\right) - \frac{1}{2} \mathbb{E}_F \left(||\tilde{\mathbf{x}}_t - \tilde{\mathbf{x}}_t'||_2 \right),
\end{align*}
where $x_t$ is the observed $K$-dimensional realization and $\tilde{\mathbf{x}}_t$, respectively $\tilde{\mathbf{x}}_t'$ are independent random vectors distributed according to $F$ \citep{gneiting2007strictly}. A low value of the ES indicates a good probabilistic forecast. In the univariate case the energy score becomes the Continous Ranked Probability Score (CRPS).

\section{Results}\label{sec_results}
We consider a rolling window forecasting study with a window size of 1000 days which correspond to about four years of data. As we have 1956 observations in total and compute 30-steps-ahead forecasts, this leaves 926 potential starting points for the rolling window study. To reduce computational costs we sample uniformly $n=250$ data points from the 926 grid points for further evaluation. We report the multivariate distributional forecast of all considered models by providing $2^{12}= 2048$ samples from the predicted distribution. We evaluate the full spectrum of the aforementioned VECM and RW models in conjunction with an exponential smoothing model, denoted by $\text{ETS}^{\sigma}$, which we include as a benchmark. For the sake of brevity, however, we report the performance only of selected ones, in particular the best-performing models. Precisely, we report the results for the following models: $\text{VECM}^{r4,\sigma,\rho_t}$, $\text{VECM}^{r3,\sigma_t,\rho}$, $\text{VECM}^{r2,\sigma,\rho_t}_{\text{ncp}}$, $\text{VECM}^{r1,\sigma,\rho_t}_{\text{ncp}}$, $\text{VECM}^{r0,\sigma,\rho_t}$, $\text{RW}^{\sigma_t, \rho_t}_{\text{ncp}}$, $\text{RW}^{\sigma_t, \rho}_{\text{ncp}}$, $\text{RW}^{\sigma_t, \rho_t}_{\text{ncp, lev}}$, $\text{RW}^{\sigma_t, \rho}_{\text{ncp, lev}}$, $\text{RW}^{\sigma, \rho}$.

The performance of the selected models with respect to the ES can be found in Table \ref{tab:energy}. Their ES is compared to the ES of a simple random walk model with constant copula and no conditional heteroskedasticity modeling, $\text{RW}^{\sigma, \rho}$.
$\text{ES}_{1-30}^{\text{All}}$ denotes the score of the $(4 \times 30)$-dimensional multivariate ($4$-dimensional) 1-30 days ahead forecast. This general form of the ES allows us to evaluate the effects of modeling both the temporal and the cross sectional dependence.
$\text{ES}_{1-30}^{\text{EUA}}$ denotes the score of the univariate 30-dimensional 1-30 days ahead forecast for the EUA prices (analogous for oil, natural gas and coal prices). In this formulation of the ES, only temporal modeling of the time series is evaluated.
$\text{ES}_1^{\text{All}}$ denotes the score of the 4-dimensional multivariate one day-ahead forecast (analogously $\text{ES}_5^{\text{All}}$ denotes the score of the multivariate 5 day ahead forecast). This form of the ES evaluates the effect of modeling the cross-sectional dependence.

\begin{table}[!h]
    \caption{Improvement in energy score of selected models relative to $\textrm{RW}^{\sigma, \rho}$ in \% (higher = better). Colored according to the test statistic of a DM-Test comparing to $\textrm{RW}^{\sigma, \rho}$  (greener means lower test statistic i.e., better performance compared to $\textrm{RW}^{\sigma, \rho}$). \label{tab:energy}}
    \centering
    \resizebox{\linewidth}{!}{
        \begin{tabular}[t]{l>{}r>{}r>{}r>{}r>{}r>{}r>{}r>{}r}
            \toprule
            Model                                       & $\text{ES}^{\text{All}}_{1-30}$         & $\text{ES}^{\text{EUA}}_{1-30}$         & $\text{ES}^{\text{Oil}}_{1-30}$         & $\text{ES}^{\text{NGas}}_{1-30}$        & $\text{ES}^{\text{Coal}}_{1-30}$        & $\text{ES}^{\text{All}}_{1}$            & $\text{ES}^{\text{All}}_{5}$            & $\text{ES}^{\text{All}}_{30}$           \\
            \midrule
            $\textrm{RW}^{\sigma, \rho}$                & \cellcolor[HTML]{FFFFFF}{1.11}          & \cellcolor[HTML]{FFFFFF}{0.37}          & \cellcolor[HTML]{FFFFFF}{0.39}          & \cellcolor[HTML]{FFFFFF}{0.71}          & \cellcolor[HTML]{FFFFFF}{0.47}          & \cellcolor[HTML]{FFFFFF}{0.04}          & \cellcolor[HTML]{FFFFFF}{0.11}          & \cellcolor[HTML]{FFFFFF}{0.28}          \\
            $\textrm{RW}^{\sigma_t, \rho}_{lev, ncp}$   & \cellcolor[HTML]{B5D462}{\textbf{0.76}} & \cellcolor[HTML]{C3D961}{0.59}          & \cellcolor[HTML]{7AC168}{3.83}          & \cellcolor[HTML]{E6E55C}{\textbf{0.43}} & \cellcolor[HTML]{FFDD4D}{-0.24}         & \cellcolor[HTML]{9ACC65}{2.12}          & \cellcolor[HTML]{97CB65}{1.94}          & \cellcolor[HTML]{AFD363}{0.71}          \\
            $\textrm{RW}^{\sigma_t, \rho_t}_{lev, ncp}$ & \cellcolor[HTML]{BDD761}{0.71}          & \cellcolor[HTML]{A3CF64}{0.88}          & \cellcolor[HTML]{7EC268}{\textbf{4.15}} & \cellcolor[HTML]{F3E95A}{0.20}          & \cellcolor[HTML]{FFDD4D}{-0.26}         & \cellcolor[HTML]{97CB65}{\textbf{2.27}} & \cellcolor[HTML]{96CA66}{\textbf{1.96}} & \cellcolor[HTML]{AED263}{\textbf{0.80}} \\
            $\textrm{RW}^{\sigma_t, \rho}_{ncp}$        & \cellcolor[HTML]{E7E55C}{0.25}          & \cellcolor[HTML]{D9E05E}{0.38}          & \cellcolor[HTML]{7AC068}{3.32}          & \cellcolor[HTML]{FFE553}{-0.20}         & \cellcolor[HTML]{FFDE4D}{-0.23}         & \cellcolor[HTML]{9CCC65}{2.24}          & \cellcolor[HTML]{9FCD65}{1.85}          & \cellcolor[HTML]{FDEC58}{0.02}          \\
            $\textrm{RW}^{\sigma_t, \rho_t}_{ncp}$      & \cellcolor[HTML]{EBE65B}{0.22}          & \cellcolor[HTML]{CFDD5F}{0.47}          & \cellcolor[HTML]{82C368}{3.22}          & \cellcolor[HTML]{FFE452}{-0.24}         & \cellcolor[HTML]{FFE754}{-0.10}         & \cellcolor[HTML]{A0CE65}{2.14}          & \cellcolor[HTML]{A0CE65}{1.86}          & \cellcolor[HTML]{F5EA5A}{0.10}          \\
            $\textrm{VECM}^{r0, \sigma, \rho_t}$        & \cellcolor[HTML]{EAE65B}{0.07}          & \cellcolor[HTML]{D0DD5F}{0.52}          & \cellcolor[HTML]{FFB833}{-0.37}         & \cellcolor[HTML]{E7E55C}{0.16}          & \cellcolor[HTML]{81C368}{\textbf{0.54}} & \cellcolor[HTML]{FFCF43}{-0.63}         & \cellcolor[HTML]{FFC43C}{-0.48}         & \cellcolor[HTML]{DEE25D}{0.11}          \\
            $\textrm{VECM}^{r1, \sigma, \rho_t}_{ncp}$  & \cellcolor[HTML]{F36A49}{-4.20}         & \cellcolor[HTML]{FFB531}{-3.67}         & \cellcolor[HTML]{F9893C}{-2.36}         & \cellcolor[HTML]{FB9037}{-4.92}         & \cellcolor[HTML]{FFEB57}{-0.07}         & \cellcolor[HTML]{FFB430}{-1.28}         & \cellcolor[HTML]{FFB12E}{-1.40}         & \cellcolor[HTML]{F05B4E}{-6.16}         \\
            $\textrm{VECM}^{r2, \sigma, \rho_t}_{ncp}$  & \cellcolor[HTML]{F46D48}{-4.63}         & \cellcolor[HTML]{F2E95A}{0.67}          & \cellcolor[HTML]{FFCA40}{-1.91}         & \cellcolor[HTML]{F77D41}{-7.06}         & \cellcolor[HTML]{FFDB4C}{-0.65}         & \cellcolor[HTML]{FFD547}{-0.61}         & \cellcolor[HTML]{FFBE37}{-1.32}         & \cellcolor[HTML]{F46D48}{-6.38}         \\
            $\textrm{VECM}^{r3, \sigma_t, \rho}$        & \cellcolor[HTML]{F9853D}{-6.35}         & \cellcolor[HTML]{FFD145}{-2.32}         & \cellcolor[HTML]{FFD548}{-1.53}         & \cellcolor[HTML]{FC9634}{-9.51}         & \cellcolor[HTML]{FFDD4D}{-0.91}         & \cellcolor[HTML]{EDE75B}{0.38}          & \cellcolor[HTML]{FFE955}{-0.15}         & \cellcolor[HTML]{F77B42}{-9.95}         \\
            $\textrm{VECM}^{r4, \sigma, \rho_t}$        & \cellcolor[HTML]{F46D48}{-9.61}         & \cellcolor[HTML]{FFBB36}{-4.88}         & \cellcolor[HTML]{FFBE37}{-3.83}         & \cellcolor[HTML]{F67844}{-14.67}        & \cellcolor[HTML]{FFCA40}{-2.80}         & \cellcolor[HTML]{FFC83E}{-1.00}         & \cellcolor[HTML]{FFBA35}{-2.42}         & \cellcolor[HTML]{F05C4D}{-14.60}        \\
            $\textrm{ETS}^{\sigma}$                     & \cellcolor[HTML]{EF564F}{-5.21}         & \cellcolor[HTML]{82C368}{\textbf{2.88}} & \cellcolor[HTML]{FFE956}{-0.04}         & \cellcolor[HTML]{FD9C2F}{-2.88}         & \cellcolor[HTML]{FFE855}{-0.40}         & \cellcolor[HTML]{EE5250}{-3.74}         & \cellcolor[HTML]{EE5250}{-5.50}         & \cellcolor[HTML]{EE5250}{-9.26}         \\
            \bottomrule
            \multicolumn{9}{l}{\rule{0pt}{1em}Coloring w.r.t. test statistic: \colorbox[HTML]{66BA6A}{<-5} \colorbox[HTML]{7CC168}{-4} \colorbox[HTML]{91C866}{-3} \colorbox[HTML]{B0D363}{-2} \colorbox[HTML]{D8E05E}{-1} \colorbox[HTML]{FFED58}{0} \colorbox[HTML]{FFD145}{1} \colorbox[HTML]{FFB531}{2} \colorbox[HTML]{FC9733}{3} \colorbox[HTML]{F67744}{4} \colorbox[HTML]{EE5250}{>5} }         \\
        \end{tabular}}
\end{table}

The results of Table \ref{tab:energy} can be summarized as follows. The random walk models perform systematically better than the VECM models on the data in log-levels. Furthermore, the performance of the VECM models deteriorates with increasing cointegrating rank. This is expected, as there is no cointegration present, hence the models are overparameterized. In this case it suffices to consider the time series as random walks and to model them in first differences. The comparison between the random walk models reveals that accounting for time-varying dependence parameters and conditional heteroskedasticity enhances their performance, with the improvement effect being more pronounced for the latter. Furthermore, including the leverage effects in the GARCH process and the non-centrality parameter in the generalized $t$-distribution for the innovations improve the forecast, indicating that asymmetry is an important feature of the data. Comparing the ES of the RW models with the ES of the $\text{ETS}^{\sigma}$ model reveals that the exponential smoothing outperforms the other models with regards to temporal modeling of the EUA time series. However, when considering the overall ES and the ES evaluating the cross-sectional dependence modeling it becomes evident that appropriate dependence modeling is crucial for multivariate forecasting.

Table \ref{tab:crps} presents the CRPS of the selected models for every univariate time series. H1, H5 and H30 denote the score for the 1-, 5- and 30-day-ahead probabilistic forecast, respectively. The CRPS gives information only about the quality of the univariate modeling, as modeling the cross-sectional dependence has no effect on the univariate density which is obtained by quantile transforming the samples generated from the copula. Even though the samples from the copula have a dependence structure imposed by the copula, the marginal distributions are distributed uniformly. The results in Table \ref{tab:crps} align with those of the ES evaluation: in general, the RW models perform best also in terms of modeling the marginal distributions in the short-term. For natural gas prices, the RW models outperform the rest at any forecast horizon, and accounting for leverage effects in the GARCH-part of the model is particularly beneficial for the long-term (30-days-ahead) forecast. In contrast, for long-term forecasts of oil and coal prices, a RW model with constant conditional variance, but accounting for time-varying correlation in the copula, seems to work best. Focusing on the EUA prices, a VECM model with cointegrating rank 2 emerges as a good choice one only when considering a 5-days-ahead forecast; in the long-term the ETS model delivers the best forecasts.

\begin{table}[!h]
    \caption{Improvement in CRPS of selected models relative to $\textrm{RW}^{\sigma, \rho}$ in \% (higher = better). Colored according to the test statistic of a DM-Test comparing to $\textrm{RW}^{\sigma, \rho}$  (greener means lower test statistic i.e., better performance compared to $\textrm{RW}^{\sigma, \rho}$). \label{tab:crps}}
    \centering
    \resizebox{\linewidth}{!}{
        \begin{tabular}[t]{l>{}r>{}r>{}r>{}r>{}r>{}r>{}r>{}r>{}r>{}r>{}r>{}r}
            \toprule
            \multicolumn{1}{c}{CRPS}                    & \multicolumn{3}{c}{EUA}                 & \multicolumn{3}{c}{NGas}                & \multicolumn{3}{c}{Oil}                 & \multicolumn{3}{c}{Coal}                                                                                \\
            \cmidrule(l{3pt}r{3pt}){1-4} \cmidrule(l{3pt}r{3pt}){5-7} \cmidrule(l{3pt}r{3pt}){8-10} \cmidrule(l{3pt}r{3pt}){11-13}
            Model                                       & H1                                      & H5                                      & H30                                     & H1                                      & H5                                      & H30                                     & H1                                      & H5                                      & H30                                     & H1                                      & H5                                      & H30                                     \\
            \midrule
            $\textrm{RW}^{\sigma, \rho}$                & \cellcolor[HTML]{FFFFFF}{0.01}          & \cellcolor[HTML]{FFFFFF}{0.04}          & \cellcolor[HTML]{FFFFFF}{0.08}          & \cellcolor[HTML]{FFFFFF}{0.01}          & \cellcolor[HTML]{FFFFFF}{0.03}          & \cellcolor[HTML]{FFFFFF}{0.09}          & \cellcolor[HTML]{FFFFFF}{0.03}          & \cellcolor[HTML]{FFFFFF}{0.07}          & \cellcolor[HTML]{FFFFFF}{0.17}          & \cellcolor[HTML]{FFFFFF}{0.01}          & \cellcolor[HTML]{FFFFFF}{0.03}          & \cellcolor[HTML]{FFFFFF}{0.12}          \\
            $\textrm{RW}^{\sigma_t, \rho}_{lev, ncp}$   & \cellcolor[HTML]{DFE25D}{0.67}          & \cellcolor[HTML]{F7EA59}{0.14}          & \cellcolor[HTML]{97CB65}{1.33}          & \cellcolor[HTML]{9CCC65}{3.62}          & \cellcolor[HTML]{7EC268}{5.74}          & \cellcolor[HTML]{7FC268}{4.48}          & \cellcolor[HTML]{AFD363}{3.09}          & \cellcolor[HTML]{B5D462}{\textbf{2.42}} & \cellcolor[HTML]{FFEB57}{-0.04}         & \cellcolor[HTML]{DFE25D}{1.04}          & \cellcolor[HTML]{D9E05E}{\textbf{1.06}} & \cellcolor[HTML]{FFAB2A}{-0.75}         \\
            $\textrm{RW}^{\sigma_t, \rho_t}_{lev, ncp}$ & \cellcolor[HTML]{C7DA60}{1.17}          & \cellcolor[HTML]{E1E35D}{0.58}          & \cellcolor[HTML]{89C667}{1.67}          & \cellcolor[HTML]{9DCD65}{3.65}          & \cellcolor[HTML]{84C467}{5.56}          & \cellcolor[HTML]{75BF69}{\textbf{5.65}} & \cellcolor[HTML]{AAD164}{\textbf{3.34}} & \cellcolor[HTML]{B3D463}{2.42}          & \cellcolor[HTML]{FFE351}{-0.25}         & \cellcolor[HTML]{DEE25D}{\textbf{1.25}} & \cellcolor[HTML]{E1E35D}{0.90}          & \cellcolor[HTML]{FFB631}{-0.74}         \\
            $\textrm{RW}^{\sigma_t, \rho}_{ncp}$        & \cellcolor[HTML]{BAD662}{\textbf{1.36}} & \cellcolor[HTML]{F2E95A}{0.21}          & \cellcolor[HTML]{C9DB60}{0.69}          & \cellcolor[HTML]{96CA66}{\textbf{4.01}} & \cellcolor[HTML]{7AC168}{\textbf{5.75}} & \cellcolor[HTML]{94C966}{2.78}          & \cellcolor[HTML]{B6D562}{3.13}          & \cellcolor[HTML]{C1D861}{2.20}          & \cellcolor[HTML]{FFD245}{-0.75}         & \cellcolor[HTML]{E0E35D}{1.09}          & \cellcolor[HTML]{DCE15D}{0.99}          & \cellcolor[HTML]{FFB732}{-0.64}         \\
            $\textrm{RW}^{\sigma_t, \rho_t}_{ncp}$      & \cellcolor[HTML]{C2D961}{1.19}          & \cellcolor[HTML]{F0E85A}{0.25}          & \cellcolor[HTML]{BED761}{0.82}          & \cellcolor[HTML]{98CB65}{3.88}          & \cellcolor[HTML]{80C368}{5.70}          & \cellcolor[HTML]{93C966}{2.92}          & \cellcolor[HTML]{B4D462}{3.18}          & \cellcolor[HTML]{BED761}{2.27}          & \cellcolor[HTML]{FFD648}{-0.64}         & \cellcolor[HTML]{E5E45C}{0.99}          & \cellcolor[HTML]{E2E35D}{0.87}          & \cellcolor[HTML]{FFC33B}{-0.56}         \\
            $\textrm{VECM}^{r0, \sigma, \rho_t}$        & \cellcolor[HTML]{FFEC58}{-0.03}         & \cellcolor[HTML]{C6DA60}{1.00}          & \cellcolor[HTML]{DDE25D}{0.39}          & \cellcolor[HTML]{FFB934}{-1.53}         & \cellcolor[HTML]{FFB02D}{-0.79}         & \cellcolor[HTML]{FFB22F}{-0.49}         & \cellcolor[HTML]{FFD749}{-0.91}         & \cellcolor[HTML]{FFCF43}{-0.60}         & \cellcolor[HTML]{ECE75B}{\textbf{0.13}} & \cellcolor[HTML]{FFE150}{-0.51}         & \cellcolor[HTML]{FFDE4E}{-0.26}         & \cellcolor[HTML]{7EC268}{\textbf{0.67}} \\
            $\textrm{VECM}^{r1, \sigma, \rho_t}_{ncp}$  & \cellcolor[HTML]{FFE553}{-0.45}         & \cellcolor[HTML]{FFE452}{-0.43}         & \cellcolor[HTML]{FEA02C}{-7.89}         & \cellcolor[HTML]{FFDD4D}{-0.49}         & \cellcolor[HTML]{FFDB4B}{-0.38}         & \cellcolor[HTML]{F1614C}{-4.69}         & \cellcolor[HTML]{FFC139}{-1.71}         & \cellcolor[HTML]{FFBE37}{-1.77}         & \cellcolor[HTML]{FA893B}{-7.81}         & \cellcolor[HTML]{FFD94A}{-1.14}         & \cellcolor[HTML]{FFEC57}{-0.04}         & \cellcolor[HTML]{F3E95A}{0.53}          \\
            $\textrm{VECM}^{r2, \sigma, \rho_t}_{ncp}$  & \cellcolor[HTML]{FFED58}{-0.01}         & \cellcolor[HTML]{D3DE5F}{\textbf{1.74}} & \cellcolor[HTML]{FFEA56}{-0.32}         & \cellcolor[HTML]{FFCE42}{-1.12}         & \cellcolor[HTML]{FFCC41}{-1.32}         & \cellcolor[HTML]{FFD245}{-2.27}         & \cellcolor[HTML]{FFDE4D}{-0.72}         & \cellcolor[HTML]{FFC43C}{-1.87}         & \cellcolor[HTML]{F36C48}{-11.92}        & \cellcolor[HTML]{FFE150}{-0.59}         & \cellcolor[HTML]{FFE754}{-0.20}         & \cellcolor[HTML]{FFDB4B}{-1.00}         \\
            $\textrm{VECM}^{r3, \sigma_t, \rho}$        & \cellcolor[HTML]{FDEC58}{0.07}          & \cellcolor[HTML]{DBE15E}{1.55}          & \cellcolor[HTML]{FFBE37}{-6.70}         & \cellcolor[HTML]{FCEC59}{0.09}          & \cellcolor[HTML]{FFE452}{-0.41}         & \cellcolor[HTML]{FFE251}{-1.00}         & \cellcolor[HTML]{F5EA5A}{0.46}          & \cellcolor[HTML]{FFED58}{0.02}          & \cellcolor[HTML]{F9843E}{-17.96}        & \cellcolor[HTML]{EDE75B}{0.63}          & \cellcolor[HTML]{FFEA56}{-0.15}         & \cellcolor[HTML]{FFE351}{-0.89}         \\
            $\textrm{VECM}^{r4, \sigma, \rho_t}$        & \cellcolor[HTML]{F6EA5A}{0.43}          & \cellcolor[HTML]{FFD548}{-1.57}         & \cellcolor[HTML]{FFB531}{-9.05}         & \cellcolor[HTML]{FFC93F}{-1.49}         & \cellcolor[HTML]{FFD044}{-1.81}         & \cellcolor[HTML]{FFC83E}{-4.37}         & \cellcolor[HTML]{FFD648}{-1.13}         & \cellcolor[HTML]{FFC53C}{-3.26}         & \cellcolor[HTML]{F15F4D}{-25.93}        & \cellcolor[HTML]{FFE552}{-0.41}         & \cellcolor[HTML]{F6EA5A}{0.49}          & \cellcolor[HTML]{FFCC41}{-3.99}         \\
            $\textrm{ETS}^{\sigma}$                     & \cellcolor[HTML]{FAEB59}{0.07}          & \cellcolor[HTML]{D1DE5F}{0.69}          & \cellcolor[HTML]{85C467}{\textbf{4.92}} & \cellcolor[HTML]{FFAB2A}{-2.18}         & \cellcolor[HTML]{FFE14F}{-0.35}         & \cellcolor[HTML]{F1E85A}{0.19}          & \cellcolor[HTML]{91C866}{1.69}          & \cellcolor[HTML]{FFCF43}{-0.70}         & \cellcolor[HTML]{FC9534}{-4.04}         & \cellcolor[HTML]{FFEC57}{-0.06}         & \cellcolor[HTML]{FFDF4E}{-1.02}         & \cellcolor[HTML]{FFD648}{-2.79}         \\
            \bottomrule
            \multicolumn{13}{l}{\rule{0pt}{1em}Coloring w.r.t. test statistic: \colorbox[HTML]{66BA6A}{<-5} \colorbox[HTML]{7CC168}{-4} \colorbox[HTML]{91C866}{-3} \colorbox[HTML]{B0D363}{-2} \colorbox[HTML]{D8E05E}{-1} \colorbox[HTML]{FFED58}{0} \colorbox[HTML]{FFD145}{1} \colorbox[HTML]{FFB531}{2} \colorbox[HTML]{FC9733}{3} \colorbox[HTML]{F67744}{4} \colorbox[HTML]{EE5250}{>5} }                                                                                \\
        \end{tabular}}
\end{table}

Table \ref{tab:rmse} presents the point forecast performance in terms of the RMSE. All RW models perform approximately equally well. There are no significant differences as indicated by the DM statistic. This is not surprising since distributional modeling has no influence on the point forecast. Rather, the differences emerge from random fluctuations in the simulations.

\begin{table}[!h]

    \caption{Improvement in RMSE score of selected models relative to $\textrm{RW}^{\sigma, \rho}$ in \% (higher = better). Colored according to the test statistic of a DM-Test comparing to $\textrm{RW}^{\sigma, \rho}$  (greener means lower test statistic i.e., better performance compared to $\textrm{RW}^{\sigma, \rho}$). \label{tab:rmse}}
    \centering
    \resizebox{\linewidth}{!}{
        \begin{tabular}[t]{l>{}r>{}r>{}r>{}r>{}r>{}r>{}r>{}r>{}r>{}r>{}r>{}r}
            \toprule
            \multicolumn{1}{c}{RMSE}                    & \multicolumn{3}{c}{EUA}                 & \multicolumn{3}{c}{NGas}                & \multicolumn{3}{c}{Oil}                 & \multicolumn{3}{c}{Coal}                                                                                \\
            \cmidrule(l{3pt}r{3pt}){1-4} \cmidrule(l{3pt}r{3pt}){5-7} \cmidrule(l{3pt}r{3pt}){8-10} \cmidrule(l{3pt}r{3pt}){11-13}
            Model                                       & H1                                      & H5                                      & H30                                     & H1                                      & H5                                      & H30                                     & H1                                      & H5                                      & H30                                     & H1                                      & H5                                      & H30                                     \\
            \midrule
            $\textrm{RW}^{\sigma, \rho}$                & \cellcolor[HTML]{FFFFFF}{0.03}          & \cellcolor[HTML]{FFFFFF}{0.07}          & \cellcolor[HTML]{FFFFFF}{0.14}          & \cellcolor[HTML]{FFFFFF}{\textbf{0.03}} & \cellcolor[HTML]{FFFFFF}{0.07}          & \cellcolor[HTML]{FFFFFF}{0.18}          & \cellcolor[HTML]{FFFFFF}{0.05}          & \cellcolor[HTML]{FFFFFF}{\textbf{0.13}} & \cellcolor[HTML]{FFFFFF}{0.31}          & \cellcolor[HTML]{FFFFFF}{0.03}          & \cellcolor[HTML]{FFFFFF}{0.07}          & \cellcolor[HTML]{FFFFFF}{0.22}          \\
            $\textrm{RW}^{\sigma_t, \rho}_{lev, ncp}$   & \cellcolor[HTML]{FFD84A}{-0.17}         & \cellcolor[HTML]{FAEB59}{0.02}          & \cellcolor[HTML]{F3E95A}{0.07}          & \cellcolor[HTML]{FFAB2A}{-0.51}         & \cellcolor[HTML]{FFE351}{-0.07}         & \cellcolor[HTML]{FDEC58}{0.02}          & \cellcolor[HTML]{FFDA4B}{-0.15}         & \cellcolor[HTML]{FFBD37}{-0.32}         & \cellcolor[HTML]{E3E45C}{\textbf{0.13}} & \cellcolor[HTML]{F5EA5A}{0.05}          & \cellcolor[HTML]{FCEC58}{0.02}          & \cellcolor[HTML]{DAE15E}{0.12}          \\
            $\textrm{RW}^{\sigma_t, \rho_t}_{lev, ncp}$ & \cellcolor[HTML]{E0E35D}{0.18}          & \cellcolor[HTML]{E6E55C}{0.12}          & \cellcolor[HTML]{FFEB57}{-0.01}         & \cellcolor[HTML]{FFE150}{-0.14}         & \cellcolor[HTML]{FFED58}{0.00}          & \cellcolor[HTML]{FFE451}{-0.13}         & \cellcolor[HTML]{FFDF4E}{-0.09}         & \cellcolor[HTML]{FFD144}{-0.53}         & \cellcolor[HTML]{FFCF43}{-0.30}         & \cellcolor[HTML]{E1E35D}{0.33}          & \cellcolor[HTML]{FFD748}{-0.12}         & \cellcolor[HTML]{FCEC59}{0.02}          \\
            $\textrm{RW}^{\sigma_t, \rho}_{ncp}$        & \cellcolor[HTML]{D3DE5F}{0.27}          & \cellcolor[HTML]{F7EA59}{0.03}          & \cellcolor[HTML]{E4E45C}{0.15}          & \cellcolor[HTML]{FFD346}{-0.21}         & \cellcolor[HTML]{FFE452}{-0.07}         & \cellcolor[HTML]{D6DF5E}{0.21}          & \cellcolor[HTML]{FFDB4B}{-0.11}         & \cellcolor[HTML]{FFE452}{-0.10}         & \cellcolor[HTML]{FFDB4B}{-0.12}         & \cellcolor[HTML]{FFE754}{-0.05}         & \cellcolor[HTML]{D1DE5F}{0.30}          & \cellcolor[HTML]{E7E55C}{0.12}          \\
            $\textrm{RW}^{\sigma_t, \rho_t}_{ncp}$      & \cellcolor[HTML]{FFE553}{-0.08}         & \cellcolor[HTML]{DCE15D}{0.15}          & \cellcolor[HTML]{E4E45C}{0.15}          & \cellcolor[HTML]{FFD548}{-0.26}         & \cellcolor[HTML]{FFE552}{-0.10}         & \cellcolor[HTML]{FFDC4C}{-0.15}         & \cellcolor[HTML]{FFCD42}{-0.27}         & \cellcolor[HTML]{FFE04F}{-0.15}         & \cellcolor[HTML]{FFE855}{-0.03}         & \cellcolor[HTML]{FFB933}{-0.35}         & \cellcolor[HTML]{EBE65B}{0.15}          & \cellcolor[HTML]{FFE04F}{-0.07}         \\
            $\textrm{VECM}^{r0, \sigma, \rho_t}$        & \cellcolor[HTML]{FFE854}{-0.26}         & \cellcolor[HTML]{CCDC60}{0.69}          & \cellcolor[HTML]{CFDD5F}{0.65}          & \cellcolor[HTML]{FFB833}{-1.46}         & \cellcolor[HTML]{FFB330}{-0.56}         & \cellcolor[HTML]{FFE855}{-0.03}         & \cellcolor[HTML]{FFED58}{0.01}          & \cellcolor[HTML]{FFC139}{-1.03}         & \cellcolor[HTML]{FFCE42}{-0.32}         & \cellcolor[HTML]{FDEC58}{0.05}          & \cellcolor[HTML]{FFD749}{-0.40}         & \cellcolor[HTML]{F8EB59}{0.03}          \\
            $\textrm{VECM}^{r1, \sigma, \rho_t}_{ncp}$  & \cellcolor[HTML]{FFD044}{-1.45}         & \cellcolor[HTML]{FFDD4D}{-0.66}         & \cellcolor[HTML]{FEA02C}{-7.84}         & \cellcolor[HTML]{FFD447}{-0.78}         & \cellcolor[HTML]{FFB02D}{-0.85}         & \cellcolor[HTML]{EE5250}{-3.81}         & \cellcolor[HTML]{FFE553}{-0.34}         & \cellcolor[HTML]{FFC13A}{-1.65}         & \cellcolor[HTML]{FEA329}{-5.48}         & \cellcolor[HTML]{FFDF4E}{-0.64}         & \cellcolor[HTML]{FFE452}{-0.32}         & \cellcolor[HTML]{CEDD5F}{2.59}          \\
            $\textrm{VECM}^{r2, \sigma, \rho_t}_{ncp}$  & \cellcolor[HTML]{FBEC59}{0.17}          & \cellcolor[HTML]{DDE25D}{\textbf{1.25}} & \cellcolor[HTML]{FFE854}{-0.66}         & \cellcolor[HTML]{FFC63D}{-1.33}         & \cellcolor[HTML]{FFA727}{-2.03}         & \cellcolor[HTML]{FFC039}{-2.72}         & \cellcolor[HTML]{EDE75B}{\textbf{0.82}} & \cellcolor[HTML]{FFC33B}{-1.95}         & \cellcolor[HTML]{FA8A3B}{-8.55}         & \cellcolor[HTML]{FFE854}{-0.17}         & \cellcolor[HTML]{F7EA59}{0.16}          & \cellcolor[HTML]{E3E45C}{1.34}          \\
            $\textrm{VECM}^{r3, \sigma_t, \rho}$        & \cellcolor[HTML]{FFE351}{-0.51}         & \cellcolor[HTML]{ECE75B}{0.70}          & \cellcolor[HTML]{FFB934}{-8.22}         & \cellcolor[HTML]{FFBB35}{-1.16}         & \cellcolor[HTML]{FFBE37}{-1.59}         & \cellcolor[HTML]{FFDB4B}{-1.14}         & \cellcolor[HTML]{FFE955}{-0.18}         & \cellcolor[HTML]{FFDB4C}{-1.46}         & \cellcolor[HTML]{FEA329}{-13.24}        & \cellcolor[HTML]{FFED58}{0.01}          & \cellcolor[HTML]{FFDF4E}{-0.58}         & \cellcolor[HTML]{E2E35D}{1.74}          \\
            $\textrm{VECM}^{r4, \sigma, \rho_t}$        & \cellcolor[HTML]{FFE653}{-0.41}         & \cellcolor[HTML]{FFED58}{0.00}          & \cellcolor[HTML]{FFB32F}{-10.01}        & \cellcolor[HTML]{FFBC36}{-1.54}         & \cellcolor[HTML]{FFBE37}{-2.58}         & \cellcolor[HTML]{FFC139}{-3.39}         & \cellcolor[HTML]{F2E95A}{0.63}          & \cellcolor[HTML]{FFD044}{-2.40}         & \cellcolor[HTML]{FB9137}{-16.85}        & \cellcolor[HTML]{FFEC57}{-0.03}         & \cellcolor[HTML]{EAE65B}{\textbf{0.88}} & \cellcolor[HTML]{E1E35D}{\textbf{2.90}} \\
            $\textrm{ETS}^{\sigma}$                     & \cellcolor[HTML]{D7E05E}{\textbf{0.52}} & \cellcolor[HTML]{ACD163}{0.90}          & \cellcolor[HTML]{8BC667}{\textbf{4.70}} & \cellcolor[HTML]{FFD044}{-0.69}         & \cellcolor[HTML]{CEDD5F}{\textbf{0.65}} & \cellcolor[HTML]{D6DF5E}{\textbf{0.48}} & \cellcolor[HTML]{EBE65B}{0.20}          & \cellcolor[HTML]{FFC73E}{-1.15}         & \cellcolor[HTML]{FFBB35}{-2.91}         & \cellcolor[HTML]{E6E55C}{\textbf{0.60}} & \cellcolor[HTML]{FFC33B}{-4.18}         & \cellcolor[HTML]{FEA22A}{-11.68}        \\
            \bottomrule
            \multicolumn{13}{l}{\rule{0pt}{1em}Coloring w.r.t. test statistic: \colorbox[HTML]{66BA6A}{<-5} \colorbox[HTML]{7CC168}{-4} \colorbox[HTML]{91C866}{-3} \colorbox[HTML]{B0D363}{-2} \colorbox[HTML]{D8E05E}{-1} \colorbox[HTML]{FFED58}{0} \colorbox[HTML]{FFD145}{1} \colorbox[HTML]{FFB531}{2} \colorbox[HTML]{FC9733}{3} \colorbox[HTML]{F67744}{4} \colorbox[HTML]{EE5250}{>5} }                                                                                \\
        \end{tabular}}
\end{table}

In summary, the RW model with time-varying dependence structure and leverage-GARCH volatility modeling, $\text{RW}_{ncp, lev}^{\sigma_t, \rho_t}$ performs best in terms of 1-30 day-ahead probabilistic forecasting.
Figure \ref{r99_sigma} displays the temporal evolution of the  volatilities estimated by this model. It can be observed that the coal price series continously has the highest conditional variance, which remains mostly stable throughout the sample, increasing only after the beginning of the pandemic. The volatilities of the EUA, oil and natural gas prices fluctuate considerably more. Prior to 2020 EUA price volatility seems to co-move more closely with the volatility of natural gas prices and less so with that of oil and coal prices, while in 2020 the volatility spillover from the oil market seems to increase. The volatilities of all four commodity prices increase and co-move strongly after the Russian invasion of Ukraine (marked by the dashed line in Figure \ref{r99_sigma}). This result suggests rising uncertainty on the energy market following the invasion of Ukraine. However, similarly high volatility levels have been present in the emissions and energy markets before. This holds most noticable for EUA and oil prices during the spring of 2020 due to the COVID-19 pandemic, and for natural gas and oil prices in autumn 2021 due to various factors like the suspension of the certification procedure for Nordstream 2.

\let\pgfimageWithoutPath\pgfimage
\renewcommand{\pgfimage}[2][]{\pgfimageWithoutPath[#1]{anc/#2}}
\begin{figure}[!h]
    \centering
    \resizebox{\textwidth}{!}{%
        \fbox{\input{anc/r99_lev1_ncp1_csig0_crho0/plot_sigma.tex}}
    }
    \caption{Temporal evolution of the standard deviation as provided by Leverage GARCH modeling in the $\text{RW}_{ncp, lev}^{\sigma_t, \rho_t}$ model. The $y$ axis is in logarithmic scale for more clarity.}
    \label{r99_sigma}
\end{figure}
Finally, we examine the temporal evolution of the copula-implied pairwise correlations\footnote{It should be noted that these are the correlations between the individual innovations in Eq. \eqref{eq_vecm}, and not between the observed series.}, presented in Figure \ref{r99_rho}. The correlations between the prices of natural gas and coal, EUA and oil, and EUA and coal are approximately constant over time. The correlation between oil and coal fluctuates highly around a level of 0.25, with the frequency of large negative values increasing after 2020. By far, the highest on average and most volatile correlation is that between EUA and natural gas prices. Since 2019 both its mean level and volatility seem to increase, suggesting a structural shift in the link between the corresponding markets. However, immediately after the Russian invasion of the Ukraine, the correlation between EUA and gas prices changes dramatically, dropping from a level of about 0.4 to nearly -0.3. Since the initial shock it has been increasing again, but it has not yet achieved its previous level by the end of our sample. The following market mechanism provides a possible explanation: Following the drastic increase in natural gas prices, induced by a shortage due to the sanctions against Russia, less natural gas is purchased, which reduces the need for EUA. The decreased demand for EUA then lowers its price.

\begin{figure}[!h]
    \centering
    \resizebox{\textwidth}{!}{%
        \fbox{\input{anc/r99_lev1_ncp1_csig0_crho0/plot_rho.tex}}
    }
    \caption{Temporal evolution of the linear dependence as provided by time varying dependence parameter modeling in the $\text{RW}_{ncp, lev}^{\sigma_t, \rho_t}$ model.}
    \label{r99_rho}
\end{figure}

\section{Conclusion}\label{sec_conclusion}
In this paper we propose a VECM-Copula-GARCH model to jointly model EUA, oil, coal and natural gas prices with the aim of producing short-term probabilistic forecasts. Allowing for time-varying dependence parameters in the copula and for leverage effects in the model of the conditional heteroscedasticity allows us not only to utilize important data characteristics to improve the forecasts, but also to analyze the temporal evolution of the estimated parameters, providing insight into the linkages between EUA and fossil fuels markets. An extensive rolling-window forecasting study compares the proposed model with relevant nested and benchmark models in terms of energy score, continuous rank probability score and root mean squared error. Overall, the best performing model is a random walk model with both leverage effects in the GARCH model part and time-varying dependence parameters. This is expected, since no cointegration is found in the log-transformed data. As cointegration has been reported in some earlier studies, and our own preliminary results suggest that the original untransformed data may be cointegrated, it would be interesting to investigate whether the VECM-part of our proposed model would contribute to improved forecasts in this case. More sophisticated dependence models, considering e.g. vine copulas, could be explored as well. We leave these issues for future research.

\clearpage


text:
\section{Supplementary material}

\subsection{5D supersymmetric Lagrangian}

The gauge fields consist of one 5D vector $V^M$ ($M=1,\dots 5$), one 5D spinor $(\lambda, \lambda')$ and a real scalar $\sigma$, transforming as the adjoint irreducible representation (irrep) of $\text{E}_6$. They can be expressed in terms of the following $\mathcal{N}=1$ superfields:
\begin{equation}
    W^\alpha_{78} = \{ V^\mu,\ \lambda \}\,, \quad \Phi_{78} = \{  \Sigma,\ \lambda'\}\,,
\end{equation}
where $\Sigma = (\sigma + i V^5)$. Here, $W^\alpha$ is a vector superfield, while $\Phi$ is a chiral one. Together, they form a 5D hypermultiplet, which resembles that of $\mathcal{N}=2$ supersymmetry in 4D. The 5D action can be written as:
\begin{equation}
    S_{5D} = \frac{1}{4kg^2} \ \mbox{Tr}\; \left[ \left( \frac{1}{4} W^\alpha W_\alpha\ \delta^2 (\bar \theta) + \mbox{h.c.} \right) + \left( e^{-2gV} \nabla_y e^{2gV}\right)^2 \right]\,,
\end{equation}
where
\begin{eqnarray}
    W_\alpha &=& -\frac{1}{4} \bar{D}^2 e^{-2gV} D_\alpha e^{2gV}\,, \\
    \nabla_y e^{2gV} &=& \partial_y e^{2gV} - g\bar{\Phi} e^{2gV} - e^{2gV} g\Phi\,.
\end{eqnarray}

A 5D matter multiplet in the irrep $R$ consists of a Dirac spinor $\psi_R$ and two scalar states $\phi_{R1}$ and $\phi_{R2}$. Writing the Dirac spinor in terms of 4D Weyl spinors as
\begin{equation}
    \psi_R = \begin{pmatrix} \chi_R \\ \bar{\xi}_R \end{pmatrix}\,,
\end{equation}
the fields can be organised into two chiral superfields
\begin{equation}
    \Phi_R = \{ \phi_{R1},\ \chi_R\} \,, \qquad  \Phi_R^c = \{ \phi^\dagger_{R2},\ \xi_R\} \,,
\end{equation}
where the second one transforms as conjugate irrep of $\text{E}_6$. In other words, $\Phi$ contains the left-handed spinor, while $\Phi^c$ contains the charge-conjugate right-handed one.
The 5D action, including interactions with the gauge fields, can be written as
\begin{equation}
    S_{5D} = \bar{\Phi}_R e^{2gV} \Phi_R + \bar{\Phi}_R^c e^{-2gV} \Phi_R^c + \left(\Phi_R^c (\nabla_y + M) \Phi \ \delta(\bar{\theta})+ \mbox{h.c.} \right)\,.
\end{equation}
The third term above takes care of the 5D Lorentz invariance by introducing a coupling between the two chiral supermultiplets. Also, the above Lagrangian has a global U(1) symmetry under which $\Phi_R \to e^{i Q_R \theta}\ \Phi_R$ and $\Phi_R^c \to e^{-i Q_R \theta} \Phi_R^c$\,.

\subsection{$\text{E}_6$ breaking and orbifold boundary conditions}

In general, the components of an $\text{E}_6$ multiplet will receive different parity assignments. Let us consider a parity $P$ that breaks $\text{E}_6 \to \mathcal{H}$, so that an irrep $R$ of $\text{E}_6$ decomposes in components $\chi_j$. The parity assignments of each component $\chi_j$ will read:
\begin{equation}
    P (\chi_j) = \eta_R \ \mathcal{P}_j\,,
\end{equation}
where $\eta_R$ is an overall sign that can be chosen different for different fields in the irrep $R$, while the relative parities $\mathcal{P}_j$ only depend on the $\text{E}_6$ irrep (and are the same for all fields).

We consider the following breaking patterns:
\begin{eqnarray}
\text{A}:\; \text{E}_6 &\to& \text{SO}(10) \times \text{U}(1)_\psi\,, \\
\text{B}:\; \text{E}_6 &\to& \text{SU}(6)_\text{L} \times \text{SU}(2)_\text{R}\,, \\
\text{C}:\; \text{E}_6 &\to& \text{SU}(6)_\text{R} \times \text{SU}(2)_\text{L}\,.
\end{eqnarray}
The intersection of any pair of them leaves the Pati-Salam group unbroken, $\text{SU}(4)\times \text{SU}(2)_\text{L} \times \text{SU}(2)_\text{R}$, plus $\text{U}(1)_\psi$.

For the fundamental, ${\bf 27}$, we have the following decompositions:
\begin{eqnarray}
\begin{split}
\text{E}_6\rightarrow \text{SO}(10)_{Q_\psi} &\rightarrow (\text{SU}(4)\times \text{SU}(2)_\text{L} \times \text{SU}(2)_\text{R})_{Q_\psi} \\
(16)_{1} \,&=\, \left(4,2,1\right)_1\, + \,\left(\bar{4},1,2\right)_1 \\
(10)_{-2} \,&=\, \left(6,1,1\right)_{-2}\, + \,\left(1,2,2\right)_{-2}\\
(1)_{4}\,&=\, \left(1,1,1\right)_4\,;
\end{split}
\end{eqnarray}
\begin{eqnarray}
\begin{split}
\text{E}_6\rightarrow \text{SU}(6)_\text{L}\times\text{SU}(2)_\text{R} &\rightarrow (\text{SU}(4)\times \text{SU}(2)_\text{L} \times \text{SU}(2)_\text{R})_{Q_\psi} \\
\left(\bar{6},2\right) \,&=\, \left(\bar{4},1,2\right)_1\, + \,\left(1,2,2\right)_{-2} \\
\left(15,1\right) \,&=\, \left(4,2,1\right)_1\, + \,\left(6,1,1\right)_{-2}\, +\,\left(1,1,1\right)_4\,;
\end{split}
\end{eqnarray}
\begin{eqnarray}
\begin{split}
\text{E}_6\rightarrow \text{SU}(6)_\text{R}\times\text{SU}(2)_\text{L} &\rightarrow (\text{SU}(4)\times \text{SU}(2)_\text{L} \times \text{SU}(2)_\text{R})_{Q_\psi}\\
\left(6,2\right) \,&=\, \left(4,2,1\right)_1\, + \,\left(1,2,2\right)_{-2} \\
\left(\bar{15},1\right) \,&=\, \left(\bar{4},1,2\right)_1\, + \,\left(6,1,1\right)_{-2}\, +\,\left(1,1,1\right)_4\,.
\end{split}
\end{eqnarray}
The internal parities $\mathcal{P}_j$ are listed in Table~\ref{tab:27}.

\begin{table}[htb]
\centering
\begin{tabular}{l|c|c|c|}
\phantom{\big(}$\bf 27$ & $\text{SO}(10)\times\text{U}(1)_\psi$ & $\text{SU}(6)_\text{L}\times \text{SU}(2)_\text{R}$ & $\text{SU}(6)_\text{R}\times \text{SU}(2)_\text{L}$ \\ \hline\hline
\phantom{\big(}$(4,2,1)_1$ & \textcolor{blue}{even} & \textcolor{blue}{even} & \textcolor{red}{odd} \\ \hline
\phantom{\big(}$(\bar{4},1,2)_1$ & \textcolor{blue}{even} & \textcolor{red}{odd} & \textcolor{blue}{even} \\ \hline
\phantom{\big(}$(6,1,1)_{-2}$ & \textcolor{red}{odd} & \textcolor{blue}{even} & \textcolor{blue}{even} \\ \hline
\phantom{\big(}$(1,2,2)_{-2}$ & \textcolor{red}{odd} & \textcolor{red}{odd} & \textcolor{red}{odd} \\ \hline
\phantom{\big(}$(1,1,1)_4$ & \textcolor{green}{odd} & \textcolor{blue}{even} & \textcolor{blue}{even} \\ \hline
\end{tabular}
\caption{\label{tab:27} Intrinsic parities for the fundamental $\bf 27$ of $\text{E}_6$. The colours label PS fields in the same irrep of $\mathcal{H}$.}
\end{table}

For the adjoint, $\bf 78$, we have the following decompositions:
\begin{eqnarray}
\begin{split}
\text{E}_6\rightarrow \text{SO}(10)_{Q_\psi} &\rightarrow (\text{SU}(4)\times \text{SU}(2)_\text{L} \times \text{SU}(2)_\text{R})_{Q_\psi} \\
(45)_{0} \,&=\, \left(15,1,1\right)_0\,+\,\left(1,3,1\right)_0\,+\,\left(1,1,3\right)_0\,+\,\left(6,2,2\right)_0 \\
(16)_{-3} \,&=\, \left(4,2,1\right)_{-3}\, + \,\left(\bar{4},1,2\right)_{-3} \\
(\bar{16})_{3} \,&=\, \left(\bar{4},2,1\right)_3\, + \,\left(4,1,2\right)_3 \\
(1)_{0}\,&=\, \left(1,1,1\right)_0\,;
\end{split}
\end{eqnarray}
\begin{eqnarray}
\begin{split}
\text{E}_6\rightarrow \text{SU}(6)_\text{L}\times\text{SU}(2)_\text{R} &\rightarrow (\text{SU}(4)\times \text{SU}(2)_\text{L} \times \text{SU}(2)_\text{R})_{Q_\psi} \\
\left(35,1\right) \,&=\, \left(15,1,1\right)_0\,+\, \left(1,3,1\right)_0\,+\,\left(4,2,1\right)_{-3}\,+\,\left(\bar{4},2,1\right)_{3}\,+\, \left(1,1,1\right)_0 \\
\left(1,3\right) \,&=\, \left(1,1,3\right)_0 \\
\left(20,2\right) \,&=\, \left(6,2,2\right)_0\, + \,\left(\bar{4},1,2\right)_{-3}\,+\,\left(4,1,2\right)_{3}\,;
\end{split}
\end{eqnarray}
\begin{eqnarray}
\begin{split}
\text{E}_6\rightarrow \text{SU}(6)_\text{R}\times\text{SU}(2)_\text{L} &\rightarrow (\text{SU}(4)\times \text{SU}(2)_\text{L} \times \text{SU}(2)_\text{R})_{Q_\psi} \\
\left(35,1\right) \,&=\, \left(15,1,1\right)_0\,+\, \left(1,1,3\right)_0\,+\,\left(4,1,2\right)_3\,+\,\left(\bar{4},1,2\right)_{-3}\,+\, \left(1,1,1\right)_0 \\
\left(3,1\right) \,&=\, \left(1,3,1\right)_0 \\
\left(20,2\right) \,&=\, \left(6,2,2\right)_0\, + \,\left(\bar{4},2,1\right)_3\,+\,\left(4,2,1\right)_{-3}\,.
\end{split}
\end{eqnarray}

The internal parities $\mathcal{P}_j$ are listed in Table~\ref{tab:adj}. Note that for the gauge fields, the overall signs are fixed to be $\eta=+1$ for $W^\alpha$ and $\eta=-1$ for $\Phi$ for all parities. Hence, for the choice of parities B--C, it is not possible to obtain zero modes in the fourplets of $\text{SU}(4)$ in $\phi$, which could play the role of gaugino matter fields in the SM.
Instead, for the choice A--C, the gauginos contain the left-handed $\text{SU}(2)_\text{L}$ doublets, hence they cannot be used for the breaking of the PS gauge symmetry via the Scherk-Schwarz mechanism.

\begin{table}[htb]
\centering
\begin{tabular}{l|c|c|c|}
\phantom{\big(}$\bf 78$ & $\text{SO}(10)\times\text{U}(1)_\psi$ & $\text{SU}(6)_\text{L}\times \text{SU}(2)_\text{R}$ & $\text{SU}(6)_\text{R}\times \text{SU}(2)_\text{L}$ \\ \hline\hline
\phantom{\big(}$(15,1,1)_0$ & \textcolor{blue}{even} & \textcolor{blue}{even} & \textcolor{blue}{even}\\ \hline
\phantom{\big(}$(1,3,1)_0$ & \textcolor{blue}{even} & \textcolor{blue}{even} & \textcolor{green}{even}\\ \hline
\phantom{\big(}$(1,1,3)_0$ & \textcolor{blue}{even} & \textcolor{green}{even} & \textcolor{blue}{even}\\ \hline
\phantom{\big(}$(6,2,2)_0$ & \textcolor{blue}{even} & \textcolor{red}{odd} & \textcolor{red}{odd}\\ \hline
\phantom{\big(}$(4,2,1)_{-3}$ & \textcolor{red}{odd} & \textcolor{blue}{even} & \textcolor{red}{odd}\\ \hline
\phantom{\big(}$(\bar{4},1,2)_{-3}$ & \textcolor{red}{odd} & \textcolor{red}{odd} & \textcolor{blue}{even}\\ \hline
\phantom{\big(}$(\bar{4},2,1)_3$ & \textcolor{green}{odd} & \textcolor{blue}{even} & \textcolor{red}{odd}\\ \hline
\phantom{\big(}$(4,1,2)_3$ & \textcolor{green}{odd} & \textcolor{red}{odd} & \textcolor{blue}{even}\\ \hline
\phantom{\big(}$(1,1,1)_0$ & even & \textcolor{blue}{even} & \textcolor{blue}{even}\\ \hline
\end{tabular}
\caption{\label{tab:adj} Intrinsic parities for the adjoint $\bf 78$ of $\text{E}_6$. The colours label PS fields in the same irrep of $\mathcal{H}$.}
\end{table}

A detailed list of the bulk field components with their assigned orbifold parities can be found in Tables~\ref{tab:5Dmatter}, \ref{tab:5DmatterX} and \ref{tab:5Dgauge}, for the choice A--B, which is presented in the main text.

\subsection{Baryon number is conserved}

At the level of the SM gauge symmetry, the model has 5 U(1) factors: $\text{U}(1)_\text{B-L}$ in $\text{SU}(4)$, $\text{U}(1)_\text{R}$ in $\text{SU}(2)_\text{R}$, $\text{U}(1)_\psi$, $\text{U}(1)_{27}$ and $\text{U}(1)_{27'}$. The two last ones are global, while the first 3 are gauged.

Two combinations are broken by the vacuum expectation values that break PS$\times U(1)_\psi$ down to the SM:
\begin{eqnarray}
    \langle ({\bf 1},{\bf 1},{\bf 1})_{-4} \rangle\;\; \mbox{in}\;\; {\bf 27'} &\Rightarrow& -4 Q_\psi - Q_{27'} = 0 \\
     \langle ({\bf 4},{\bf 1},{\bf 2})_{3} \rangle\;\; \mbox{in}\;\; {\bf 78} &\Rightarrow& Q_\text{R}-Q_\text{B-L}+3 Q_\psi= 0
\end{eqnarray}
where we assume that the PS breaking is due to the neutrino singlets zero modes contained in the adjoint (via the Scherk-Scharz mechanism).
We can further define the hypercharge as follows:
\begin{equation}
    Q_\text{Y} = \frac{1}{2} (Q_\text{B-L} + Q_\text{R})\,.
\end{equation}
Henceforth, it remains the freedom of defining two unbroken global symmetries. There is a unique charge which vanishes on both Higgs doublets, hence it will not be broken by the Higgs vacuum expectation values:
\begin{equation}
Q_\text{B} = \frac{1}{4} Q_\text{B-L} + \frac{1}{6} Q_{27} + \frac{1}{12} Q_\psi - \frac{1}{3} Q_{27'} \,,
\end{equation}
where the normalisation matches the baryon number of quarks.
The second combination must be non-zero on at least one of the Higgs doublets: choosing to normalise the charge to be 1 on $\varphi_{h2}$, we obtain
\begin{equation}
Q_\text{K} = \frac{1}{2} Q_\text{R}  + \frac{1}{6} Q_{27} - \frac{1}{6} Q_\psi + \frac{2}{3} Q_{27'} \,.
\end{equation}
The charge assignments are listed in Tables~\ref{tab:5Dmatter}, \ref{tab:5DmatterX} and~\ref{tab:5Dgauge}.

Note that the PS symmetry may also be broken by the $(\bar{\bf 4},{\bf 1},{\bf 2})_1$ in ${\bf 27'}$, however this leads to the same charge definition as above.

\subsection{Localised superpotential on the SO(10) boundary}

In the model with two localised generations on the SO$(10)$ boundary, mixings with the bulk third generation are not possible, as they are forbidden by the $Q_\psi$ charges. As a solution, one can introduce a minimal set of localised fields, such as a set of localised Higgs states in the $10_0$, which allows for
\begin{equation}
    \tilde{y}^i\ \Phi^i_{16_{-1}} \left. \Phi_{27}\right|_{16_1} \Phi_{10_0}\,,
\end{equation}
hence generating a mixing in the left-handed sector of the theory. Note that no additional gauge anomalies are generated by such localised Higgs states. A mixing between the localised and bulk Higgs fields could be generated by SUSY-breaking quartic couplings, or by introducing additional singlets charged under U$(1)_\psi$ yielding the following superpotential:
\begin{equation}
    a\ \Phi_{1_2} \Phi_{1_{-2}} + b\ \Phi_{10_0} \left. \Phi_{27}^c \right|_{10_2} \Phi_{1_{-2}} + c\ \Phi_{1_{2}} \Phi_{1_2} \left. \Phi_{27'}^c \right|_{1_{-4}}\,.
\end{equation}

\newpage

\begin{table}[htb]
\centering\begin{tabular}{|c|c|c|c|c|c|c||c|c|c|c|}
\hline
 $\text{E}_6$ fields &  PS fields & $(\mathbb{Z}_2,\mathbb{Z}_2')$ & SM & name &$Q_\text{B}$&$Q_\text{K}$&$Q_\text{R}$&$Q_\text{B-L}$&$Q_\psi$&$Q_{27}$ \\
\hline
$\Phi_{27}$ & $(\bf{4}, \bf{2}, \bf{1})_{1}$ & $(+,+)$ & $(3,2)_{1/6}$ & $q_L$ & $1/3$ & $0$ &$0$&$1/3$&$1$&$1$ \\
& & & $(1,2)_{-1/2}$ & $l_L$ & $0$ & $0$ &$0$&$-1$&$1$&$1$ \\ \cline{2-7}
& $(\bar{\bf 4}, {\bf 1}, {\bf 2})_{1}$ & $(+,-)$ & $(\bar{3},1)_{1/3}$ & & $1/6$ & $1/2$ &$1$&$-1/3$&$1$&$1$ \\
&&& $(\bar{3},1)_{-2/3}$ & & $1/6$ & $-1/2$ &$-1$&$-1/3$&$1$&$1$  \\
&&& $(1,1)_{1}$ & & $1/2$ & $1/2$ &$1$&$1$&$1$&$1$ \\
&&& $(1,1)_{0}$ & & $1/2$ & $-1/2$  &$-1$&$1$&$1$&$1$ \\\cline{2-7}
& $(\bf{6}, \bf{1}, \bf{1} )_{-2}$ & $(-,+)$ & $(\bar{3},1)_{1/3}$ & & $1/6$ & $1/2$ &$0$&$2/3$&$-2$&$1$ \\
&&& $(3,1)_{-1/3}$ & & $-1/6$ & $1/2$ &$0$&$-2/3$&$-2$&$1$ \\ \cline{2-7}
& $(\bf{1}, \bf{2}, \bf{2} )_{-2}$ & $(-,-)$ & $(1,2)_{1/2}$ & & $0$ & $1$ &$1$&$0$&$-2$&$1$ \\
&&& $(1,2)_{-1/2}$ & & $0$ & $0$ &$-1$&$0$&$-2$&$1$ \\ \cline{2-7}
& $(\bf{1},\bf{1},\bf{1})_{4}$ & $(-,+)$ & $(1,1)_0$ & & $1/2$ & $-1/2$ &$0$&$0$&$4$&$1$ \\
\hline
$\Phi^c_{27}$ & $(\bar{\bf{4}}, \bf{2}, \bf{1})_{-1}$ & $(-,-)$ & $(\bar{3},1)_{-1/6}$ &  & $-1/3$ & $0$ &$0$&$-1/3$&$-1$&$-1$ \\
& & & $(1,2)_{1/2}$ &  & $0$ & $0$ &$0$&$1$&$-1$&$-1$ \\ \cline{2-7}
& $(\bf 4, {\bf 1}, {\bf 2})_{-1}$ & $(-,+)$ & $(3,1)_{-1/3}$ & & $-1/6$ & $-1/2$ &$-1$&$1/3$&$-1$&$-1$ \\
&&& $(3,1)_{2/3}$ & & $-1/6$ & $1/2$ &$1$&$-1/3$&$-1$&$-1$ \\
&&& $(1,1)_{-1}$ & & $-1/2$ & $-1/2$ &$-1$&$-1$&$-1$&$-1$ \\
&&& $(1,1)_{0}$ & & $-1/2$ & $1/2$ &$1$&$-1$&$-1$&$-1$ \\\cline{2-7}
& $(\bf{6}, \bf{1}, \bf{1} )_{2}$ & $(+,-)$ & $(\bar{3},1)_{1/3}$ & & $1/6$ & $-1/2$ &$0$&$2/3$&$2$&$-1$ \\
&&& $(3,1)_{-1/3}$ & & $-1/6$ & $-1/2$ &$0$&$-2/3$&$2$&$-1$ \\ \cline{2-7}
& $(\bf{1}, \bf{2}, \bf{2} )_{2}$ & $(+,+)$ & $(1,2)_{1/2}$ & $\varphi_{h1}$ & $0$ & $0$ &$1$&$0$&$2$&$-1$ \\
&&& $(1,2)_{-1/2}$ & $\varphi_{h2}^\ast$ & $0$ & $-1$ &$-1$&$0$&$2$&$-2$ \\ \cline{2-7}
& $(\bf{1},\bf{1},\bf{1})_{-4}$ & $(+,-)$ & $(1,1)_0$ & & $-1/2$ & $1/2$ &$0$&$0$&$-4$&$-1$ \\
\hline
\end{tabular}
\caption{\label{tab:5Dmatter} Field decomposition for the Matter multiplet $R={\bf 27}$.}
\end{table}

\begin{table}[htb]
\centering\begin{tabular}{|c|c|c|c|c|c|c||c|c|c|c|}
\hline
 $\text{E}_6$ fields &  PS fields & $(\mathbb{Z}_2,\mathbb{Z}_2')$ & SM & name & $Q_\text{B}$&$Q_\text{K}$ & $Q_\text{R}$ & $Q_\text{B-L}$ & $Q_{\psi}$ & $Q_{27'}$ \\
\hline
$\Phi_{27'}$ & $(\bf{4}, \bf{2}, \bf{1})_{1}$ & $(+,-)$ & $(3,2)_{1/6}$ &  & $-1/6$ & $1/2$ &$0$& $1/3$&$1$&$1$\\
& & & $(1,2)_{-1/2}$ &  & $-1/2$ & $-1/2$ &$0$&$-1$&$1$&$1$\\ \cline{2-7}
& $(\bar{\bf 4}, {\bf 1}, {\bf 2})_{1}$ & $(+,+)$ & $(\bar{3},1)_{1/3}$ & ${b'}_R^c$ & $-1/3$ & $1$ &$1$&$-1/3$&$1$&$1$ \\
&&& $(\bar{3},1)_{-2/3}$ & ${t'}_R^c$ & $-1/3$
& $0$ &$-1$&$-1/3$&$1$&$1$  \\
&&& $(1,1)_{1}$ & ${\tau'}_R^c$ & $0$ & $1$ &$1$&$1$&$1$&$1$ \\
&&& $(1,1)_{0}$ & ${\nu'}_R^c$ & $0$ & $0$  &$-1$&$1$&$1$&$1$ \\\cline{2-7}
& $(\bf{6}, \bf{1}, \bf{1} )_{-2}$ & $(-,-)$ & $(\bar{3},1)_{1/3}$ & & $-1/3$ & $1$ &$0$&$2/3$&$-2$&$1$ \\
&&& $(3,1)_{-1/3}$ & & $1/3$ & $1$ &$0$&$-2/3$&$-2$&$1$ \\ \cline{2-7}
& $(\bf{1}, \bf{2}, \bf{2} )_{-2}$ & $(-,+)$ & $(1,2)_{1/2}$ & & $-1/2$ & $3/2$ &$1$&$0$&$-2$&$1$ \\
&&& $(1,2)_{-1/2}$ & & $-1/2$ & $1/2$ &$-1$&$0$&$-2$&$1$ \\ \cline{2-7}
& $(\bf{1},\bf{1},\bf{1})_{4}$ & $(-,-)$ & $(1,1)_0$ & & $0$ & $0$ &$0$&$0$&$4$&$1$ \\
\hline
$\Phi^c_{27'}$ & $(\bar{\bf{4}}, \bf{2}, \bf{1})_{-1}$ & $(-,+)$ & $(\bar{3},2)_{-1/6}$ &  & $1/6$ & $-1/2$ &$0$&$-1/3$&$-1$&$-1$ \\
& & & $(1,2)_{1/2}$ &  & $1/2$ & $1/2$ &$0$&$1$&$-1$&$-1$ \\ \cline{2-7}
& $(\bf 4, {\bf 1}, {\bf 2})_{-1}$ & $(-,-)$ & $(3,1)_{-1/3}$ & & $1/3$ & $-1$ &$-1$&$1/3$&$-1$&$-1$ \\
&&& $(3,1)_{2/3}$ & & $1/3$ & $0$ &$1$&$1/3$&$-1$&$-1$ \\
&&& $(1,1)_{-1}$ & & $0$ & $-1$ &$-1$&$-1$&$-1$&$-1$ \\
&&& $(1,1)_{0}$ & & $0$ & $0$ &$1$&$-1$&$-1$&$-1$ \\\cline{2-7}
& $(\bf{6}, \bf{1}, \bf{1} )_{2}$ & $(+,+)$ & $(\bar{3},1)_{1/3}$ & $B_R^c$ & $-1/3$ & $-1$ &$0$&$2/3$&$2$&$-1$ \\
&&& $(3,1)_{-1/3}$ & $B_R$ & $1/3$ & $-1$ &$0$&$-2/3$&$2$&$-1$ \\ \cline{2-7}
& $(\bf{1}, \bf{2}, \bf{2} )_{2}$ & $(+,-)$ & $(1,2)_{1/2}$ & & $1/2$ & $-1/2$ &$1$&$0$&$2$&$-1$ \\
&&& $(1,2)_{-1/2}$ & & $1/2$ & $-3/2$ &$-1$&$0$&$2$&$-1$ \\ \cline{2-7}
& $(\bf{1},\bf{1},\bf{1})_{-4}$ & $(+,+)$ & $(1,1)_0$ & $\varphi_\psi$ & $0$ & $0$ &$0$&$0$&$-4$&$-1$ \\
\hline
\end{tabular}
\caption{\label{tab:5DmatterX} Field decomposition for the Matter multiplet $R={\bf 27'}$.}
\end{table}

\begin{table}[htb]
\centering\begin{tabular}{|c|c|c|c|c|c|c||c|c|c|}
\hline
 $\text{E}_6$ fields &  PS fields & $(\mathbb{Z}_2,\mathbb{Z}_2')$ & SM & name & $Q_\text{B}$&$Q_\text{K}$ & $Q_\text{R}$ & $Q_\text{B-L}$ & $Q_{\psi}$ \\
\hline
$W^\alpha_{78}$/$\Phi_{78}$ & $(\bf{15}, \bf{1}, \bf{1})_{0}$ & $(+,+)$ & $(8,1)_{0}$ & $G_\mu$ & $0$ & $0$ &$0$& $0$&$0$\\
& & & $(1,1)_{0}$ & $A^{B-L}_\mu$ & $0$ & $0$ &$0$&$0$&$0$ \\
& & & $(3,1)_{2/3}$ & $X_\mu$ & $1/3$ & $0$ &$0$&$4/3$&$0$\\
& & & $(\bar{3},1)_{-2/3}$ & $X^\dagger_\mu$ & $-1/3$ & $0$ &$0$&$-4/3$&$0$\\\cline{2-7}
&$(\bf{1},\bf{3},\bf{1})_0$&$(+,+)$&$(1,3)_0$&$W_\mu$& $0$ & $0$ &$0$&$0$&$0$\\ \cline{2-7}
&$(\bf{1},\bf{1},\bf{3})_0$&$(+,+)$&$(1,1)_1$&$W^{R+}_\mu$& $0$ & $0$ &$2$&$0$&$0$\\
&&&$(1,1)_{-1}$&$W_\mu^{R-}$& $0$ & $0$ &$-2$&$0$&$0$\\
&&&$(1,1)_{0}$&$A^R_\mu$& $0$ & $0$ &$0$&$0$&$0$\\ \cline{2-7}
&$(\bf{1},\bf{1},\bf{1})_0$&$(+,+)$&$(1,1)_0$&$A^\psi_\mu$& $0$ & $0$ &$0$&$0$&$0$\\ \cline{2-7}
&$(\bf{6},\bf{2},\bf{2})_0$&$(+,-)$&$(\bar{3},2)_{5/6}$&& $1/6$ & $1/2$ &$1$&$2/3$&$0$\\
&&&$(\bar{3},2)_{-1/6}$&& $1/6$ & $-1/2$ &$-1$&$2/3$&$0$\\
&&&$(3,2)_{1/6}$&& $-1/6$ & $1/2$ &$1$&$-2/3$&$0$\\
&&&$(\bar{3},2)_{-5/6}$&& $-1/6$ & $-1/2$ &$-1$&$-2/3$&$0$\\ \cline{2-7}
&$(\bf{4},\bf{2},\bf{1})_{-3}$&$(-,+)$&$(3,2)_{1/6}$& & $-1/6$ & $1/2$ &$0$&$1/3$&$-3$\\
&&&$(1,2)_{-1/2}$& & $-1/2$ & $1/2$ &$0$&$-1$&$-3$\\ \cline{2-7}
&$(\bar{\bf{4}},\bf{2},\bf{1})_{3}$&$(-,+)$&$(\bar{3},2)_{-1/6}$& & $1/6$ & $-1/2$ &$0$&$-1/3$&$3$\\
&&&$(1,2)_{1/2}$& & $1/2$ & $-1/2$ &$0$&$1$&$3$\\ \cline{2-7}
&$(\bar{\bf{4}},\bf{1},\bf{2})_{-3}$&$(-,-)$&$(\bar{3},1)_{1/3}$& $b_R^c$ & $-1/3$ & $1$ &$1$&$-1/3$&$-3$\\
&&&$(\bar{3},1)_{-2/3}$& $t_R^c$ & $-1/3$ & $0$ &$-1$&$-1/3$&$-3$\\
&&&$(1,1)_{1}$& $\tau_R^c$ & $0$ & $1$ &$1$&$1$&$-3$\\
&&&$(1,1)_{0}$& $\nu_R^c$ & $0$ & $0$ &$-1$&$1$&$-3$\\\cline{2-7}
&$(\bf{4},\bf{1},\bf{2})_{3}$&$(-,-)$&$(3,1)_{-1/3}$& $b'_L$ & $1/3$ & $-1$ &$-1$&$1/3$&$3$\\
&&&$(3,1)_{2/3}$& $t'_L$ & $1/3$ & $0$ &$1$&$1/3$&$3$\\
&&&$(1,1)_{-1}$& $\tau'_L$ & $0$ & $-1$ &$-1$&$-1$&$3$\\
&&&$(1,1)_{0}$& $\nu'_L$ & $0$ & $0$ &$1$&$-1$&$3$\\
\hline
\end{tabular}
\caption{\label{tab:5Dgauge} Field decomposition for the gauge multiplet $R={\bf 78}$.
Fields with orbifold parities $(+,+)$ have a zero mode in $W^\alpha_{78}$, while $(-,-)$ entails a zero mode in $\Phi_{78}$.}
\end{table}

\end{document}

redpajama_set_name: RedPajamaArXiv


text:
\section{Introduction}

Empirical evidence has shown that the dynamics of inter-event
times driven by human actions may not be random and not well
approximated by Poisson processes~\cite{paxflo96,masmon03,sca06}.
Based on this, Barab\'{a}si~\cite{bar05} developed a very
interesting model of human activity where the distribution of the
inter-event time is a consequence of a decision queue processes.
He considers that among the most relevant protocols for driving
human dynamics, e.g. first-in-first-out protocol, random protocol
and a protocol based on the execution of the high priority item,
this later protocol seems to be the most important. In this
protocol, while high priority tasks are executed as soon as they
are added to the list, low priority tasks wait for a long time
until all high priority tasks are executed, i.e., the instants of
execution of low priority tasks are separated by long times of
inactivity. Using this assumption, it is numerically~\cite{bar05}
(analytically~\cite{vaz05})  shown that the distribution of
inter-event times follows a power law.

Two interesting contributions were introduced by~\cite{grilin06}.
First the authors map the variable length priority model
considered above onto a model of biased diffusion deriving
asymptotic distributions for the inter-event times. Second, in
order to investigate the arising of power laws in more general
situations, they generalize the fixed length model queue to
contain tasks with a priority label and with a class label where
there is always an active class and an inactive class. If the
highest priority task of the inactive class exceeds that of the
active class by at least a fixed switching cost, the inactive
class becomes active and the active class becomes inactive.

An interesting discussion is considered in~\cite{ken06,baroli06}
where it is argued that other mechanisms contribute for the
distributions of waiting times such as deadlines, time dependence
of priorities and the social context of the problem. In line with
this debate, \cite{bla06} relaxes the assumption that the
priorities of tasks do not change over time and studies queueing
systems where deadlines are assigned to the incoming tasks and the
urgency to attend a task increases with time showing that only in
the former
 model fat tails arise naturally as consequence of the scheduling
 rule.

In this letter, we investigate the assumption that people execute tasks on
a protocol that execute firstly the high priority item. In
particular, we suppose that people assign priorities to the tasks
on their lists in order to minimize some cost index, i.e., a cost
associated to the fact of not processing a given collection of
tasks in a given time step. Therefore, based on this assumption
and inspired on~\cite{bar05,vaz05,grilin06}, we have built a
discounted stochastic dynamic programming model with two types of tasks
(low and high priority tasks) and a cost per stage for keeping a
number of low and high priority tasks without processing.

This is not the first time that a kind of optimization principle
is used to understand the structure and dynamics of complex
systems. In~\cite{rodrin92,caj05}, for instance, it is shown that
complex networks may arise from optimization principles.

It is also important to stress that although there is a large
literature dealing with control of queue
discipline~\cite{cra77,kitryk95,sen99} which this work is related,
the model presented in this letter is neither an extension nor a
particular case of any of these results.

We have found that the type of protocol used to execute tasks is
strongly dependent on the kind of instantaneous cost of keeping a
task in the queue for an additional stage. When linear costs are
used the protocol of executing preferentially the high priority
costs always is the best solution. However, this does not happen
when quadratic costs are considered. In this case different types
of protocol are considered. Furthermore, depending on the
parameters of the system, the protocol considered generates
complex one dimensional dynamics.

\section{Setup of the problem} We consider that there are two queues waiting for a service on a single server.
Let $g(x_L,x_H)$ be the current cost of having state $(x_L,x_H)$
which is the state of the system, $x_{L}$ ($x_H$) is the number of
tasks in the first (second) queue. We say that the first queue is
a low priority queue (or the second queue is a high priority
queue) if $\frac{\partial g(x_L,x_H)}{\partial
x_L}|_{x_L=x_H}<\frac{\partial g(x_L,x_H)}{\partial
x_H}|_{x_L=x_H}$. We assume that this is the case. The dynamics of
these queues are modeled as follows: At each discrete time step
with probability $\lambda\rho$ a new task arrives in the queue
formed by high priority tasks and with probability $\lambda
(1-\rho)$ a new task arrives in the queue formed by low priority
tasks. Within each of the queues the tasks are executed on a
First-In, First-Out basis. With probability $\mu u(x_L,x_H)$ the
first task of the high priority queue is executed and with
probability $\mu (1-u(x_L,x_H))$ the first task of the low
priority queue is executed. We assume here that $u(x_L,x_H)$ is a
state dependent control variable that the agent will choose in
order to minimize the total cost function
$J_u(x_L,x_H)=E_{x_L,x_H}^{u}[\sum_{t=1}^{\infty}\alpha^t
g((x_L(t),x_H(t)))]$, where $E_{x_L,x_H}^{u}[\cdot]$ is the
expected value conditioned to the current state $(x_L,x_H)$ and
the state control variable $u$ and $\alpha$ is the discount
factor.

Due to the principle of optimality~\cite{bel57,ber01} and the
Banach fixed point theorem, if the minimum cost function
$J(x_L,x_H)=\min_{u(x_L,x_H)\in [0,1]}J_u(x_L,x_H)$ exists, it
must be given by the unique solution of the Bellman equation, that
may be written as
\begin{equation}J(x_L,x_H)=F(x_L,x_H)+\min_{u(x_L,x_H)\in [0,1]}u(x_L,x_H)G(x_L,x_H)\label{bellmanEquation} \end{equation}
where
\begin{eqnarray}F(x_L,x_H)&=&g(x_L,x_H)+ \lambda  \rho  (1-\mu)  [\alpha
  J(x_L,x_H+1)]\nonumber\\
  &+&\lambda  (1-\rho)  (1-\mu)  [\alpha J(x_L+1,x_H)]\nonumber\\
  &+&(1-\lambda)  \mu  [\alpha  J(x_L-1,x_H)]\nonumber\\&+&\rho  \lambda  \mu  [\alpha  J(x_L-1,x_H+1)]\nonumber\\
  &+& (1-\rho)  \lambda  \mu  [\alpha  J(x_L,x_H)]\nonumber\\
  &+& (1-\lambda)  (1-\mu)  [\alpha  J(x_L,x_H)] \end{eqnarray}

and

\begin{eqnarray}&&G(x_L,x_H)=\nonumber\\&&(1-\lambda)  \mu
[\alpha
    (J(x_L,x_H-1)-J(x_L-1,x_H))]\nonumber\\
    &+&  \rho  \lambda  \mu  [\alpha
    (J(x_L,x_H)-J(x_L-1,x_H+1))]\nonumber\\
    &+& (1-\rho)  \lambda  \mu  [\alpha
    (J(x_L+1,x_H-1)-J(x_L,x_H))] \label{eq:G}\end{eqnarray}

Since the optimization problem (\ref{bellmanEquation}) is a linear
programming problem, the optimal control $u(x_L,x_H)$ in each state
$(x_L,x_H)$ will depend explicitly on the signal of $G(x_L,x_H)$.
If $G(x_L,x_H)>0$, then $u(x_L,x_H)=0$. If $G(x_L,x_H)<0$, then
$u(x_L,x_H)=1$. Finally, if $G(x_L,x_H)=0$, $u(x_L,x_H)$ is a
mixed strategy that may present any value in the interval $[0,1]$.
It is quite intuitive this result. Indeed, one may note that the terms in square brackets defined in
$G(x_L,x_H)$, equation (\ref{eq:G}), comprises the variations in the cost function due to changes in the states of the queue related to the execution of one of the tasks.

Since the properties of the solution of $J(x_L,x_H)$ of the
Bellman equation (\ref{bellmanEquation}) are strongly dependent on
choice of the cost per stage $g(x_L,x_H)$, in the next sections,
two different choices for $g(x_L,x_H)$ are investigated.

\section{Linear costs}

In this section, we assume that $g(x_L,x_H)=h_L x_L+h_H x_H$, for
$0<h_L<h_H$, i.e., the current cost of having one additional high
priority task in the queue is larger than having one additional
low priority task in the queue.

Since the space of polynomials of degree 1 with sup-norm is a
Banach space, one can show inductively, making recursive
iterations of the dynamic programming mapping, that $J(x_L,x_H)$
is also linear. Therefore, for $x_L>0$ and $x_H>0$, guessing this
form, one may easily solve the Bellman equation
(\ref{bellmanEquation}) and show that the cost function is given
by \footnote{The constants are given by
$c_L=\frac{h_L}{1-\alpha}$, $c_H=\frac{h_H}{1-\alpha}$ and
$c=\frac{\alpha}{(1-\alpha)^2} (\lambda (1-\rho)
h_L+(\lambda\rho-\mu)h_H )$.}
\begin{equation}J(x_L,x_H)=c+c_L x_L +c_H x_H\end{equation}

\vspace{2mm}

Furthermore,
\begin{equation}G(x_L,x_H)=\mu \frac{\alpha}{1-\alpha} (h_L - h_H) \end{equation}
is always negative implying that $u(x_L,x_H)=u=1$ for every state
$(x_L,x_H)$. Therefore, if linear costs are considered, the
protocol to be considered is the one based on the execution of the
high priority task whenever there is at least one item in this
queue, i.e., $x_H>0$. This kind of protocol was very well studied
in~\cite{bar05,vaz05,grilin06} where analytic results for the
emerging of power laws may be found. In the next section, a much
wealthier situation happens where the optimal policy is not only
limited to execute the high priority item in the queue, but the
optimal policy is state-dependent.

\section{Quadratic costs} Now, we assume that $g(x_L,x_H)=h_L x_L^2+h_H
x_H^2$, for
$0<h_L<h_H$. Following the same reasoning already presented before for the
linear cost case, one may conclude a quadratic form for the cost
function.

Solving the Bellman equation, one may show that the solution of
the problem depends explicitly on the signal of the function
$G(x_L,x_H)$, defined in (\ref{eq:G}), in the state $(x_L,x_H)$.
In fact, three different regions will arise. We will call region
$A$ the domain of $(x_L,x_H)$ where $G(x_L,x_H)>0$, region $B$ the
domain of $(x_L,x_H)$ where $G(x_L,x_H)=0$ and region $C$ the
domain of $(x_L,x_H)$ where $G(x_L,x_H)<0$.

We have found that the minimum cost function is given by

\begin{equation}J(x_L,x_H)=\left\{\begin{array}{c}
                                    J^A(x_L,x_H)\;\mathrm{if}\; (x_L,x_H)\in A \\
                                    J^B(x_L,x_H)\;\mathrm{if}\; (x_L,x_H)\in B \\
                                    J^C(x_L,x_H)\;\mathrm{if}\; (x_L,x_H)\in C \\
                                  \end{array}\right.\label{J} \end{equation}
and the optimal control is given by

\begin{equation}u(x_L,x_H)=\left\{\begin{array}{c}
                                    0\;\mathrm{if}\; (x_L,x_H)\in A \\
                                    u\in [0,1]\;\mathrm{if}\; (x_L,x_H)\in B \\
                                    1\;\mathrm{if}\; (x_L,x_H)\in C \\
                                  \end{array}\right.\label{u} \end{equation}

where for $i=A,B,C$
\begin{eqnarray}J^i(x_L,x_H)&=&c^{i}+c_{L1}^{i}x_L + c_{H1}^i x_H +
c_{L2} x_{L}^{2} + c_{H2}x_{H}^{2}\nonumber\\\end{eqnarray}

\begin{eqnarray}G^i(x_L,x_H)&=&\mu[c_{L1}^{i} -c_{H1}^{i}+2 (c_{L2} x_l -c_{H2} x_H)\nonumber\\
&+& c_{L2} (2\lambda (1-\rho) -1)
+c_{H2}(1-2\lambda\rho)]\nonumber\\
\end{eqnarray}
 and $G(x_L,x_H)=G^i(x_L,x_H)$, if $(x_L,x_H)\in i$~\footnote{The constants are given by $c_{L2}=\frac{h_L}{1-\alpha}$,
$c_{H2}=\frac{h_H}{1-\alpha}$, $c_{L1}^{A}=\frac{2\alpha h_L (-\mu
+ \lambda (1-\rho))}{(1-\alpha)^2}$,
$c_{L1}^{B}(\delta)=\frac{1}{(h_L+h_H)}\{\frac{h_{L}^{2}}{(1-\alpha)}[1-
2\lambda(1-\rho)]
 +\frac{h_L h_H}{(1-\alpha)^2} [2\alpha(\lambda-\mu)
+(2\lambda\rho-1+2\delta)(1-\alpha)]$, $c_{L1}^{C}=\frac{2\alpha
h_L \lambda (1-\rho)}{(1-\alpha)^2} $, $c_{H1}^{A}=\frac{2\alpha
h_H \lambda \rho}{(1-\alpha)^2}$,
$c_{H1}^{B}(\delta)=\frac{1}{(h_L+h_H)}\{\frac{h_{H}^{2}}{(1-\alpha)}[1-
2\lambda\rho-2\delta] +\frac{h_L h_H}{(1-\alpha)^2}
[2(\lambda-\mu\alpha) -(2\lambda\rho+1)(1-\alpha)]\}$,
$c_{H1}^{C}=\frac{2\alpha h_H (-\mu+\lambda\rho)}{(1-\alpha)^2} $,
$c^A=\frac{\alpha}{(1-\alpha)^3}\{[(1-\alpha)(\mu+(1-\rho)\lambda))+2\alpha(\lambda^2+\mu^2
-\mu\lambda)+2\rho\lambda(\mu(1+\alpha)+\lambda\alpha(\rho-2))-2\mu\lambda]h_L
+ \lambda\rho[1-\alpha+2\lambda\rho\alpha]h_H$,
$c^B=\frac{\alpha}{(h_L+h_H)(1-\alpha)^3}\{[2\lambda\alpha(1-\alpha)(-\rho^2
\lambda+2\lambda\rho-1-\lambda-\rho) ]h_{L}^2
+[\delta(2\lambda\rho\alpha^2+4\lambda\alpha-4\lambda\rho\alpha+1-2\lambda-2\alpha+\alpha^2-2\lambda\alpha^2+2\lambda\rho)+2\lambda(-\mu\alpha+2\lambda\rho^2\alpha^2+2\lambda\rho\alpha
-2\lambda\rho^2\alpha-\alpha^2
\mu-2\lambda\alpha^2\rho+\lambda\alpha^2)+2\mu(\alpha-\alpha^2
+\mu\alpha^2)]h_L h_H
+[\delta(-4\lambda\rho\alpha+2\delta-\alpha^2+2\delta\alpha^2+2\lambda\rho\alpha^2
-1
-4\alpha\delta+2\alpha+2\lambda\rho)+2\lambda\rho(\lambda\rho\alpha^2-\lambda\rho\alpha-\alpha^2+\alpha)]h_{H}^2
\}$ and $c^C=\frac{\alpha}{(1-\alpha)^3}\{[\lambda
(1-\alpha)(1-\rho)+2\lambda^2\alpha (1-\rho)^2]h_L
+[(1-\alpha)(\mu+\lambda\rho)+2\mu^2\alpha +2\lambda\rho
(\lambda\rho\alpha -\mu\alpha -\mu]h_H \}$.}.

\vspace{2mm}

Furthermore, the parameter
$\delta\in[\underline\delta,\overline\delta]$, defines the set of
points of $\Re^2$ such that
$G(x_L,(h_L/h_H)x_L+\delta)=0$~\footnote{The constants
\underline{$\delta$} $\,$ and $\overline\delta$ are respectively
given by
\underline{$\delta$}$=\frac{1}{1-\alpha}\{[\frac{1-\alpha}{2}-\lambda\rho]+
\frac{h_L}{h_H}[-\frac{1-\alpha}{2}+\lambda(1-\rho)-\alpha\mu]\}$
and
$\overline{\delta}=\frac{1}{1-\alpha}\{[\frac{1-\alpha}{2}-\lambda\rho+\alpha\mu]
+ \frac{h_L}{h_H}[-\frac{1-\alpha}{2}+\lambda(1-\rho)]\}$.}.

Indeed,
\begin{equation}\delta\rightarrow \underline{\delta} \Rightarrow \left(J^B(\delta)\rightarrow J^A,\;G^B(\delta) \rightarrow G^A \right)\end{equation}
and
\begin{equation}\delta\rightarrow \overline{\delta} \Rightarrow \left(J^B(\delta)\rightarrow J^C,\;G^B(\delta) \rightarrow G^C \right)\end{equation}

In order to understand the intuition behind this solution, one is
invited to consider a particular case where $h_L/h_H\approx 1$,
$\lambda=\mu=1$ and $\rho=1/2$. In this case, the region $B$ is
defined by the set of points where $|x_H-x_L|<\alpha/(1-\alpha)$.
Therefore, the set of points where the decision maker can use any
strategy depends strictly on the discount factor. If the discount
factor is large, the decision maker may keep queues with  a large
difference between their sizes. On the other hand, if the discount
factor is small this situation is not accepted as a solution
anymore.

\begin{figure}
\begin{picture}(100,100)(-30,0)
\thicklines \unitlength=0.6mm
\put(40,10){\vector(0,1){80}} \put(40,10){\vector(1,0){90}}
 \put(135,10){$x_L$}
\put(30,95){$x_H$} \put(40,33){\line(3,1){75}}
\put(30,33){$\overline{\delta}$}\put(92,61){$x_H=(h_L/h_H)x_L+\overline{\delta}$}
\put(70,10){\line(3,1){75}}
\put(122,40){$x_H=(h_L/h_H)x_L+\underline{\delta}$}
\put(112,17){$A$} \put(80,25){$B$} \put(50,50){$C$}
\put(60,0){$-(h_H/h_L)\underline{\delta}$}
\end{picture}

\caption{The regions $A$, $B$ and $C$ in the plane $x_L - x_H$.}\label{fig0}
\end{figure}

Differently from the linear costs case, several types of protocol
are possible. Region $C$ considers a protocol based on the
execution of the high priority task. Region $A$ considers a
protocol based on the execution of the low priority task. It
occurs in order to avoid that the size of the queue of the low
priority tasks do not increase too much. ``Too much'' here is
measured by the ratio $h_L/h_H$. Region $B$ does not determine a
protocol. It can be a random protocol (mixed strategy) or simply a
protocol such the one considered in region $C$ or region $A$.
Figure \ref{fig0} shows the geometry of these regions in the plane
$x_L - x_H$.

It is not difficult to show that the expected value of the state
obeys the following dynamics

\begin{eqnarray}E_t[x(t&+&1)]=E_t\left[\begin{array}{c} x_L(t+1) \\ x_H(t+1) \\\end{array}\right]=\left[\begin{array}{c} x_L(t) \\ x_H(t) \\\end{array}\right]\nonumber\\&+&\left[\begin{array}{c} \lambda (1-\rho)-\mu (1-u(x_L(t),x_H(t))) \\ \lambda \rho-\mu u(x_L(t),x_H(t)) \\\end{array}\right] \end{eqnarray}

which has infinite fixed points if and only if $\lambda=\mu$ and
$u(x_L,x_H)=u=\rho$.

We will analyze only the most interesting situation which is the
fixed-length-queue, i.e., $\lambda=\mu$. Therefore, assuming that
$\lambda=\mu$, $u^B=\rho+\epsilon$ and $\epsilon>0$, then the
expected value of the system is governed by
$E_t[x(t+1)]=x(t)+\lambda\epsilon e$ if it is in region $B$ and by
$E_t[x(t+1)]=x(t)-\lambda\rho e$ if it is in region $A$ (if the
state is in region $C$, the expected state will certainly come to
region $B$ and not come back to this region), where $e=(1,-1)'$.
Therefore, the dynamics takes place in the line passing by $x(0)$
and following the direction $e$. Thus, if the expected state is in
region $B$ it goes into the direction of region $A$ and viceversa.
This dynamics is equivalent to the one dimensional system

\begin{equation}y(t+1)=\left\{\begin{array}{c}
                                y(t)+t^{+}\;\mathrm{if}\; y(t)\le 0 \\
                                y(t)-t^{-}\;\mathrm{if}\; y(t)>0 \\
                              \end{array} \right. \label{dynamics}\end{equation}
defined on the interval $(-t^{-},t^{+}]$, where
$t^{+}=\lambda\epsilon$ and
$t^{-}=\lambda\rho$

\begin{figure}
\begin{center}
\includegraphics[width=8cm,height=8cm]{fig1.eps}
\caption{The evolution of $y(t)$ for $y_0=-0.2$, $\lambda=0.5$,
$\rho=0.5$ and $\epsilon=0.3$.} \label{fig1}
\end{center}
\end{figure}

The dynamics of this system is plotted in figure \ref{fig1} for
the case of $t^{-}=0.25$ and $t^{+}=0.15$. The dynamics defined in
(\ref{dynamics}) is topologically conjugate with the translation
in the circle~\cite{dem93}. Therefore, if $t^+/t^-=p/q$, where
$p/q$ is a irreducible ratio representation of rational number,
then this system follows a limit cycle with period $p+q$.
Otherwise, the $\omega$-limit of any point in the interval is a
dense subset of it. Therefore, we can conclude that the stochastic
process that defines the length of each queue is not stationary.
Moreover, the dynamics of the expected value of the length of the
queue exhibits a complex behavior: infinitely many cycles or a
$\omega$-limit set being a dense subset in the interval. The
intuition behind that complex dynamics is quite reasonable. In the
region close to the frontier $x_H=(h_L/h_H)x_L
+\underline{\delta}$ that separate $A$ and $B$ (see figure
\ref{fig0}), we can observe the following: if the expected state
is in $A$, its dynamics moves toward region B, since the priority
is of $L$ . Once the expected state is in $B$, the dynamics takes
it back to the region A, since in this case in average the
priority is of $H$ (due to the condition $u^B=\rho+\epsilon$ and
$\epsilon>0$). Because the frequency of tasks arriving is equal to
that of attending them, a cyclical or complex dynamics emerges
close to the referred frontier. Figure \ref{fig1} shows the case
where this system is a limit cycle. A similar situation involving
regions $B$ and $C$ arises in the case of $\epsilon>0$ and
$\rho=u^B+\epsilon$. In these situations, the protocol is ruled by
the protocols considered in regions $A$ and $B$ in the former case
and by the protocols considered in regions $B$ and $C$ in the
later case.

For $\lambda\neq \mu$, either the expected value goes to infinite,
converges to 0, to axis $x_L=0$ or to axis $x_H=0$, following
different routes. Furthermore, different kinds of protocols
are possible.
\section{Final Remarks}In the human dynamics of the tasks execution decisions the priority of
one task is not always defined as being the most important current
task. Actually, the dynamics of the work executions depends on the
cumulated tasks of short run priorities, the importance of each
kind of task and the intertemporal discount factor. In this
letter, we provide a stochastic dynamic programming model
containing all those elements and analyze the dynamics of the
execution of tasks, shedding new light to the discursion
considered in~\cite{ken06,baroli06}. In this setting, we have
found that the dynamics of the expected state of the system may be
complex, exhibiting cycles of any order or with limit set being a
dense subset of the interval depending on the parameter values of
the model. This is a contribution to a better understanding of how
human dynamics may evolve in this type of problem. Finally, it is
worth noting that complex dynamics in the solution of dynamic
programming problems are usually obtained for low discount
factors~\cite{monsor96}. However, in our quadratic case, complex
dynamics arises for discount factors of any size.

\section{Acknowledgment} The authors are indebt to the Brazilian
agency CNPQ for financial support. The first author is also indebt
to Professor G. Grinstein (at IBM Watson Research Center) who
helped him by means of a personal communication to understand some
details of his paper~\cite{grilin06}.