Realized Volatility Forecasting in the Spanish Electricity Market During the 2021–2025 Energy Crisis

1. Introduction

The Iberian Electricity Market (MIBEL), established in 2007, has undergone significant structural transformation driven by market liberalization [1], rapid renewable penetration [2,3], and increased exposure to geopolitical and macroeconomic shocks [3,4]. The 2021–2025 period is particularly characterized by the European energy crisis, triggered by the Russia–Ukraine conflict, extreme fuel price fluctuations, and the accelerated integration of wind and solar generation in Spain [5]. These factors have substantially altered price formation in the Spanish electricity spot market, increasing volatility, amplifying extreme price movements, and intensifying the occurrence of zero and negative prices. A key manifestation of this transition is the “duck curve” phenomenon, where periods of high renewable output lead to sharp price collapses [6].

Volatility is a statistical measure of the variability of an asset’s returns over time. It is usually interpreted as an indicator of risk associated with price movements, and is a latent variable [7,8]. Therefore is a central object of analysis in electricity markets, as it directly affects risk management, hedging strategies [9], investment decisions, and regulatory oversight. Unlike traditional financial assets, electricity cannot be stored, and supply and demand must be balanced in real time [10,11]. Prices are driven by weather conditions, renewable intermittency, and regulatory interventions. This generates intraday and seasonal patterns, long-memory dependence, and abrupt jumps, which challenge standard econometric models.

The initial literature on volatility modeling is based on conditional heteroskedasticity models such as ARCH [12] and GARCH [13], along with extensions designed to capture asymmetries, structural breaks, and exogenous effects, including Markov-Switching GARCH [14], GARCH-X, and Taylor–Schwert-type specifications [15,16]. While used in energy markets, these models rely on low-frequency information and squared returns, limiting their ability to fully exploit high-frequency dynamics, especially under strong seasonality and jump behavior.

The generalization of intraday data facilitated the creation of realized volatility metrics, which offer non-parametric estimates of underlying volatility through high-frequency returns [17]. This procedure led to robust estimators designed to address jumps and microstructure noise, including Bipower Variation, Realized Kernel estimators [18,19] and Median Realized Volatility [20]. Within this literature, the Heterogeneous Autoregressive (HAR) model [21] offers a parsimonious representation of long-memory volatility dynamics through the aggregation of daily, weekly, and monthly components. In electricity markets, HAR-type models consistently show strong empirical performance relative to GARCH-based specifications, particularly when augmented with jump and calendar effects [22,23].

Despite these advances, several challenges remain. The presence of zero and negative prices prevents the direct application of logarithmic transformations, requiring appropriate preprocessing such as price shifting. Strong deterministic intraday and seasonal patterns also necessitate demeaning procedures, such as those proposed by [10]. In addition, realized volatility measures are subject to measurement error and jump contamination, which complicates estimation and forecast evaluation. As pointed out by [8], very short-term predictions using realized measures may be biased. Therefore, it is necessary to verify afterward if the models are biased.

Forecast comparison in this context is non-trivial due to model uncertainty, the risk of data-snooping and the complexity of analyzing a latent variable such as volatility. To address this, this study employs a loss based evaluation strategy using the QLIKE function [24], combined with the Model Confidence Set procedure [25], which identifies a set of statistically superior models with equal predictive ability. Additional diagnostic tools, including out-of-sample R² and Mincer–Zarnowitz regressions, are used to assess forecast properties such as bias and efficiency without affecting model selection.

This paper contributes to the literature by providing a comprehensive out-of-sample comparison of volatility forecasting models for the Spanish electricity spot market during the 2021–2025. The analysis considers a wide range of specifications, including GARCH-type models and HAR-based extensions with jump components, functional transformations, and calendar effects, using high-frequency data and robust realized volatility estimators.

The remainder of the paper is organized as follows. Section 2 describes the data and the main characteristics of the Spanish electricity market. Section 3 presents the volatility estimators, model specifications, and the forecasting evaluation. Section 4 reports and discusses the empirical results of the study. Section 5 concludes this study by summarising the main empirical findings and highlighting their implications for volatility modelling in electricity markets.

2. Data and Spanish Market Characteristics

This work uses hourly price data from January of 2021 to December of 2026 of the Spanish day-ahead electricity market, which is operated by the Spanish Market Operator (OMIE) as part of the Iberian Electricity Market (MIBEL). This market was established through a 2004 agreement between Spain and Portugal and has been functional since 2007. The market was created due to European liberalization reforms introduced by Directive 1996/92/EC and reinforced by Directive 2003/54/EC. These directives required the separation of generation, distribution, and supply activities. In MIBEL, OMIE manages the spot market, the Portuguese Market Operator (OMIP) oversees derivatives trading, and Red Eléctrica de España takes care of system operations.

Despite attempts to integrate, Spain has limited connections with Europe [26]. It is often called an “energy island. which leads to specific market dynamics.

An important structural feature during this period is the rapid growth of renewable energy. According to Red Eléctrica de España [27], renewables made up about 56% of total production in 2025. Wind power remains the main technology, while solar photovoltaic energy has grown rapidly. The increase in renewable generation has changed the market dynamics [2].

Since 2023, the market has seen more zero and negative price hours—events rare earlier—mainly during times of high renewable output and low demand, especially in spring and early summer. This shift has increased intraday volatility, caused abrupt price jumps, and occasional price collapses, making standard volatility models less effective. There is a body of literature suggesting that a higher share of renewable energy leads to an increase in volatility [28,29,30].

Electricity prices also exhibit distinct patterns, including marked intraday, weekly, and seasonal cycles driven by demand and solar profiles; persistent and high volatility and long memory; and exogenous jumps related to weather, renewable intermittency, and policy actions and nonlinear dynamics in prices, limiting traditional modeling approaches [10,31].

To capture these features, intraday prices are aggregated to produce daily realized volatility, based on hourly returns. This approach keeps the information from high-frequency data while supporting daily risk forecasts. Realized volatility is especially suitable because squared daily returns miss latent volatility due to jumps and intraday effects. The high renewable share also increases jump activity, making jump-robust measures like Median Realized Volatility essential for separating persistent volatility from temporary shocks.

The analysis covers the period from 2021 to 2025. The Figure 1, Figure 2, Figure 3 and Figure 4 illustrate the corresponding price dynamics.

Figure 1 shows hourly prices from 2021 to 2025. Prices were stable early in 2021, between 40 and 100 €/MWh. After mid-2021, both levels and variability sharply increased, with spikes over 300–400 €/MWh and occasional drops near 0 €/MWh, indicating increased volatility. In Spain, although renewable penetration mitigated part of the upward pressure, electricity prices remained closely linked to gas-fired generation through the marginal pricing mechanism [32].

In 2022, Figure 2, prices spiked due to the outbreak of the Russia–Ukraine conflict, which triggered a severe supply shock in European gas markets. The sharp reduction in Russian gas flows led to unprecedented increases in gas prices, directly translating into extreme electricity price spikes across EU power markets [3,32]. In OMIE, this period is characterized by abrupt price surges and exceptionally high volatility, with peak prices reflecting both elevated gas marginal costs and heightened risk premia

Figure 3 depicts a partial normalisation in 2023, with prices trending downward. In March and April, prices frequently neared 0 €, driven by more renewables and the gas cap introduced by the Spanish government [33,34]. The moderation observed has also been documented in other European Markets during the same period [35].

This pattern became more pronounced in 2024 (Figure 4), with zero-price episodes increasing, especially in spring. Although prices recovered drops to zero remained common.

Figure 4 is the hourly prices of 2024, where the frequency of zero-price episodes rose significantly, particularly in spring. Despite some price recovery in the second half of the year, the market continued to exhibit recurrent drops to zero.

By 2025, Figure 5, the downtrend continued, with zero-price hours increasing and negative prices appearing more frequently. Prices stabilized around 100 €/MWh toward year’s end, but high variance and recurring zeros persisted. High volatility during the 2024–2025 could be attributed to geopolitical tensions possibly disruptions of the production and growing share of renewable energy production [9,35,36].

Overall, the Spanish electricity market from 2021 to 2025 presents a highly complex environment for volatility modeling due to extreme price swings, renewable growth, regulatory actions, and negative prices, offering a rigorous test for forecasting methods amid elevated uncertainty. To better understand the evolution of these market dynamics, Figure 6 illustrates two defining patters in price behavior.

Figure 6 shows average hourly spot electricity prices (€/MWh) by hour of the day for the period 2021–2025, highlighting a marked transformation in intraday price patterns in the Spanish electricity market. A clear “duck curve” shape emerges in the most recent years, particularly from 2024 onward.

During 2021, the intraday profile is relatively flat, with moderate night–day price differentials and limited midday depression, reflecting pre crisis conditions and lower renewable penetration. In 2022, prices increase sharply across all hours, and the duck curve structure is largely obscured by crisis driven price levels. Elevated gas prices during the peak of the European energy crisis dominate price formation, compressing intraday variability upward.

From 2023 onward, and especially in 2024–2025, the intraday pattern becomes increasingly pronounced. Prices fall sharply during midday hours and recover rapidly toward the evening peak. This corresponds to the canonical duck curve, where high solar generation depresses midday prices (the “belly” of the duck), followed by steep price ramps in the late afternoon and evening as solar output declines and demand rises. The depth of the midday trough and the steepness of the evening ramp intensify over time, indicating a growing influence of variable renewable generation on price formation.

Overall, Figure 6 illustrates a shift from crisis-dominated price levels in 2022 to a post-crisis regime characterised by structural intraday heterogeneity driven by renewable intermittency. While average price levels moderate after the energy crisis, the amplification of the duck curve pattern implies greater intraday volatility and sharper ramps, reinforcing the need for volatility models capable of capturing time of day effects, jumps, and persistent heteroskedasticity in modern electricity markets.

3. Volatility Modeling Methodology

This section describes the methodology used to model and forecast volatility in the Spanish electricity market. The proposed procedure is organized around four key components. First, classic volatility models are introduced as benchmark specifications. Second, realized volatility measures constructed from intraday data are defined to approximate latent volatility. Third, specific data pre-processing steps are applied to address features unique to electricity prices, such as zero or negative values and strong intraday seasonality. Finally, a rich class of HAR-type realized-volatility models is presented, incorporating heterogeneous persistence, jumps, robust estimators, functional transformations, and calendar effects.

3.1. GARCH Models

Autoregressive Conditional Heteroskedasticy (ARCH) [12] and Generalized Autoregressive Conditional Heteroskedasticty (GARCH) [13] volatility models constitute a common benchmark for modeling and forecasting conditional variance in financial and commodity markets. These models treat volatility as a latent process driven by past returns and innovations, without explicitly exploiting intraday price information. The seminal contribution of Engle [12] introduced the Autoregressive Conditional Heteroskedasticity (ARCH), later generalized by Bollerslev [13] through the GARCH model to allow for a more parsimonious representation of volatility persistence.

We consider the daily log-returns defined as $r_t=\mathrm{l}\mathrm{o}\mathrm{g}\left({\displaystyle \frac{P_t}{P_{t-1}}}\right)$. We assume the mean follows an ARMA(0,0) process, it remains constant over time (${\mu }_t$). This differs from an ARMA-GARCH model, where the conditional mean varies dynamically. In this paper, we focus solely on volatility dynamics, considering the mean as fixed intercept:

$${\mathrm{r}}_{\mathrm{t}} ={ \mu }_{\mathrm{t}} +{ \mathrm{e}}_{\mathrm{t}}$$

(1)

where $\mu$ is a constant and the innovation $e_t$ follows

$$e_t=\sqrt{h_t} z_t,$$

(2)

with $z_t$ being an i.i.d. sequence with zero mean and unit variance, typically assumed to follow a Gaussian or Student-$t$ distribution.

3.1.1. Classic GARCH Model

The conditional variance ${\sigma }_t^2=\mathrm{V}\mathrm{a}\mathrm{r}(e_t\mid {\mathcal{F}}_{t-1})$, given the information set ${\mathcal{F}}_{t-1}$, evolves according to:

$${\sigma }_t^2={\alpha }_0+{\displaystyle \sum _{i=1}^s}{\alpha }_i{e}_{\left(t-i\right)}^2+\sum _{\left[j=1\right]}^r{\beta }_j{\sigma }_{\left(t-j\right)}^2$$

(3)

${\alpha }_0>0$: It constitutes the constant term, symbolizing the baseline or average level of long-term volatility.

${\alpha }_i\ge 0$: They are the ARCH coefficients, which model the immediate effect of the residuals.

${\beta }_j\ge 0$: These are the GARCH coefficients, which incorporate the moving average component for the conditional variance.

$\left({\sum }_{i=1}^s\left({\alpha }_i\right)+{\sum }_{j=1}^r\left(B_j\right)\right)<1$ This is the variance persistence. If it equals 1, it indicates infinite persistence in variance, suggesting an IGARCH [37] model should be used. However, if the values are close to 1, it could also imply a change in the mean value, requiring the use of different models [38].

Equation (2) presents the conditional variance. However, [7] explains, the classical GARCH model assumes symmetric effects of shocks, which is unrealistic because empirically most cases involve asymmetric effects.

3.1.2. Absolute Value GARCH (AVGARCH)

Absolute Value GARCH, also called Taylor-Schwert GARCH proposed by [15,16] is a variant of the GARCH model defined as:

$${\sigma }_t={\alpha }_0+{\displaystyle \sum _{i=1}^p}{\alpha }_i\left|\left(e_{\mathrm{t}-\mathrm{i}}\right)\right|+{\displaystyle \sum _{j=1}^q}{\beta }_j{\sigma }_{\mathrm{t}-\mathrm{j}}$$

(4)

This GARCH model characterizes the conditional standard deviation utilizing the absolute value of residuals. The principal objective of this model is to reduce the influence of extreme or atypical observations [39].

Using the same restrictions as the classical GARCH model.

${\alpha }_0>0$: It constitutes the constant term, symbolizing the baseline or average level of long-term volatility.

${\alpha }_i\ge 0$: They are the ARCH coefficients, which model the immediate effect of the residuals.

${\beta }_j\ge 0$: These are the GARCH coefficients, which incorporate the moving average component for the conditional variance.

3.1.3. GARCH-X

The GARCH-X model is an extension of the classical GARCH that incorporates exogenous variables (X_t) into the variance equations as follows. While standard GARCH relies solely on internal history, GARCH-X allows external factors to influence volatility predictions. So the GARCH-X (s,r,k) is defined as:

$${\sigma }_t^2={\alpha }_0+{\displaystyle \sum _{i=1}^s}{\alpha }_i{e}_{\left(t-i\right)}^2+{\displaystyle \sum _{\left[j=1\right]}^r}{\beta }_j{\sigma }_{\left(t-j\right)}^2+{\displaystyle \sum _{p=1}^k}\gamma X_{t-p}$$

(5)

${\alpha }_0>0$: It constitutes the constant term, symbolizing the baseline or average level of long-term volatility.

${\alpha }_i\ge 0$: They are the ARCH coefficients, which model the immediate effect of the residuals.

${\beta }_j\ge 0$: These are the GARCH coefficients, which incorporate the moving average component for the conditional variance.

$\gamma \ge 0$ it captures the impact of exogenous variables.

While the GARCH-X can accommodate a vector of k covariates, this study focuses on a single exogenous variable: Realized Volatility (RV). This choice is driven by the ‘Information Hypothesis,’ as RV leverages high-frequency intra-day data to offer a more accurate and observable measure of price variation, unlike squared daily returns, which are noisier proxies for volatility.

GARCH models have been extensively applied in energy and electricity markets to capture time-varying volatility and volatility clustering (see, among others [30,37]). However, empirical evidence suggests that their forecasting performance is often constrained by the exclusive reliance on daily returns, which discard a substantial amount of information contained in intraday price movements. This limitation is particularly relevant in electricity markets, where prices are characterized by strong intraday seasonality, frequent jumps, and extreme observations.

These shortcomings are especially pronounced in the Spanish electricity market during the 2021–2025 period, marked by high renewable penetration, increased jump activity, and the emergence of zero and negative prices. As a result, daily return-based volatility estimates may provide a noisy and incomplete representation of the underlying volatility process. For this reason, GARCH-type models are employed in this study primarily as benchmark specifications, against which the performance of realized-volatility-based and HAR-type models is evaluated within a unified out-of-sample comparison procedure.

While return-based GARCH models rely exclusively on daily price changes and treat volatility as a latent process, realized-volatility-based approaches explicitly exploit intraday information, allowing volatility to be measured rather than inferred. This informational advantage motivates the use of HAR-type models built on realized volatility measures, which are introduced in the next section.

3.2. Realized Volatility Measures

Intraday data availability across financial markets has grown, exposing volatility patterns that classical models like GARCH or stochastic volatility often miss, especially when trying to capture the persistent and long-memory traits seen in realized volatility [17]. To overcome this, the Heterogeneous Autoregressive model of Realized Volatility (HAR-RV) was introduced by [21].

$$d{p}_t=u_{t}dt+{\sigma }_{t}d{W}_t+k_{t}d{q}_t$$

(6)

Let the logarithmic price of an asset follow a jump-diffusion process over time, with the dynamics driven by both a continuous diffusion component and a discrete jump component. In this model, the drift term is represented by u_t. ${\sigma }_t$ is a stochastic volatility process independent of the $W_t$ standard Brownian motion, $k_t$ and denotes the jump size in log prices.

Realized Volatility exhibits some unique traits, including notable autocorrelations at long lags, indicating long-range dependence. It also shows a leverage effect, where volatility bursts are more likely following negative returns and jumps [3].

The volatility of the price process at day t is measured by quadratic variation ($\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{g}\mathrm{r}\mathrm{a}\mathrm{t}\mathrm{e}\mathrm{d} \mathrm{v}\mathrm{a}\mathrm{r}\mathrm{i}\mathrm{a}\mathrm{n}\mathrm{c}\mathrm{e} \mathrm{a}\mathrm{s}\mathrm{s}\mathrm{o}\mathrm{c}\mathrm{i}\mathrm{a}\mathrm{t}\mathrm{e}\mathrm{d}, {IV}_t)$:

$${IV}_t={\int }_t^{t+1}{\sigma }_s^{2}ds+\sum _{t<s<t+1}{k}_s^2$$

(7)

Being ${\sigma }_s^2$ the contribution of the continuous sample path. And $\sum _{t<s<t+1}{k}_s^2$ is the contributions of the jumps.

However, it is not directly observable, so an estimator is needed [41]. If there were no noise due to microstructure, the Realized Variance would be the best estimator. In the presence of noise, several robust estimators are available. Therefore, in this work, we have estimated the classic Realized Variance and also used robust estimators, as described below.

3.2.1. Realized Volatility

Realized volatility is then as follows:

$$R{V}_t={\displaystyle \sum _{j=1}^M}{r}_{t,j}^2$$

(8)

As the sampling frequency increases, $R{V}_t$ converges in probability to the integrated variance of the continuous-time price process [17]. Empirically, RV has been shown to substantially outperform squared daily returns as a proxy for latent volatility. However, RV is sensitive to jumps and microstructure noise [18,20], which may dominate volatility dynamics in electricity markets characterized by abrupt price movements and extreme observations.

3.2.2. Noise-Robust Realized Volatility Measures

The Kernel Realized Variance, proposed in [19], is defined as:

$$R{K}_t={\gamma }_{0,t}+{\sum }_{h=1}^{H}K({\displaystyle \frac{h-1}{H}})({\gamma }_{h,t}+{\gamma }_{-h,t})$$

(9)

where ${\gamma }_{h,t}$ denotes intraday return autocovariances of order $h$ on day t:

$${\gamma }_{h,t}= {\displaystyle \sum _{j=h+1}^M}{r}_{t,j}\ast r_{t,j-h}$$

(10)

$K(\cdot )$ is a kernel weighting function, and $H$ is a bandwidth parameter controlling the number of autocovariances included in the estimation. By aggregating both positive and negative return autocovariances through a spectral weighting scheme, the Kernel Realized Variance estimator corrects for the impact of microstructure noise contaminating raw high-frequency returns.

The Kernel Realized Variance estimator provides a robust benchmark. Its inclusion allows us to assess whether volatility dynamics and forecasting performance are sensitive to alternative noise-robust volatility proxies. In this work, we only utilize the Bartlett Kernel, defined by K(x) = 1 − x for $0\le x\le 1$.

3.2.3. Jump-Robust Realized Volatility Measures

To address the sensitivity of standard Realized Variance to discontinuities in the price process, the literature has developed alternative estimators that are robust to jumps. Among these, Median Realized Volatility (MedRV) is used due to its robustness and ease of implementation [20].

Median Realized Volatility is a jump-robust estimator defined as

$${MedRV}_t={\displaystyle \frac{\pi }{6-4\sqrt{3}+\pi }}{\displaystyle \frac{M}{M-2}}{\displaystyle \sum _{i=2}^{M-1}}[med(\mid r_{t,i-1}\mid ,\mid r_{t,i}\mid ,\mid r_{t,i+1}\mid ){]}^2$$

(11)

where $\mathrm{m}\mathrm{e}\mathrm{d}(\cdot )$ denotes the median operator. [20] show that MedRV is consistent for the integrated variance under a wide class of jump-diffusion processes and effectively filters extreme intraday returns associated with discontinuous price changes.

By construction, MedRV isolates the continuous component of quadratic variation, thereby providing a more reliable measure of persistent volatility dynamics in the presence of jumps. This property is particularly important in electricity markets, where price formation is frequently affected by abrupt supply-demand imbalances, renewable intermittency, and system constraints.

In the Spanish electricity market over the 2021–2025 period, high renewable penetration has been associated with an increased incidence of price spikes and outliers. In such an environment, jump-robust estimators such as MedRV offer a more accurate representation of the underlying continuous volatility component that drives persistence and predictability.

3.3. Data Pre-Processing: Price Shift and Deaming Returns

Electricity markets occasionally exhibit zero or negative prices, which makes impossible todo the logarithm. We apply a price displacement (shift) to ensure that all values are positive

$$P_{t,j}^{shift}=P_{t,j}+\lambda$$

(12)

where $\lambda =1-\mathrm{m}\mathrm{i}\mathrm{n}\left(P\right)$ ensures that the minimum price is equal to 1. This linear transformation allows applying the log-return calculations, While the price displacement ensures that $P_{t,j}^{shift}>0$, allowing for the application of logarithmic return calculations. It is important to acknowledge that this linear transformation could introduce a specific degree of measurement bias. By shifting the price level, the returns are altered, particularly during periods of low absolute prices where $\lambda$ the constant represents a larger proportion of the total value. However, in the context of the Spanish Electricity Market, this bias is considered a necessary trade-off to address the numerical instability caused by zero or negative prices, which would otherwise lead to undefined log-returns and extreme values. To get the daily prices, we calculated the average for each day.

In addition, predictable intraday and seasonal patterns driven by demand cycles and renewable production may distort volatility estimates. To isolate stochastic price movements, intraday returns are demeaned following [10].

In the work of [10] it is used to demean the month of the year, day of the week, and half-hour of the day. In the MIBEL market, spot prices are recorded hourly until October 2025. Therefore, instead of using half-hour periods, the data is adjusted using hourly values [42].

$$\widehat{{\mu }_{t,j}}=med\left(r_{mm,dy,hr}\right)$$

(13)

We calculate the median return for the day t in month mm, on day of the week dy and hour j.

Then the demeaning returns are:

$$r_{t,j}^{*}=r_{t,j}-\widehat{{\mu }_{t,j}}$$

(14)

All realized volatility measures are computed using the adjusted returns $r_{t,j}^{*}$, ensuring that volatility is constructed on a stable and economically meaningful return series and reflects genuine uncertainty rather than deterministic price patterns. These preprocessed volatility measures are used as inputs to HAR-type models capturing heterogeneous volatility persistence across multiple time horizons.

3.4. HAR Volatility Models

To model the realized volatility, we employ Heterogenous Autoregressive (HAR) specifications, which characterize volatility persistence through linear combinations of realized volatility components aggregated over multiple time horizons. Introduced by Corsi [21], the Heterogeneous Autoregressive (HAR) model provides a parsimonious linear approximation to long-memory behavior by combining daily, weekly, and monthly realized volatility measures, without imposing explicit fractional integration.

Let $R{V}_t$ denote a realized volatility measure for day $t$, constructed from preprocessed intraday returns as described in Section 3.2. The baseline HAR-RV model serves as the reference specification and captures heterogeneous volatility persistence across short-, medium-, and long-run horizons.

Building on the general jump-diffusion representation of the price process, the HAR model can be extended to explicitly account for discontinuous price movements, giving rise to HAR-J models that separate persistent volatility dynamics from transitory jump-induced variation. In addition, since realized volatility measures are subject to time-varying measurement error, the HAR specification is further augmented using realized quarticity, leading to HAR-Q models. These extensions are jointly combined in the HAR-QJ specification, which simultaneously captures heterogeneous persistence, jump dynamics, and measurement uncertainty.

Finally, electricity markets exhibit pronounced calendar effects, particularly day-of-the-week patterns in volatility arising from systematic demand cycles, generation constraints, and institutional market arrangements. To account for residual seasonal heteroskedasticity, HAR-type models are further augmented with weekday dummy variables, allowing calendar effects to be modeled jointly with volatility persistence, jump behavior, and time-varying measurement precision within a unified and internally consistent framework.

3.4.1. HAR-RV Model

The baseline HAR-RV model is specified as

$${RV}_t={\beta }_0+{\beta }_1{RV}_{t-1}^d+{\beta }_2{RV}_{t-1}^w+{\beta }_3{RV}_{t-1}^m+{\epsilon }_t$$

(15)

where $R{V}_t^d=R{V}_{t-1}$ denotes daily realized volatility. The weekly and monthly components are defined as

$${RV}_t^w={\displaystyle \frac{1}{7}}{\displaystyle \sum _{i=0}^6}{RV}_{t-i}$$

(16)

$${RV}_t^m={\displaystyle \frac{1}{30}}{\displaystyle \sum _{i=0}^{29}}{RV}_{t-i}$$

(17)

Corresponding to 7-day and 30-day moving averages of realized volatility.

The coefficients ${\beta }_1$, ${\beta }_2$, and ${\beta }_3$ measure volatility persistence across short-, medium-, and long-run horizons, respectively. By construction, the HAR model provides a parsimonious linear approximation to long-memory behavior by combining realized volatility components observed at multiple aggregation scales, without imposing explicit fractional integration. This representation is well suited to environments characterized by persistent volatility dynamics and recurring regime changes.

3.4.2. HAR-Q Model

Let $I{V}_t$ denote the latent daily volatility given by quadratic variation (Equation (7)). The realized variance $R{V}_t$, constructed from intraday returns, provides a consistent but noisy estimator of $I{V}_t$, which can be expressed as:

$$R{V}_t=I{V}_t+{\eta }_t$$

(18)

where ${\eta }_t$ denotes the measurement error. In finite samples, this error arises from discretization, heavy-tailed intraday returns, and extreme price movements, and its variance is time-varying, being larger during periods of intense intraday activity.

The conditional variability of the measurement error depends on the fourth moment of intraday returns. This motivates the use of Realized Quarticity (RQ) as a proxy for time-varying measurement uncertainty, defined as:

$${RQ}_t= {\displaystyle \frac{M}{3}}{\displaystyle \sum _{i=1}^M} r_{t,i}^4$$

(19)

where $r_{t,i}$ denotes intraday returns and $M$ is the number of intraday observations.

The HAR-Q model [43] incorporates realized quarticity to adjust the predictive contribution of realized volatility components according to their measurement precision. It is specified as:

$${RV}_t= {\beta }_0+({\beta }_1+{\beta }_{1Q}{RQ}_{t-1}^{(d){\displaystyle \frac{1}{2}}}){RV}_{t-1}^d+ {(\beta }_2{+\beta }_{2Q}{RQ}_{t-1}^{(w){\displaystyle \frac{1}{2}}}){RV}_t^w+({\beta }_3{+\beta }_{3Q}{RQ}_{t-1}^{(m){\displaystyle \frac{1}{2}}}){RV}_{t-1}^m+ {\epsilon }_t$$

(20)

where $R{Q}_{t-1}^{\left(d\right)}$, $R{Q}_{t-1}^{\left(w\right)}$, and $R{Q}_{t-1}^{\left(m\right)}$ denote daily, weekly, and monthly aggregated quarticity measures. High quarticity values correspond to low-precision volatility estimates and therefore lead to an endogenous down-weighting of realized volatility components, improving estimation efficiency and forecasting accuracy [41,42].

3.4.3. HAR-J Model

Electricity prices are characterized by frequent discontinuous movements, renewable intermittency, weather shocks, and regulatory or market interventions [2,28,36]. From a continuous-time perspective, these features imply that the price process cannot be adequately described by a pure diffusion and must allow for jumps. Consequently, the total daily quadratic variation of prices contains both a continuous component and a jump component:

$$R{V}_t=C_t+J_t$$

(21)

Continuous component ($C_t)$: Let $R{V}_t$ denote the realized variance obtained from intraday returns. In the presence of jumps, $R{V}_t$ captures both continuous and discontinuous variation. To isolate the continuous component of quadratic variation, we use the Bipower Variation (BPV) estimator proposed by [18], defined as:

$${BPV}_t= {\displaystyle \frac{\pi }{2}}\left({\displaystyle \frac{M}{M-i}}\right){\displaystyle \sum _{i=1}^{M-i}} \left|r_{t,i}\right|\left|r_{t,i+1}\right|$$

(22)

where $r_{t,i}$ are intraday returns, $M$ is the number of intraday observations, and $i$ is a lag parameter.

Under standard regularity conditions, $BP{V}_t$ is a consistent estimator of the integrated variance associated with the continuous diffusion component, even in the presence of jumps. Intuitively, the product of adjacent absolute returns eliminates the dominant impact of isolated large jumps, which tend to inflate squared returns disproportionately

Jump identification: Although both $R{V}_t$ and $BP{V}_t$ are observable, a statistical procedure is required to determine whether their difference reflects jump activity or sampling variability. This is achieved through the $Z_t$ proposed by [46]:

$$Z_t= {\displaystyle \frac{{(RV}_t-{BPV}_t{)/RV}_t}{\sqrt{\left({\left(\frac{\pi }{2}\right)}^2+\pi -5\right)\frac{1}{M}max\left(1,\frac{{TPQ}_t}{{BPV}_t^2}\right)}}}$$

(23)

Being TPQ the tripower quarticity for day t:

$${TPQ}_t= {\displaystyle \frac{1}{4}}{\left[{\displaystyle \frac{\mathsf{\Gamma }(7/6)}{\mathsf{\Gamma }(1/2)}}\right]}^{-3}{\displaystyle \frac{M^2}{M-2i}}{\displaystyle \sum _{j=1}^{M-2i}}{\left|r_{t,j}\right|}^{4/3}{\left|r_{t,j+i}\right|}^{4/3}{\left|r_{t,j+2i}\right|}^{4/3}$$

(24)

Which provides a robust estimator of the integrated fourth moment of intraday returns. Tripower Quarticity accounts for time-varying volatility and heavy-tailed return distributions and is used to consistently estimate the asymptotic variance of the bipower variation. Here, $\mathsf{\Gamma }(\cdot )$ denotes the Gamma function.

Under the null hypothesis of no jumps and as $M\to \mathrm{\infty }$, the standardized test statistic satisfies [44]

$$Z_t\stackrel{d}{\to }\mathcal{N}\left(\mathrm{0,1}\right),$$

(25)

Significant jumps are identified when the $Z_t$ statistic exceeds a standard normal critical value ${\Phi }_{\alpha }$ [44]. In this study, following the recommendations of [10]. In our empirical pipeline, we set alpha = 0.01, representing a 99% confidence level, with a critical value of ${\Phi }_{0.01}$ = 2.326. Consequently, the daily jump component is defined as:

$$J_t= I(Z_t>{\Phi }_{\alpha })({RV}_t-{BPV}_t)$$

(26)

where $I(\cdot )$ denotes the indicator function. This construction ensures that only statistically significant deviations between realized variance and bipower variation are attributed to jumps, while small discrepancies are treated as estimation noise.

Following [10] we increase the lag (i = 3) for the Bipower variation and Tripower quarticity to improve the performance of the test. In our case we use the Jumps for 1 day.

The HAR-J model is defined as:

$${RV}_t= {\beta }_0+ {\beta }_1{RV}_{t-1}^d+ {\beta }_2{RV}_{t-1}^w+ {\beta }_3{RV}_{t-1}^m+{\beta }_4{J}_{t-1}+ {\epsilon }_t$$

(27)

3.4.4. HAR-QJ Model

The HAR-QJ model combines the two previous extensions by jointly accounting for measurement error and jump dynamics. This specification is particularly relevant in electricity markets, where realized volatility is affected simultaneously by intraday noise and frequent discontinuities.

The HAR-QJ model jointly accounts for time-varying measurement error and jump dynamics

$$\begin{array}{cc}R{V}_t=& {\beta }_0+({\beta }_1+{\beta }_{1Q}\sqrt{R{Q}_{t-1}^{\left(d\right)}})R{V}_{t-1}^d+({\beta }_2+{\beta }_{2Q}\sqrt{R{Q}_{t-1}^{\left(w\right)}})R{V}_{t-1}^w\\ & +({\beta }_3+{\beta }_{3Q}\sqrt{R{Q}_{t-1}^{\left(m\right)}})R{V}_{t-1}^m+{\beta }_4{J}_{t-1}+{\epsilon }_t.\end{array}$$

(28)

By construction, HAR-QJ allows volatility persistence to be estimated conditionally on both time-varying measurement precision and jump activity, yielding a more robust representation of realized volatility dynamics.

3.4.5. HAR-C and HAR-CJ Models

Following [11], we decompose the realized variance into its continuous and jump components to assess their individual predictive powers. We define the jump-adjusted continuous component C_t as RV_t when no significant jump is detected, and as BPV otherwise. This leads to the HAR-C (Continuous) and HAR-CJ (Continuous-Jump) models:

$$C_t=\{\begin{array}{cc}R{V}_t,& \mathit{if} \mathit{no} \mathit{jump} \mathit{is} \mathit{detected},\\ BP{V}_t,& \mathrm{otherwise}.\end{array}$$

(29)

This leads to HAR-C model,

$${RV}_t={\beta }_0+{\beta }_1{C}_{t-1}^d+{\beta }_2{C}_{t-1}^w+{\beta }_3{C}_{t-1}^m+{\epsilon }_t$$

(30)

And the HAR-CJ model,

$${RV}_t={\beta }_0+{\beta }_1{C}_{t-1}^d+{\beta }_2{C}_{t-1}^w+{\beta }_3{C}_{t-1}^m+{\beta }_4{J}_{t-1}+{\epsilon }_t$$

(31)

3.4.6. HAR Median and HAR Kernel Model and Transformations

To assess the role of robustness in volatility forecasting, the HAR model is further extended by replacing standard realized variance with jump-robust and noise-robust volatility estimators.

$$Es{t}_t\in \{{\mathit{MedRV}}_t, {\mathit{RK}}_t\}$$

(32)

where MedRV denotes Median Realized Volatility and RK, Kernel Realized Volatility.

HAR-M (Median) and HAR-K (Kernel):

$${RV}_t= {\beta }_0+ {\beta }_1{Est}_{t-1}^d+ {\beta }_2{Est}_{t-1}^w+ {\beta }_3{Est}_{t-1}^m+ {\epsilon }_t$$

(33)

where $Es{t}_t^d$, $Es{t}_t^w$, and $Es{t}_t^m$ denote daily, weekly, and monthly aggregated components of the alternative realized volatility estimator

HAR-MJ and HAR-KJ:

$${RV}_t= {\beta }_0+ {\beta }_1{Est}_{t-1}^d+ {\beta }_2{Est}_{t-1}^w+ {\beta }_3{Est}_{t-1}^m+ {\beta }_4{J}_{t-1}+{\epsilon }_t$$

(34)

where $J_t$ denotes the jump component identified using bipower variation and jump tests.

Given the heavy-tailed nature of electricity price volatility and the frequent occurrence of extreme observations, all HAR-type models are estimated under alternative functional transformations of realized volatility. Specifically, model estimation is performed using $R{V}_t$, $\mathrm{l}\mathrm{o}\mathrm{g}\left(R{V}_t\right)$, and $\sqrt{R{V}_t}$, corresponding to level, logarithmic, and square-root transformations, respectively. These transformations are employed to stabilize variance, reduce the influence of outliers, and enhance numerical stability and out-of-sample forecasting performance.

3.4.7. HAR Models with Day-of-the-Week Effects

Electricity markets are heavily influenced by the day of the week, as both consumption and prices exhibit a pronounced weekly seasonal pattern. Following the approach in [22], we incorporate day-of-the-week dummy variables to examine their specific impact on realized variance. This is a crucial step because, while the demeaning of returns (Equations (13) and (14)) accounts for deterministic patterns in the mean price level, it does not address the seasonal heteroskedasticity inherent in the variance process.

$${RV}_t=HA{R}_{model}+{\displaystyle \sum _{i=2}^7}{\phi }_i{D}_{i,t}+{\epsilon }_t$$

(35)

where HAR model denotes any of the HAR specifications introduced in Section 3.4.1, Section 3.4.2, Section 3.4.3, Section 3.4.4, Section 3.4.5 and Section 3.4.6.

Let $D_{i,t}$ denote indicator variables defined as:

$$D_{i,t}=\{\begin{array}{cc}1,& \mathrm{if} \mathrm{day} t \mathrm{corresponds} \mathrm{to} \mathrm{weekday} i,\\ 0,& \mathrm{otherwise},\end{array} ; i=2,\dots ,7,$$

(36)

Monday serves as the reference category. Since the electricity market operates continuously (24/7), we include six dummy variables to represent the remaining days of the week. This seasonal adjustment is implemented across all HAR specifications and functional transformations previously described.

The coefficients ϕ measure systematic deviations in realized volatility on weekday I relative to Monday, conditional on the information contained in the HAR dynamics. A statistically significant ϕ indicates the presence of weekly seasonal heteroskedasticity.

Importantly, the inclusion of weekday dummies does not alter the dynamic structure of the HAR model but augments it with deterministic regressors that capture residual calendar regularities in the variance process. This allows volatility persistence, jump dynamics, measurement error, and calendar effects to be modeled jointly within a unified linear framework.

3.5. Out of Sample Evaluation and Statistical Model Selection

The models are estimated by using an expanding window approach to evaluate their predictive accuracy. The initial estimation subsample (In-Sample) spans from 2021-01-02 to 2024-12-31, while the entire year of 2025 is reserved for out-of-sample (OOS) evaluation. One day-ahead forecasts are generated iteratively, with the estimation window expanding after each prediction to incorporate the most recent observation. Model parameters are re-estimated at each step, ensuring that forecasts are based on all available information up to the prediction date.

Figure 7 provides an overview of the out-of-sample evaluation and statistical model selection procedure. Out-of-sample volatility forecasts are evaluated and compared using a multi-stage procedure designed to ensure statistical coherence and robustness. Model selection is based exclusively on a loss-based criterion within the Model Confidence Set (MCS). Specifically, forecasts are evaluated using the QLIKE loss function, and the MCS procedure yields a Superior Set of Models (SSM) satisfying the Equal Predictive Ability hypothesis [25]. Conditional on SSM membership, additional diagnostics—namely the out-of-sample R² and the Mincer–Zarnowitz regression—are employed to assess forecast accuracy and bias. These diagnostics are purely descriptive and do not influence the loss-based model selection.

3.5.1. Loss-Based Model Comparison: MCS with QLIKE

Model Confidence Set (MCS)

The Model Confidence Set proposed by [25] provides a statistically approach for jointly comparing multiple forecasting models and identifying a subset of specifications that are not statistically dominated under a common loss function. Let $\{L_{i,t}{\}}_{t=1}^T$ denote the sequence of out-of-sample losses generated by model $i\in {\mathcal{M}}_0$, where ${\mathcal{M}}_0$ is the initial model set. For each pair of models $\left(i\right|j)$, define the loss differential as:

$$d_{ij,t}=L_{i,t}-L_{j,t}$$

(37)

with expected value $\mathbb{E}\left[d_{ij,t}\right]$. The null hypothesis of Equal Predictive Ability (EPA) for the model set $\mathcal{M}\subseteq {\mathcal{M}}_0$ is given by

$$H_{0,\mathcal{M}}: \mathbb{E}{d}_{ij,t}=0,\forall i,j\in \mathcal{M}.$$

(38)

Rejection of the null hypothesis implies that at least one model in the set exhibits inferior predictive performance. The MCS procedure sequentially eliminates the model associated with the largest standardized average loss until the EPA hypothesis can no longer be rejected.

$$T_{\mathit{MAX}}=\underset{i\in \mathcal{M}}{max}\left({\displaystyle \frac{{\bar{d}}_i}{\sqrt{\widehat{Var}\left({\bar{d}}_i\right)}}}\right),$$

(39)

where

$${\bar{d}}_i={\displaystyle \frac{1}{\mid \mathcal{M}\mid -1}}\sum _{\begin{array}{c}j\in M\\ j\ne i\end{array}}{\bar{d}}_{ij}$$

(40)

Denotes the average differential loss of model $i$ relative to all remaining competitors, and $\widehat{\mathrm{V}\mathrm{a}\mathrm{r}}\left({\bar{d}}_i\right)$ is a consistent estimator of its variance. Statistical inference is conducted using a moving block bootstrap to account for serial dependence in the loss differentials. The output of the procedure is the Superior Set of Models (SSM), defined as a subset of models for which the EPA hypothesis cannot be rejected at the chosen significance level. In this application, the Model Confidence Set enables a statistically coherent comparison of volatility forecasting models estimated under different information sets, functional transformations, and volatility proxies, while avoiding ad-hoc pairwise testing and data-snooping biases.

QLIKE Loss Function

The error function we use is QLIKE, an asymmetrical metric that penalizes underprediction more than overprediction [24,45]. We use realized volatility as a proxy because it is less noisy than squared returns, although QLIKE itself is a robust error measure that handles noisy proxies well [24].

$$\mathrm{Q}\mathrm{L}\mathrm{I}\mathrm{K}\mathrm{E}=\ln \left(h_t\right)+{\displaystyle \frac{{RV}_t}{h_t}}$$

(41)

The choice of QLIKE is motivated by properties in volatility forecasting. First, QLIKE is robust to measurement error in volatility proxies and remains consistent for ranking competing forecasts even when the realized measure is noisy or affected by jumps and microstructure effects [24].

Second, QLIKE is intrinsically asymmetric, imposing a larger penalty on volatility underprediction than on over-prediction. This feature has a clear economic interpretation in electricity markets, where underestimating volatility entails asymmetric costs related to risk management, such as insufficient collateral requirements, mispricing of derivatives, and exposure to extreme price movements. Given the empirical characteristics of electricity prices, frequent jumps, heavy tails, and pronounced volatility clustering, this asymmetry renders QLIKE particularly appropriate.

Importantly, the use of QLIKE ensures coherence between the forecast evaluation metric and the MCS testing, as the loss differentials entering the EPA hypothesis are computed using a strictly consistent volatility scoring rule. Consequently, model selection decisions within the MCS are driven by robust loss differences.

3.5.2. Forecast Evaluation Within the Superior Set of Models

The MCS procedure yields a Superior Set of Models (SSM) containing volatility forecasting specifications that are statistically indistinguishable in terms of expected QLIKE loss. Inclusion in the SSM reflects loss-based equivalence rather than absolute forecast dominance. Consequently, additional performance measures are employed conditionally on SSM membership, with the objective of characterizing forecast accuracy and statistical optimality rather than overriding the loss-based selection mechanism.

R² Out-of-Sample

The MCS process yields a set of models that are predictively the best, although statistically they are indistinguishable. Hence, we use the Out-of-sample R² ($R_{OOS}^2)$ and the Mincer-Zarnowitz regression to identify the best, unbiased model from this set.

R² Out-Of-Sample proposed by [48] is:

$$R_{OOS}^2=1-{\displaystyle \frac{{\displaystyle {\sum }_{t=1}^T}{\left({RV}_t-h_t\right)}^2}{{\displaystyle {\sum }_{t=1}^T}{\left({RV}_t{-\overline{RV}}_{mean,t}\right)}^2}}$$

(42)

where $h_t$ is the prediction of the model for the day t and ${\overline{RV}}_{mean,t}$ is the mean value at the day t. Therefore, with this R², we assess if the prediction outperforms the average value. This measure assesses economic forecast gains relative to a naïve benchmark and is used exclusively to compare models within the SSM, not to determine statistical dominance.

Mincer-Zarnowitz Regression

We use the proposed Mincer-Zarnowitz regression [49]. Using this particular regression we can check if the predictions are unbiased and efficient or not [50]. The Mincer-Zarnowitz regression is defined as:

$${RV}_t= \alpha +{\beta }_1{h}_t+{\epsilon }_t$$

(43)

where $h_t$ is the prediction of the model for the day t.

We do a Joint Hypothesis Test H₀: $\alpha =0$,${\beta }_1=1$: This test determines if the forecast is perfectly unbiased. $\alpha =0$ ensures no systematic constant bias (intercept at the origin ${ \beta }_1=1$. Ensures the forecast reacts correctly to the scale of volatility. As with the $R_{\mathrm{OOS}}^2$, the Mincer–Zarnowitz regression is applied conditionally on SSM membership, ensuring that forecast interpretation remains consistent with the loss-based model selection strategy.

4. Results and Discussion

This section reports the results of the out-of-sample volatility forecasting exercise based on the hierarchical evaluation procedure described in Section 3.5. Model selection is first carried out using the Model Confidence Set (MCS) methodology under the QLIKE loss function, which identifies the Superior Set of Models (SSM) satisfying Equal Predictive Ability. Conditional on SSM membership, forecast performance is further characterized using the out-of-sample coefficient of determination $R_{\mathrm{OOS} }^2$ and Mincer–Zarnowitz regressions, which assess forecast accuracy, unbiasedness, and scale efficiency without affecting the loss-based selection.

Figure 8 illustrates the impact of the price displacement (shift) transformation. As observed, the transformation strictly performs a linear shift of the price levels, while preserving the series’ temporal structure and original dynamics. The daily returns series illustrates a progressive increase in market volatility over time. A surge is observed starting in 2023, with fluctuations becoming particularly pronounced throughout 2024 and 2025. This trend can be attributed to the increased penetration of renewable energy in the total power mix, as well as heightened geopolitical and political uncertainty.

This volatility could be closely linked to the European Union’s (EU) recent efforts to accelerate its decarbonization agenda. By implementing stricter emissions regulations and shifting toward decentralized energy sources, the EU has triggered a structural transformation of the energy market. While these policies aim for long-term sustainability, the transition phase has introduced significant short-term instability and price sensitivity across the sector.

And the realized volatility tells a story of two distinct regimes. Prior to 2022, the series had only occasional spikes. The energy shock of 2022 marks a clear structural break causing more spikes and a higher mean volatility level.

Table 1 shows that although the minimum value changes from −15 to 1, the standard deviation stays the same as the original series. This is important because variance, which is key to our modeling, does not change with linear shifts. Therefore, this adjustment does not create substantial bias in the volatility analysis.

Additionally, alternative approaches like excluding negative observations—which make up 1.84% of our sample—were not used because removing data points could break the time-series continuity and may distort the long-memory features of volatility. Although more complex non-linear transformations are available (such as those suggested by [51]), the selected adjustment preserves interpretability and numerical stability, ensuring that the estimated realized variance accurately captures the market’s original fluctuations without losing information. Applying the price displacement (shift) removes negative and near-zero values, allowing for stable calculation of logarithmic returns.

Table 2 shows descriptive statistics for two series: Raw, which includes intraday log-returns with shifted prices, and Ullrich, which uses the demeaning method outlined in Equation (11).

Consistent with existing research on high-frequency energy data, neither series follows a normal distribution, with high kurtosis and notable skewness. When comparing them, the Ullrich adjustment slightly raises the mean, while the median remains exactly zero. This is an important result, as it shows that the demeaning process effectively centered the returns, isolating idiosyncratic price changes from hourly and seasonal patterns.

Additionally, the Ullrich-adjusted shows a marked decrease in standard deviation, as evidenced by the narrower gap between its minimum and maximum values compared to the Raw series. This lowered dispersion indicates that much of the ‘raw’ volatility was due to predictable seasonal patterns. By eliminating these deterministic factors, the resulting realized volatility better reflects the actual stochastic uncertainty in the MIBEL market.

Table 3 presents the descriptive statistics for the daily realized variance estimators Realized Variance (RV), Median Realized Variance (MedRV), and Kernel Realized Variance, using Bartlett estimator (KRV).

The mean and median values differ considerably across the three measures. The MedRV exhibits the lowest central tendency, which can be attributed to its robustness against outliers and jumps. While the standard RV is inflated by price spikes (jumps), the Median transformation “smooths” the series, providing a more stable estimate of the continuous volatility component.

In contrast, the KRV using the Bartlett kernel shows the highest mean and median values. This indicates that the kernel does not reduce central tendency; instead, it captures a higher level of total variance by accounting for the microstructure noise and autocovariances inherent in hourly price data.

Regarding dispersion, the MedRV successfully reduces the standard deviation, whereas the KRV shows an increase, indicating higher variability in its daily estimates. However, a significant improvement is observed in the “shape” of the distributions:

In terms of skewness, the MedRV shows only a slight improvement over the raw series. However, the KRV significantly reduces skewness, approaching a more symmetric distribution.

The Kernel Realized Volatility kurtosis achieves the lowest excess kurtosis among all estimators, moving closer to a more mesokurtic distribution.

In summary, we observe a technical trade-off. The Kernel estimator shows a “worsening” in terms of central tendency and dispersion (higher and more volatile values), but a clear “improvement” in its distributional shape, producing a series that behaves more like a Normal distribution. On the other hand, the Median Realized Volatility acts as a global smoother, reducing noise and the impact of jumps across all statistical fields.

Figure 9 shows the three estimators over the full sample period. The top panel displays relatively stable and low values during 2021, with isolated spikes appearing in 2022. From 2023 onwards, both the frequency and magnitude of these spikes increase, reaching the highest values in the sample around early 2025.

The middle panel follows the same general pattern but with noticeably smoother fluctuations. The extreme values visible in the RV series are attenuated here, which reflects the median estimator’s ability to reduce the influence of outlying observations. This makes it easier to identify the underlying trend in volatility over time.

The bottom panel tends to record higher values than both RV and MedRV throughout the sample, and its spikes are more pronounced. This is consistent with the descriptive statistics reported earlier, where KRV showed the highest mean and maximum values across the sample. The difference in magnitude relative to the other two estimators becomes more evident from 2023 onward, when market conditions appear to have become more volatile.

The divergence across estimators is worth noting. While all three capture the same broad volatility dynamics, their sensitivity to extreme observations differs considerably. Realized volatility and Kernel are more reactive to large individual returns, whereas Median provides a more stable signal. In practice, this suggests that the choice of estimator may have non-trivial implications for downstream applications such as volatility forecasting or risk measurement, particularly during turbulent market periods.

Table 4 reports the Superior Set of Models (SSM) identified by the Model Confidence Set procedure based on the QLIKE loss function. All specifications retained in the SSM are statistically indistinguishable in terms of expected loss and therefore satisfy the Equal Predictive Ability criterion. Consequently, differences in subsequent performance measures should be interpreted conditionally within this statistically equivalent set, rather than as evidence of loss-based dominance.

Within the SSM, Out-of-Sample $R_{\mathrm{OOS}}^2$ values are used to quantify the magnitude of forecast gains relative to the historical mean benchmark, while Mincer–Zarnowitz regressions assess forecast unbiasedness and scale efficiency. These diagnostics serve to characterize the economic relevance and statistical optimality of the competing forecasts without affecting the loss-based selection. Among the models belonging to the SSM, the HAR-MJ-D specification combines the highest $R_{\mathrm{OOS} }^2$ with statistically unbiased and efficient forecasts and is therefore retained as the preferred model for subsequent analysis.

Several empirical regularities emerge from the loss-based selection results reported in Table 4. First, all GARCH-type specifications are excluded from the Superior Set of Models, indicating systematically inferior predictive performance relative to HAR-type formulations when evaluated under the QLIKE loss. This result is consistent with the inability of classical conditional variance models to fully exploit the information contained in high-frequency price data, particularly in electricity markets characterized by pronounced long-memory dynamics, jump behavior, and strong calendar effects.

Figure 10 provides a comparative analysis between two GARCH specifications, the optimal HAR model (HAR-M-J with day-of-the-week), and the empirical realized volatility. The results indicate that the GARCH models consistently underestimate volatility throughout most of the sample period. While the AVGARCH model shows a slight improvement over the sGARCH in terms of bias, both models remain significantly decoupled from actual volatility levels. Consequently, they are heavily penalized for this systematic underestimation and are subsequently excluded from the Superior Set of Models (SSM).

Second, the Superior Set is dominated by HAR-type specifications estimated on realized volatility measures, with no logarithmic transformations surviving the MCS procedure. This finding suggests that, once deterministic intraday and seasonal patterns are removed via the Ullrich demeaning adjustment, additional variance-stabilizing transformations provide no further gains in out-of-sample performance. In contrast, both the inclusion of jump components and day-of-the-week dummy variables remains highly relevant, underscoring the importance of capturing discontinuous price movements and residual weekly seasonality in volatility forecasting for the Spanish electricity market.

Within the SSM, the HAR-MJ-D specification—based on Median Realized Volatility, an explicit jump component, and weekday indicators—exhibits the highest Out-of-Sample $R_{\mathrm{OOS}}^2$ while simultaneously satisfying the Mincer–Zarnowitz unbiasedness and efficiency conditions. This combination indicates that, among statistically equivalent competitors, the HAR-MJ-D model delivers the most economically meaningful forecast gains relative to a historical mean benchmark without exhibiting systematic bias or scale distortions. These results are consistent with [40], confirming that Median Realized volatility yields the highest forecasting accuracy among the estimators considered.

In electricity markets, Out-of-Sample $R_{\mathrm{OOS}}^2$ values are typically moderate in absolute magnitude. This reflects the intrinsic characteristics of power prices, which differ substantially from standard financial assets. Electricity prices exhibit frequent and sudden jumps, strong non-storability constraints, pronounced weather-driven shocks, regulatory interventions, and structural changes associated with renewable integration. These features generate a large unpredictable component in realized volatility that cannot be fully anticipated using past information alone.

To further examine the sources of this superior performance and assess the economic relevance of the model components, the following Table 5 reports the in-sample parameter estimates of the preferred HAR-MJ-D specification. These estimates provide insight into volatility persistence across heterogeneous horizons, the role of jump dynamics, and the impact of residual weekly seasonality. The results are summarized in Table 5.

Restricting attention to statistically significant regressors, the fitted HAR-MJ-D model can be written as

$$\begin{array}{cc}R{V}_t=& -0.0298+0.3284 {\mathrm{MedRV}}_t^d+0.2601 {\mathrm{MedRV}}_t^w+0.4564 {\mathrm{MedRV}}_t^m \\ & +0.2324 J_t+0.2133 D_{6,t}+0.4370 D_{7,t}+{\epsilon }_t,\end{array}$$

(40)

where $J_t$ denotes the jump component and $D_{6,t}$ and $D_{7,t}$ correspond to Saturday and Sunday, respectively, with Monday as the reference category.

The in-sample estimation of the HAR-M-J, in Table 5, (Heterogeneous Autoregressive model with MedRV and Jumps), augmented with weekday dummies, provides significant insights into the volatility structure of the series. With an Adjusted R² of 0.3101, the model captures approximately 31% of the total variation in future realized volatility, indicating a robust fit for high-frequency financial data.

The model results strongly support the Heterogeneous Market Hypothesis, where agents with different time horizons contribute to the volatility process.

Daily Component: The daily median realized volatility is highly significant, confirming a strong short-term momentum effect.

Weekly Component: This remains significant at the 1% level, reflecting the influence of medium-term market participants.

Monthly Component: The monthly component exhibits the highest coefficient. This suggests a “long-memory” property in the series, where long-term structural trends—could be linked to the geopolitical and renewable factors— influences on the current volatility regime.

The Jump component is statistically significant at the 1% level. This confirms that the continuous part of volatility is insufficient to describe the market’s behavior. Despite the Ullrich demeaning and using a robust estimator. Discontinuous price shocks (jumps) possess their own predictive power and contribute significantly to the escalation of future volatility. This indicates that sudden “bursts” of information or unexpected events are quickly priced into the volatility of subsequent periods.

To account for potential day-of-the-week effects, the model includes dummy variables, using Monday as the baseline reference.

Working days: The coefficients for Tuesday through Friday are not statistically significant. This suggests that volatility during the standard workweek does not deviate meaningfully from the levels observed on Monday.

The Weekend Effect (Sat–Sun): In contrast, the coefficients for Saturday and Sunday are both positive and highly significant. The magnitude of the Sunday coefficient is particularly striking, being the largest dummy in the model. This reveals a substantial volatility during the weekend.

Overall, the results reported in this section provide consistent evidence that volatility dynamics in the Spanish electricity market are best captured by models that explicitly account for long-memory behavior, discontinuous price movements, and residual calendar effects. The out-of-sample evaluation based on the QLIKE loss function and the Model Confidence Set highlights the superiority of HAR-type specifications over GARCH models, while conditional diagnostics further indicate that incorporating jump components and day-of-the-week effects improves forecast accuracy without introducing bias or inefficiency. Taking together, these findings suggest that reliable volatility forecasting in electricity markets requires models capable of exploiting high-frequency information and accommodating structural features unique to power prices. The implications of these results for risk management, market monitoring.

Importantly, volatility forecasting in electricity markets is fundamentally constrained by the dominance of jump-driven and exogenous shocks, implying that very high $R_{\mathrm{OOS}}^2$ values would be unexpected and potentially indicative of overfitting. Consequently, moderate out-of-sample explanatory power is fully consistent with the realized volatility forecasting.

5. Conclusions

This paper analyzes volatility forecasting in the Spanish electricity spot market over the period 2021–2025, characterized by uncertainty, frequent price jumps, and the presence of zero and negative prices. To address these features, we combine a price-shift transformation and the Ullrich demeaning procedure [10] with realized volatility measures constructed from high-frequency data, and evaluate a wide range of competing volatility models within a coherent out-of-sample forecasting procedure.

Model comparison is conducted using a loss-based approach using QLIKE loss function and the Model Confidence Set methodology. This procedure provides statistically robust inference on predictive performance while avoiding multiple-testing and data-snooping biases. The empirical results show that all return-based GARCH specifications are excluded from the Superior Set of Models, whereas HAR-type models based on realized volatility measures systematically dominate in terms of expected loss.

Conditional diagnostics within the Superior Set further reveal that models incorporating jump components and day-of-the-week effects achieve superior economic performance. In particular, the HAR-MJ-D specification based on Median Realized Volatility combines the highest out-of-sample $R_{\mathrm{OOS}}^2$ with unbiased and efficient forecasts according to Mincer–Zarnowitz tests. These results highlight the importance of modeling long-memory volatility components, discontinuous price movements, and residual weekly seasonality in modern electricity markets.

Although Out-of-Sample $R_{\mathrm{OOS}}^2$ values are moderate in absolute terms, this magnitude is expected in electricity markets, where volatility dynamics are dominated by jumps, weather shocks, and structural changes. Given the strength of the historical mean benchmark and the high intrinsic unpredictability of power prices, the reported $R_{\mathrm{OOS}}^2$ values represent substantial and economically meaningful forecast gains.

Overall, the findings underscore the relevance of coherent loss-based evaluation procedures and high-frequency information for volatility forecasting. And the importance of demeaning, using robust estimators and modeling the day-effects. Future research could extend this approach to multiple time steps, using dummy variables to account for structural changes. The effect of climatic variables—such as dummy variables for periods of drought or wind speed, among others—could be examined.