Modeling Volatility and Dependence of European Carbon and Energy Prices

Modeling Volatility and Dependence of European Carbon and Energy Prices

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 power 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: 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 chevallier2019conditional, who find a negative link between carbon prices and oil and gas prices. meier2020commodities also find that EUA prices are significantly influenced by fossil fuel prices, while zhou2020measurements show evidence that the carbon price drives the price of renewable energies. Employing an artificial neural network, 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). 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 on forecasting EUA prices, which is particularly relevant for portfolio management. paolella2008econometric look into forecastability and value-at-risk forecasting in the EU-ETS, employing an AR-GARCH model with generalized skewed student-t innovations. benz2009modeling use a Markov-switching framework including GARCH effects to model EUA prices, while 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 trabelsi2022co2, benz2009modeling, paolella2008econometric.

Another strand of the literature deals with the volatility process of EUA prices and possible spillovers to other markets. 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. hanif2021nonlinear investigate volatility spillovers between EUA prices and energy indices using copulas. Finally, 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 heteroscedasticity in the individual series as well as for \color{black}{{time-varying}} \color{black} dependence structures between them in a simultaneous yet flexible way. In contrast to earlier combination models which are estimated in stages, we utilize 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. \color{black}{{Second, we evaluate the proposed model in an extensive rolling-window forecasting study for both data in levels and after log-transformation, whose properties differ (e.g. variance, cointegration).}} \color{black} 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 \color{black}{{against competing standard benchmarks. Finally, we consider real price data, normalized by the carbon emissions of the respective commodity, to allow for scale-invariant evaluation.}} \color{black}

The remainder of this paper is organized as follows. The data are presented in Section (ref). Section (ref) 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), while Section (ref) discusses the results. Finally, Section (ref) concludes.

Data description and preliminary analysis

The data, plotted in Figure (ref), comprise 3257 daily observations of short-term futures prices covering the period from March 15, 2010, until October 14, 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 crude oil. Oil and coal prices, originally given in USD, have been converted from USD into EUR. \color{black}{{All depicted price series have been standardized by the CO$_2$ emissions of the respective commodity (for conversion factors see quaschning2016understanding). This means that the emission-adjusted prices reflect one tonne of CO$_2$ emissions for each of the three commodities. EUA prices needed no adjustment since they already reflect one tonne of CO$_2$ emissions. Furthermore, we adjusted for inflation by Eurostat's HICP, excluding energy\footnote{Source: \href{https://ec.europa.eu/eurostat/databrowser/bookmark/13937a8d-e341-406d-aacb-a6c34f93c175?lang=en}{Eurostat} (HICP - monthly data, Overall index excluding energy: TOT_X_NRG, Euro Area: EA, Base Year 2015: I15)}. Thus the general price increase over time is accounted for without dampening the commodities' price volatility which is to be explained by the model.}} \color{black}

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 \color{black}{{and logarithms of all four series, respectively,}} \color{black} and stationarity in their first differences\footnote{Results of unit root and cointegration tests are omitted for brevity but are available upon request.}. \color{black}{{As for cointegration, Johansen's likelihood ratio trace test points to two cointegrating relationships in the system in levels, which vanish in the log-transformed data due to the non-linear nature of the transformation. Subsequently, both data in levels and in logarithms have been considered for the estimation and evaluation of the proposed model.}} \color{black} This allows us to evaluate whether accounting for long-term equilibrium relationships contributes to better forecasting performance.

The VECM-Copula-GARCH model

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 \color{black}{{and cointegrated in the levels.}} \color{black} While cointegration is not formally detected in the \color{black}{{log-}} \color{black}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 \color{black}{{models as in e.g. jondeau2006copula}} \color{black}.

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 \color{black}{{ $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 $F_{\boldsymbol X_t|\mathcal{F}_{t-1}}$ of the model variables in terms of their marginal distributions $F_{X_{k,t}|\mathcal{F}_{t-1}}$ for $k=1,\ldots, K$, and the copula $C_{\boldsymbol U_{t}|\mathcal{F}_{t - 1}}$,}} \color{black} as is popular in many financial applications jondeau2006copula, hu2006dependence. Thus, for the realized values $\boldsymbol x_t= (x_{1,t},\ldots, x_{K,t})^\intercal$ it holds that \color{black}{{

}} \color{black} where $\boldsymbol u_t=(u_{1,t},\ldots, u_{K,t})^\intercal$ and $u_{k,t} = F_{X_{k,t}|\mathcal{F}_{t-1}}(x_{k,t})$, $k=1, \ldots, K$. That is, we specify $K$ marginal model components along with a copula that captures their dependency structure. In order to simplify the notation we drop the conditioning on $\mathcal{F}_{t-1}$ in equation (ref) and write henceforth $F(\boldsymbol x_t) = C(\boldsymbol u_t)$.

The model can thus be specified as follows: \color{black}{{

}} \color{black} where the \color{black}{{$(K(K-1)/2)$-dimensional}} \color{black} 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 $C$ as 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. (ref), 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 of freedom and non-centrality parameters gathered in \boldmath$\nu$ $= (\nu_1, \ldots, \nu_K)^{\intercal}$ and \boldmath$\lambda$ $= (\lambda_1, \ldots, \lambda_K)^{\intercal}$. For estimation, each $F_i$ has been parameterized so 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 \color{black}{{

}} \color{black} where $\Pi = \alpha \beta^{\intercal}$, $\alpha \in {\mathbb R}^{K \times r}$ and $\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:

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 $K(K-1)/2$ primary dependence parameters of the copula, gathered in the (scaled) correlation matrix $\Xi_t$, is modeled similarly to a GARCH process:

where $\xi_{ij,t}$ is a latent process, $z_{i,t}$ denotes the $i$-th standardized residual from time series $i$ \color{black}{{at time point $t$}} \color{black}, and $\Lambda(\cdot)$ is a link function. The latter 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 $10^{-5}$.

Relevant nested models

The proposed model (ref) is the most general model incorporating cointegration, conditional heteroscedasticity, and time-varying dependence parameters. This model is denoted by $\text{VECM}_{\text{lev}, \text{ncp}}^{r, \sigma_t, \rho_t}$. \color{black}{{Here the subscripts $\textbf{lev}$ and $\textbf{ncp}$ denote joint estimation along with the other model parameters of the leverage and non-centrality parameters, respectively,}} \color{black} 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. $\text{RW}$ will stand for a model with no cointegration such 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., $\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. \color{black}{{Models estimated on log-transformed data are signified by the additional subscript $\text{log}$.}} \color{black}

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

with the multivariate conditional density at any time point $t$ given by

Here the copula density $c$ can be derived analytically from the copula as $c(\mathbf{u}) = \partial_{\mathbf{u}} C(\mathbf{u})$.

Forecasting study

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$}_{t}^2$ are determined using the information contained in $\mathcal{F}_{t}$. Then $n$ samples are drawn from the $K$-dimensional 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$}_{t}^2, \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.

Forecast evaluation

Evaluation of the probabilistic forecast is conducted using different evaluation measures. The probabilistic multivariate forecasts will be evaluated by the energy score (ES). The energy score of a multivariate probabilistic forecast with distribution $F$ is given by

where $\mathbf{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$ gneiting2007strictly. A low value of the ES indicates a good probabilistic forecast. In the univariate case, the energy score becomes the Continuous Ranked Probability Score (CRPS). \color{black}{{Mean forecasts are evaluated by the root mean squared error (RMSE). Significant differences in forecast performances are tested with the Diebold-Mariano test diebold2002comparing with the adjustment by harvey1997testing.}} \color{black}

Results

\color{black}{{We consider a rolling-window forecasting study with a window size of 1000 days which corresponds to about four years of data. As we have 3257 observations in total and compute 30-steps-ahead forecasts, this leaves 2227 potential starting points for the rolling-window study. To reduce computational costs, we sample $n=250$ data points uniformly from those 2227 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.}} \color{black}

We evaluate the full spectrum of the aforementioned VECM and RW models in conjunction with a \color{black}{{simple}} \color{black} exponential smoothing model, denoted by $\text{ETS}^{\sigma}$, \color{black}{{and a vector exponential smoothing model \color{black}{{svetunkov2023new}} \color{black}, denoted by $\text{VES}^{\sigma}$, }} \color{black} which we include as benchmarks. \color{black}{{We choose a univariate and a multivariate ETS model as benchmarks since they are commonly used for forecasting jonsson2014exponential, abdul2019comparison, berrisch2022distributional.}} \color{black} For the sake of brevity, however, we report the performance only of selected ones, particularly the best-performing models. \color{black}{{The performance of the models estimated on log-transformed data has been evaluated by taking the exponent of the forecasts without bias correction, in accordance with the recommendation of Demetrescu2020 for persistent data.}} \color{black} \color{black}{{However, the exponent of some log forecasts evaluates to infinity. We replace these values with a large real integer ($10^5$) to enable a valid evaluation of these forecasts. This issue affects less than 0.001\textperthousand\ of all forecasts of the reported log models. We used $10^5$ as a replacement because it can be considered an extreme scenario for the prices in question. Note that we deal with inflation-adjusted, emission-normalized prices here. Evaluating the research on the social costs of carbon confirms that $10^5$ would be an extreme scenario, at least for EUA anthoff2013uncertainty, tol2019social, aldy2021keep. The latter also holds for the other commodities if other price components are negligible.}} \color{black}

The performance of the selected models with respect to the ES can be found in Table (ref). Their ES is compared to the ES of a simple random walk model with constant copula and no conditional heteroscedasticity modeling, $\text{RW}^{\sigma, \rho}$, \color{black}{{estimated on the data in levels.}} \color{black} $\text{ES}_{1-30}^{\text{All}}$ denotes the score of the $(4 \times 30)$-dimensional multivariate ($4$-dimensional) 1-30 days ahead forecast. The 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.

The results of Table (ref) can be summarized as follows. \color{black}{{No model dominates uniformly in terms of forecasting performance for all variables and all horizons. Models estimated on the log-transformed data generally perform better than their counterparts estimated on the data in levels, despite the na{\"\i}ve forecasts of the former being potentially biased. This result can be explained by the log-transformation successfully stabilizing the variance of the series}} \color{black} Lutkepohl2012. \color{black}{{As no cointegration was found in the log-transformed data, the VECM models are over-parameterized and are hence not expected to deliver superior forecasts. Indeed, considering all variables and all forecast horizons altogether, the best-performing models yield an improvement of ca. 12% in $\text{ES}^{\text{All}}_{1-30}$ over the simple RW model in levels. This is achieved by the $\text{VECM}^{r0}$ and the RW models with constant volatility and log transformation. Allowing for time-varying dependence in the copula, or modeling the non-centrality parameter does not seem to influence their results much. These models also deliver the best forecasts at long horizons, and for natural gas prices in particular.}} \color{black} \color{black}{{VECM models in levels provide slightly worse overall forecasts than the RW in logarithms but seem to perform best of all at short horizons with $r=3$. Nevertheless, accounting for cointegrating relations offers at most modest improvement over the corresponding RW models in levels. Rather, allowing for time variation in the conditional variance and dependence parameters in the levels-RW models improves the forecast significantly and comparably to considering a simple RW in logarithms for $H=1,5$. The contribution of the leverage effects and the non-centrality parameter for the models in levels is subordinate. The forecasts of the level VECMs, however, seem quite unsatisfactory for the EUA and oil price series. They are rather best forecast by RW and ETS: the EUA price in logarithms, while the oil price in levels.}} \color{black}

Table (ref) presents the CRPS of the selected models for every univariate time series. H1, H5, and H30 denote the score \color{black}{{or improvement thereof}} \color{black} for the 1-, 5- and 30-day-ahead probabilistic forecast, respectively. The CRPS gives information about the quality of the univariate modeling only, as modeling the cross-sectional dependence has no effect on the univariate density 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) align with those of the ES evaluation: \color{black}{{There is no uniformly superior model in general. The RW and $\text{VECM}^{r0}$ models seem to perform quite well also in terms of modeling the marginal distributions: in levels for oil and gas prices, in logarithms for EUA prices. For natural gas prices, the levels-VECM models outperform the rest at short horizons, and accounting for leverage effects in the GARCH-part of the model is beneficial.}} \color{black}

Table (ref) presents the \color{black}{{mean}} \color{black} forecast performance in terms of the RMSE. \color{black}{{Considering that no model significantly outperforms the level RW models we find evidence supporting the efficient markets hypothesis.}} \color{black} All RW models perform similarly as distributional modeling has no influence on the mean forecast. Rather, the differences emerge from random simulation fluctuations. \color{black}{{The models in logarithms, however, seem to perform significantly worse than the level ones, in particular for the EUA prices. This may be attributed to the few extreme trajectories that required capping, as discussed above. }} \color{black}

\let\pgfimageWithoutPath\pgfimage \color{black}{{We next evaluate the temporal evolution of the copula-implied pairwise correlations of the best-performing (in terms of ES) model $\text{VECM}^{r0, \sigma, \rho_t}_{\text{log}}$ presented in Figure (ref).}} \color{black} The correlations between the prices of natural gas and coal, and EUA and coal are approximately constant over time. \color{black}{{The correlations of oil with coal and gas fluctuate around a level of 0.2 featuring a temporary increase of up to 0.5 in the emerging Russian invasion.}} \color{black} \color{black}{{The}} \color{black} most volatile correlation is that between EUA and natural gas prices. \color{black}{{It fluctuates strongly around a level of approximately 0.3.}} \color{black} Immediately after the Russian invasion of Ukraine, the correlation between EUA and gas prices changes dramatically, dropping from a level of about \color{black}{{0.3}} \color{black} to nearly \color{black}{{-0.4.}} \color{black} \color{black}{{By mid-March 2022 the relation has stabilized.}} \color{black}