A Study on the Nature of Complexity in the Spanish Electricity Market Using a Comprehensive Methodological Framework

Abstract

The existence of chaos is particularly relevant, as the identification of a chaotic behavior in a time series could lead to reliable short-term forecasting. This paper evaluates the existence of nonlinearity and chaos in the underlying process of the spot prices of the Spanish electricity market. To this end, we used daily data spanning from 1 January 2013, to 31 March 2021 and we applied a comprehensive framework that encompassed a wide range of techniques. Nonlinearity was analyzed using the BDS method, while the existence of a chaotic structure was studied through Lyapunov exponents, recurrence plots, and quantitative recurrence analysis. While nonlinearity was detected in the underlying process, conclusive evidence supporting chaos was not found. In addition, the generalized autoregressive conditional heteroscedastic (GARCH) model accounts for part of the nonlinear structure that is unveiled in the electricity market. These findings hold substantial value for electricity market forecasters, traders, producers, and market regulators.

1. Introduction

In all European markets, the day-ahead market session is held at 12:00 CET, where prices and electrical energies for the twenty-four hours of the next day are set. The price and volume of energy at a specific hour are established in such a way that supply and demand match, following the model agreed upon and approved by all of the European markets. Buy (sell) orders are accepted in order of increasing (decreasing) prices until the total demand (supply) is met [1]. Additionally, in each system there are several sessions in the intraday market to adjust, in real time, the bids made the previous day. Despite the existence of these sessions in the intraday market, the reliable prediction of electricity prices remains a major challenge to achieve the goal of maximizing profits for both electricity producers and buyers. In the case of electricity, the intrinsic difficulties of any prediction on the evolution of its price are increased by certain specific characteristics that differentiate it from other commodities. In this regard, it is worth noting the high volatility of electricity prices, which is largely due to the impossibility of storing it and to fluctuations in its demand [2].

Given its interest and complexity, it is unsurprising that the literature review reveals many works employing numerous methods to forecast wholesale electricity market prices [3,4,5,6,7]. As a demonstration of this, in a comprehensive review, Weron [1] noted that, up to 2014, there were a total of 304 publications related to this topic in the Web of Science (WoS), with 136 being journal articles. Hence, a literature review on electricity price forecasting is outside of the scope of this paper. Other research on electricity prices has focused on other aspects such as market performance multifractality and price dynamics over time [8,9,10,11,12,13,14,15]. Interestingly, Strozzi et al. [3] applied recurrence quantification analysis on NordPool spot prices and compared the volatility found in the series with the metrics calculated from the recurrence plots. In this case the analysis focused on hourly prices. Papaioannou et al. [15] compared the efficiency of several European electricity markets from an index composed of linear and nonlinear measures that quantifies aspects like complexity, entropy, or long-run memory of the series.

In the context of forecasting, the identification of chaos is particularly significant as it allows for reliable short-term predictions, although not in the medium or long term [16,17]. Moreover, if the series possesses a chaotic structure the traditional linear forecasting methods do not work properly, and nonlinear models should be applied [18]. For this reason, numerous studies have been conducted to analyze chaotic behavior in macroeconomic and financial series, yielding mixed results regarding its presence [19]. However, this issue has received reduced interest in the electricity market literature. It is also important to highlight that drawing an overall conclusion from studies on the presence of chaos in electricity prices is challenging, given the diversity of commercial regulations worldwide [20]. For example, one of the works in this field is that of Bigdeli and Afshar [21] who studied the oscillations in the price series of the Iranian electricity market based on recurrence plots. The authors detected evidence of deterministic chaos in the system with significant short-term predictability. However, it is important to note their conclusions should not be taken as a reference for this research due to its characteristics: relying on a relatively small dataset, using a single method to detect chaos and, moreover, mainly due to the fact that it was carried out for the Iranian intraday electricity market, which has very different characteristics from those of the European markets.

The main purpose is to investigate the underlying dynamics of the Spanish electricity market price series (SPEL) and to study whether they come from a stochastic system or, on the contrary, from a deterministic chaotic system. By using daily data covering the period from 1 January 2013 to 31 March 2021, we analyzed the presence of nonlinearity and chaos. Nonlinearity analysis technique is used, as well as metric and topological methods of chaos detection that are suitable for economic and financial series. The considered methodology considers various features that characterize chaotic time series, such as nonlinearity, sensitivity on initial conditions and recurrence. Each attribute is examined through its respective technique. The presence of nonlinearities is assessed using the BDS method suggested in [22,23]. Similarly, the Lyapunov Exponent test, modified for dynamic systems with noise [24], along with the application of a recurrence plot (RP) in conjunction with Recurrence Quantification Analysis (RQA), have been employed to assess the existence of chaos.

The overall methodological framework of this work is based on the one proposed by Inglada-Pérez [25,26] which has been specifically extended for this research by incorporating the following processes, among others: (i) the possible existence of seasonality in the time series that could affect the results and conclusions of the study is analyzed; (ii) the Lagrange multiplier test is incorporated to contrast whether the presence of nonlinearity is due to the existence of heteroscedasticity; and (ii) in the stage of the detection of a chaotic regime, the methodological tool of Recurrence Quantification Analysis is introduced, which, through the determination of several ratios, apart from providing information on other characteristics of the series, due to its quantitative nature, significantly improves the power of the recurrence plots to detect the existence of chaos. In summary, the methodological process is as follows:

Initially, as some of the later-used techniques necessitate stationary time series, unit root tests are conducted to assess the stationarity of the SPEL series. Once it is confirmed that the SPEL series is I(1), i.e., the return series is stationary, the analysis proceeds on the daily return series computed as the first difference of the logarithmic price level. Recognizing the presence of weekly seasonality, a seasonal difference is applied to the return series. Linear dependence is then filtered out of the data by applying seasonal autoregressive moving average (SARMA) filters based on the Box–Jenkins methodology [27]. Subsequently, to identify nonlinearity we employed the BDS method to the residuals previously obtained. If nonlinear dependence is noticed, it suggests the presence of some form of nonlinear dependence in the return series. Given that linear structures have already been removed using the best-fit SARMA model, the nonlinearity may arise from either a nonlinear stochastic system or a nonlinear deterministic chaotic system. Stochastic nonlinearity may be attributed to the presence of a volatility cluster, i.e., large variations tend to cluster together, a common characteristic of financial time series. The Lagrange multiplier test is used to check if nonlinearity is caused by the presence of heteroscedasticity in the data. In such a scenario, a generalized autoregressive conditional heteroskedasticity model (GARCH) [28] is fitted.

Following this, the presence of chaos is inspected through a set of robust procedures, including two particularly well-suited for the analysis of time series with noise: (i) the Lyapunov exponent test using artificial neural networks algorithm to estimate the largest Lyapunov exponent, and (ii) a recurrence plot (RP) together with the tool called Recurrence Quantification Analysis (RQA) in order to clarify the interpretation of the information contained in the RP and to increase knowledge.

Hence, the novelty of this research compared to previous investigations and its main contribution to enrich the scarce existing literature on the presence of chaos in the electricity market is that, to the best of our knowledge, it is the first time that the existence of chaotic structure in a crucial variable in the economy—the electricity price—is explored from a multidimensional perspective that considers all of the properties that characterize a chaotic system and using a series of robust metric and topological procedures, including those that are suitable for noisy time series analysis, i.e., the Lyapunov exponent and recurrence plots in combination with Recurrence Quantification Analysis. We also explore other topics such as the existence and modeling of volatility clustering and the presence of seasonality. Findings indicate the presence of nonlinearity in the underlying process of the SPEL, although no conclusive evidence supporting chaos was found. Our results have substantial value by providing insights into the inherent behavior and dynamics of electricity prices. Overall, our insights should be relevant to all stakeholders involved in the power market, including, productors, brokers, regulators, forecasters, and system operators, among others.

The remainder of the paper is structured as follows: The second section outlines the methodological process. In the next section, we present the data used and the main results obtained. Additionally, we analyze and discuss the existing evidence both supporting and refuting the hypotheses of nonlinearity and chaos. Finally, the concluding section reports a series of conclusions and suggests potential avenues for future investigation.

2. Methods

2.1. Methodological Framework

The methodology used is built upon the one introduced by Inglada-Perez [25,26], encompassing these steps:

• (i). Stationarity testing:

  This was examined using the following methods:

  • ▪. Augmented Dickey–Fuller Test (ADF) [29];

  • ▪. Phillips–Perron Test (PP) [30];

  • ▪. Kwiatkowski–Phillips–Schmidt–Shin Test (KPSS) [31].

Once we determined that the series under examination was not stationary, the log-returns series (from now on referred to as returns series) was obtained as described: $ln\left(R_t\right)=ln\left(P_t\right)-ln\left(P_{t-1}\right)$ such that $P_t$ represents the time series data at time t. The return series possesses crucial properties such as stationarity, an essential property since it is a prerequisite for most of procedures used further on.

• (ii). Linear modelling:

After confirming that the return series is stationary, the analysis is conducted on this series. Upon detecting a weekly seasonal pattern, seasonal differencing is applied to the return series. Subsequently, the BDS test is employed to examine the existence of nonlinearity. To achieve this, it is essential to eliminate the linear structure of the series. Linear dependence is removed by applying SARMA filters based on the Box–Jenkins methodology [27]. The SARMA model is an extension of the ARMA (Autoregressive Moving Average) model, particularly designed for time series that exhibit a seasonal pattern, as is the case of our time series. We used a multiplicative SARMA model, or SARMA(p, q)x(P,Q)S, of the type

(1)${\varphi }_p\left(B\right){\Phi }_P\left(B^S\right)Y_t=\mu +{\theta }_q\left(B\right){\mathsf{\theta }}_Q\left(B^S\right) a_t,$

where $Y_t$ is the series obtained by applying a seasonal difference to the returns series, so that it is stationary. B is the backshift operator, S indicates the order of seasonality, ${\varphi }_p\left(B\right)=1-{\varphi }_{1}B\dots -{\varphi }_p{B}^p$ and ${\theta }_q\left(B\right)=1+{\theta }_{1}B\dots +{\theta }_q{B}^q$ denote the non-seasonal polynomials of order p and q in B respectively, ${\Phi }_P\left(B^S\right)=1-{\Phi }_{1}B\dots -{\Phi }_p{B}^{PS}$ and ${\mathsf{\theta }}_Q\left(B^S\right)=1+{\mathsf{\theta }}_{1}B\dots +{\mathsf{\theta }}_q{B}^{QS}$ denote the seasonal polynomials of order P and Q in $B^S$ respectively, $a_t$ is the error term or random disturbance and $\mu$ is a constant.

The selected SARMA model was the one that obtained the minimum value according to the Schwarz Information Criterion (SIC) [32]. After establishing the optimal linear model, the residuals are obtained. These newly acquired residuals are referred as the SARMA series.

• (iii). Nonlinearity:

Next, we applied the BDS method [22,23] to the SARMA series to examine the presence of nonlinearity. In the case that the BDS test rejects the null hypothesis of linearity, the Lagrange multiplier test (LM) for autoregressive conditional heteroskedasticity (ARCH) proposed in [33] is applied to the residual series. This is conducted to investigate whether the detected nonlinearity can be attributed to the existence of heteroscedasticity. The test fits a linear regression to the squared residuals and assesses the significance of the fitted model. The null hypothesis is whether the squared residuals are white noise, indicating homoscedasticity in the SARMA series. If the null hypothesis is rejected, it suggests evidence that the nonlinearity may be driven by the presence of heteroscedasticity in the series.

• (iv). Volatility:

If nonlinear dependence is detected in the residual series from SARMA models, it may be attributed to the presence of a volatility cluster. To verify this, the conditional variance is modeled by fitting ARCH-type models (GARCH). The specific model is determined as the SARMA model was. Afterwards, the standardized residuals are obtained from their conditional standard deviations (hereafter GARCH series). The BDS test is then applied to this new GARCH series. This step serves to examine whether the process has successfully removed all nonlinear dependence. If no dependence is detected, it implies that the data is not chaotic, as the ARCH-type model has effectively captured all nonlinearities.

• (v). Study of chaos:

  The examination of a chaotic structure was conducted using these procedures:

  • ▪. Lyapunov’s Exponent, which analyzes the system’s sensitivity to initial conditions [24].

  • ▪. The recurrence plot, a visual tool that permits an analysis of the existence of periodic patterns, among other aspects of time series [34]. In addition, Recurrence Quantification Analysis is utilized to provide clarity on the information presented by the recurrence plots.

These techniques are described below.

2.2. Reconstructing the State Space

Before delving into the main points of the employed techniques, it is worth depicting how the state space can be reconstructed. The embedding theorem proposed by [35] provides a means to recover the system dynamics, i.e., the original trajectory of an unknown dynamical system that gave rise to an observed scalar time series. This is achieved by reconstructing a new state space using successive values of the time series. Utilizing Takens’ method, a reconstructed phase space, which is topologically equivalent to the original one, is obtained from an m-dimensional dynamic system, of which only the evolution of a state variable is known, given by the series ${\left\{x_t\right\}}_{t=1}^T$. The reconstructed space is obtained using delayed versions of the time series a specific number of times. In other words, it is formed by grouping the values of the time series, thus obtaining vectors where each coordinate corresponds to the values of the original series and their respective delays. Specifically, we construct a sequence of delayed vectors by associating for each time period a vector in a reconstructed state space Rm, whose coordinates satisfy this equation:

(2)$X_t^m=\left(x_t,\dots .,x_{t+m-1}\right),$

where $\tau$ is the lag or reconstruction delay, and m is the embedding dimension of the reconstructed space.

If the embedding dimension $m$ is sufficiently larger than the attractor dimension, Takens’ theorem guarantees the existence of a diffeomorphism between the original and the reconstructed attractor. In essence, the theorem implies that both attractors represent the same dynamical system in different coordinate systems if $m\ge 2n+1$. In the embedding space, numerous invariants of the dynamics remain conserved.

2.3. Linear Analysis: BDS Test

The main hypotheses and features that characterize this method are described below. Herein the BDS method is used to analyze nonlinearity. While the BDS test is originally designed to examine independence rather than chaos or nonlinearity, various studies, including [36,37], have illustrated its capability to detect nonlinearity. Specifically, when it is applied on the residuals from a linear model (ARMA process), the rejection of the null hypothesis, which assumes that the series is independent and identically distributed (IID), indicates the existence of linear dependence. So, it can be used to identify the dependencies that may persist in the data, and also to discern the existence of nonlinear structures. The BDS test is usually applied to residuals to check whether the best-fit model for a time series is linear or nonlinear. It is constructed using the algorithm outlined in [38].

As mentioned above, in accordance with Takens’ theorem the method of delays can be used to embed a scalar time series ${\left\{x_t\right\}}_{t=1}^T$. Therefore, the reconstructed state vector or m-history is $X_t^m=\left(x_t,x_{t+\tau },\dots .,x_{t+\left(m-1\right)\tau }\right)$ where $\tau$ stands for the lag or reconstruction delay, and m is the embedding dimension of the reconstructed space. The reconstructed trajectory, X, of N m-dimensional vectors $X_t^m$ can be expressed as a N × m matrix with $N=T-\left(m-1\right)\tau$.

Under the null hypothesis, the BDS statistic is computed from the correlation integral [39], which measures the frequency at which temporal patterns are reiterated in the data. The correlation integral at the embedding dimension m is defined by

(3)$C\left(T,m,\epsilon \right)=\frac{2}{N\left(N-1\right)}{\displaystyle \sum }_{t<s}H\left\{ϵ-\Vert X_t^m-X_s^m\Vert \right\},$

where T, ‖.‖, m and H represent the size of the data sets, a norm operator, the embedding dimension and the Heaviside function. Furthermore, $N=T-\left(m-1\right)\tau$ is the number of embedded points in the m-dimensional space.

$C\left(T,m,\epsilon \right)$ measures the fraction of pairs of points $X_t^m$, t = 1; 2; …, T, whose separation is smaller than $ϵ$. If the limit of $C\left(T,m,\epsilon \right)$ as $T\to \infty$ exists for each $\epsilon$, we denote the fraction of all state vector points that are within $\epsilon$ of each other as $C\left(m,\epsilon \right)=\underset{T\to \infty }{\lim }C\left(T,m,\epsilon \right)$.

If the data originate from a strictly stationary stochastic process that is absolutely regular, then this limit exists [28]. Moreover, in situations where the process is IID and since $\mathsf{\Theta }\left(\epsilon -\Vert X_t^m-Y_t^m\Vert \right)={\displaystyle \prod }_{k=1}^m\mathsf{\Theta }\left(\epsilon -\left|x_k-y_k\right|\right)$, it follows that $C\left(m,\epsilon \right)=C^m\left(1,\epsilon \right)$.

Likewise, the difference $C\left(m,\epsilon \right)-C^m\left(1,\epsilon \right)$ follows an asymptotic normal distribution, with a zero mean and variance provided by:

(4)$\frac{{\sigma }^2\left(T,m,\epsilon \right)}{4}=m\left(m-1\right) C^{2\left(m-1\right)} \left(K-C^2\right)+K^m-C^{2m}+2{\displaystyle \sum }_{i=1}^{m-1}\left[C^{2i}\left(K^{m-1}-C^{2\left(m-1\right)}\right)-m{C}^{2\left(m-i\right)}\left(K-C^2\right)\right]$

We can consistently estimate the two constants in the expression above: C by $C\left(1,\epsilon \right)$ and K by $K\left(T,m,\epsilon \right)=\frac{6}{T\left(T-1\right)\left(T-2\right)}{{\displaystyle \sum }}_{t<s}\left[\mathsf{\Theta }\left(\epsilon -\Vert X_t^m-X_s^m\Vert \right)\mathsf{\Theta }\left(\epsilon -\Vert X_t^m-X_s^m\Vert \right)\right]$.

Under the IID null hypothesis, the BDS statistic is defined as

(5)$BDS\left(T,m,\epsilon \right)=\frac{\sqrt{T}}{\sigma }\left[C\left(m,\epsilon \right)-C^m\left(1,\epsilon \right)\right]$

This statistic follows a standard normal distribution under the null hypothesis of IID as $T\to \infty$ and its critical values are obtained using the standard normal distribution.

We use this test to investigate the nonlinearity of the price returns series, after removing the linear dependencies in the series. As in [37], the test is performed for embedding dimensions between 2 and 8 as well as for different values of $\epsilon$. Specifically, we consider ε values of 0.5σ, σ, 1.5σ, and 2σ, where σ represents the standard deviation of the residual series [36]. It is important to note that the series comprises approximately 3000 observations. This is relevant since the distribution of the BDS statistic is only demonstrated for samples of 500 or more observations [36]. Thus, to obtain a reliable analysis, large data sets are needed.

2.4. Lyapunov Exponents

The main hypotheses and features that characterize this method are described below. A distinctive property of chaotic systems is their sensitivity to initial conditions. Lyapunov exponents (LE) are a widely adopted method for quantifying this sensitivity, providing a measure of the rate of divergence of closely spaced trajectories whose initial conditions only differ by an infinitesimal amount. Within the international literature on chaos, it is observed that LE stands out as one of the most extensively employed criteria for investigating the presence of chaos in time series [40,41]. In the context of a bounded system, the presence of a positive LE is an indicator of chaos in nonlinear dynamical systems exhibiting seemingly random and unpredictable behavior. Therefore, the computation of the maximal LE is as a valuable tool for detecting a chaotic behavior [42].

As mentioned above, according to Takens’ theorem the method of delays can be used to embed a scalar time series ${\left\{x_t\right\}}_{t=1}^T$. The reconstructed state vector or m-history is $X_t^m=\left(x_t,x_{t+\tau },\dots .,x_{t+\left(m-1\right)\tau }\right)$ where $\tau$ is the lag or reconstruction delay, and m is the embedding dimension of the reconstructed space. The reconstructed trajectory, X, of N m-dimensional vectors $X_t^m$ can be expressed as a N × m matrix with $N=T-\left(m-1\right)\tau$.

The state space reconstruction procedure outlined above enables us to extract all of the invariant properties from the time series related to the unknown dynamical system that generates the data. In fact, the LE are expected to exhibit approximately the same values in both the true and the reconstructed state space. This observation empowers us to test the existence of chaos in the original unknown dynamical system by assessing the sign of the maximum LE obtained from the time series.

Literature reveals two types of methods for estimating λ from experimental data, based on state space reconstruction using delay coordinate embedding approximations. Direct methods measure the divergence growth rate between two neighboring trajectories with an infinitesimal difference in their initial conditions [43]. A notable limitation of this approach is its sensitivity to both sample size and the existence of noise in the data [44]. Furthermore, these methods do not allow for statistical inference about chaos, as the asymptotic distribution of the estimator does not exist.

The indirect approach, also known as the Jacobian-based method, overcomes the limitations of direct methods. This method can yield robust estimates of the LE despite existing noise [45,46]. Also, the asymptotic distribution of the estimator can be determined [44], enabling the formulation of formal tests and inferences about the LE. In other words, the estimation of this Jacobian will let us to test the hypothesis of chaos characterized by a positive LE value.

With the latter approach, a model is initially constructed for the data by relying on approximations of the trajectories in the reconstructed state space. Next, the Jacobian matrices of the estimated dynamical system are used to calculate the LE. However, since the map that generates the process is typically unknown, it is not possible to obtain the Jacobian matrix. Thus, the LE cannot be calculated. Hence, it is required to approximate the unknown chaotic map using some form of function. Searching for the best option, McCaffrey et al. [45] compared the results of several methodologies including neural networks, radial basis functions, projection search, and thin-plate splines to approximate this unknown map. In light of the findings obtained, it was concluded that the neural network performed as the most suitable regression method for chaotic systems in which noise is present. In fact, Bailey et al. [46] investigated a regression approach employing neural networks. The outcomes they achieved for a Henon noisy system with noise showed that the neural network regression method yielded precise estimates for the LE.

Interestingly, Hornik et al. [47] demonstrated that any standard feed-forward network with only one hidden layer can approximate any measurable Borel function from one finite-dimensional space to another for any desired level of precision, given an adequate number of hidden units. Based on this, an algorithm relying on a neural network with a single hidden layer has been considered in this research.

Following [24,48], this technique encompasses these stages: let ${\left\{x_t\right\}}_{t=1}^T$ be a time series. A chaotic system in which noise is present can be formulated in the following way:

(6)$X_t=f\left(x_{t-\tau },x_{t-2\tau },\dots , x_{t-m\tau }\right)+{\epsilon }_t,$

where t denotes the time script and m, τ, ${\epsilon }_t$, and f stand for the embedding dimension, the time delay, noise added to the series, and an unknown chaotic map, respectively.

The LE is then defined as follows:

(7)$\widehat{\lambda }=\frac{1}{2M}\ln \left(v_1\right),$

where $M\le T$ is the “block-length” or number of equally spaced evaluation points employed for computing the LE and $v_1$ is the largest eigenvalue of the matrix ${\left(T_M{U}_0\right)}^{\prime }\left(T_M{U}_0\right),$ being $T_M={\displaystyle \prod }_{t=1}^{M-1}{J}_{M-t},$ $J_t$ the Jacobian matrix and $U_0={\left(1, 0, \dots , 0\right)}^{\prime }$. As we generally lack information about the chaotic map F, it is not possible to estimate the Jacobian matrix. Hence, it is necessary to approximate the unknown chaotic map by an appropriate function. The different approaches that constitute the indirect method vary in the algorithm applied for the estimation of this function.

In our algorithm the chaotic process is estimated by the following equation:

(8)$x_t\approx {\alpha }_0+{\displaystyle {{\displaystyle \sum }}_{j=1}^q}{\alpha }_j\Psi \left({\beta }_{o,j}+{\displaystyle {{\displaystyle \sum }}_{i=1}^m}{\beta }_{i,j}{x}_{t-i\tau }\right)+{\epsilon }_t,$

where ${\alpha }_0$ is the estimated network bias from the input; q is the number of neurons (or nodes) in the single hidden layer; ${\alpha }_j$ are the estimated inter-layers connection weights from the input to the hidden layer; $\mathsf{\Psi }$ is the transfer or activation function, which in our case is the logistic function $\mathsf{\Psi }\left(u\right)=\frac{1}{1+e^{-u}}$ as usual [24], ${\beta }_{o,j}$ is the estimated network bias from the hidden layer, m is the embedding dimension, and ${\beta }_{i,j}$ are the estimated layers connection weights from the hidden layer to the output.

The matter of parameter estimation is simplified into a least squares problem. The quantity to be minimized is the following one:

(9)${\displaystyle \sum }_{t=1}^T{\left(x_t-\left[{\alpha }_0+{\displaystyle \sum }_{j=1}^q{\alpha }_j\Psi \left({\beta }_{o,j}+{\displaystyle \sum }_{i=1}^m{\beta }_{i,j}{x}_{t-i\tau }\right)\right]\right)}^2$

Shintani and Linton [44] analyzed the asymptotic properties of the nonparametric neural network estimator of the LE. They constructed a statistical framework to analyze the chaos hypothesis based on that estimator and the consistent estimator of its variance. Their findings indicated the asymptotic normality of the LE estimator:

(10)$\sqrt{M}\left({\widehat{\lambda }}_k-{\lambda }_k\right)\sim N\left(0,{\phi }_k\right),$

where ${\phi }_k$ represents the variance of the kth LE estimator.

Additionally, they proved that a consistent estimator of the variance ${\phi }_k$ is as follows:

(11)${\widehat{\phi }}_k\equiv Var\left({\widehat{\lambda }}_k\right)=\underset{M\to \infty }{\lim }Var\left(\frac{1}{\sqrt{M}}\right){\displaystyle \sum }_{t=1}^M{\eta }_{k,t}$

with M being the subsampling size, the quadratic spectral kernel function ${\eta }_{k,t}$ which is calculated as ${\eta }_{k,t}={\xi }_{k,t}-{\widehat{\lambda }}_k$ and the parameter ${\xi }_{k,t}=\frac{1}{2}log{\mu }_k\left(\left|\widehat{D}{g}^t\right|\right)-\frac{1}{2}log{\mu }_k\left(\left|\widehat{D}{g}^{t-1}\right|\right)$.

The test statistics to test the hypothesis $H_0: {\widehat{\lambda }}_k>0$ versus $H_1:{\widehat{\lambda }}_k\le 0$ can be defined in the following manner:

(12)${\widehat{t}}_k=\frac{{\widehat{\lambda }}_k}{\sqrt{{\widehat{\phi }}_k/M}}\sim N\left(0,{\widehat{\phi }}_k\right)$

We will reject $H_0$ if ${\widehat{t}}_k\le -z_{\alpha }$, where $z_{\alpha }$ is the critical value that satisfies $\mathrm{Pr}\left[Z\ge z_{\alpha }\right]=\mathsf{\alpha }$ with $Z$ being a standard normal random variable and α being the significance level.

In summary, under the null hypothesis H₀ that chaos is present, the neural net estimator ${\widehat{\lambda }}_k$ leads to asymptotically valid inferences meaning that the associated p-value follows a normal distribution on $N\left(0,{\widehat{\phi }}_k\right)$. The rejection of the null hypothesis $H_0: {\widehat{\lambda }}_k\ge 0$ provides strong evidence in favor of no chaos.

2.5. Recurrence Analysis

The main hypotheses and features that characterize these methods are described below.

2.5.1. Recurrence Plot

The recurrence plot (RP), first proposed by [24,49], is a visual analysis tool for detecting recurrent patterns in time series. Although this tool was initially designed to investigate non-stationarity in time series, it has subsequently been used for numerous applications due to being a very suitable instrument for the analysis of short, non-stationary and noisy time series [9]. This procedure unveils the hidden structures of a time series through a qualitative approach and, specifically in this research, it is used to study both the nonlinearity of the data and the chaotic dynamics of the series.

By using Takens’ theorem, it is possible to obtain the RP of the series as well as some measures that allow us to quantify the patterns observed in the plot. The plot is constructed supposing that mutual distances belong to the same trajectory in the reconstructed phase space. It is obtained from the recurrence matrix [50] whose construction is performed using the following algorithm [25,51]:

(13)$R_{i,j}\left(ϵ\right)=\mathsf{\theta }\left(ϵ-\Vert x_i^m-x_j^m\Vert \right),$

where m is the embedding dimension, $X_i^m,X_j^m$, are vectors in the reconstructed state space, $i,j=1,\dots ,N-\left(m-1\right)\tau$; $\theta (-)$ is the Heaviside function, $ϵ$ is the distance threshold value and ||–|| is a norm. Hence, the recurrence matrix is generated by comparing each embedded vector $X_i^m$ with the other $X_j^m$. A point is drawn if this comparison is less than a value ε for a specific distance. That is, if the condition is met: $\Vert X_i^m-X_j^m\Vert <\epsilon$. The recurrence of a state between $t=i$ and $t=j$ occurs when $X_j^m$ is close enough to $X_i^m$. In summary, a point i, j of the matrix takes the value 1 when the states at time i and time j are close enough and takes the value 0 when they are not. Specifically, they are said to be close when the distance between the two states is less than a certain threshold. The RP is the result of plotting this matrix formed by zeros and ones where R_i,j = 1 is shown as a black point on the plot, while R_i,j = 0 is shown as a white one. The matrix M is symmetric due to the symmetry of the distance operator, so that the RP is likewise symmetric around the bisecting line of the first quadrant.

The diagonal structures of the RP indicate the degree to which one part of the path is relatively close to another at a specific time. The periodic motion manifests as long and non-interrupted diagonals. In a deterministic chaotic system, short lines parallel to the main diagonal are detectable whereas, in stochastic systems, short line segments are absent and uniformly distributed points are noticed. This method offers also supplementary information about the attractor’s structure, as the plot maintains the temporal order of the series, revealing the place and periodicity of the periodic orbits.

Therefore, the RP is a representation that shows when a system visits a state in which it has been in the past. The graph is symmetric, its axes show time, and the diagonal is a solid line. This main diagonal, known as the line of identity or LOI, shows the recurrence of a state with itself. Its appearance makes it possible to analyze the dynamics of the system, to study the presence and influence of noise, to detect states that recur over time and to observe abrupt changes in the dynamics of the system. Moreover, they remain valid for non-stationary or noisy series and do not depend on the sample size [19]. In particular, long continuous lines parallel to the LOI are usually related to systems with periodic or quasi-periodic behavior, while short, dashed diagonal lines usually indicate that the system is chaotic [50].

Therefore, in order to make a projection of the series in the reconstructed state space, it is necessary to first determine the dimension and the embedding delay. Fraser and Swinney [52] argue that the first minimum of the Mutual Information function is a good alternative to determine the embedding delay. On the other hand, the embedding dimension is usually determined using the false nearest neighbor algorithm [53].

2.5.2. Recurrence Quantification Analysis

The recurrence plot allows us to visualize the recurrence patterns existing in the time series; however, to complement this analysis it is interesting to obtain different measures that can quantify the patterns observed in the plot. Several complexity measures have been proposed to quantify these patterns and they form what is called Recurrence Quantification Analysis (RQA). This analysis was developed by [54] and subsequently by [53].

Some studies show that RQA measures are able to detect changes in the dynamics of the system. On the one hand, they allow the detection of transitions from chaotic periods to ordered periods, i.e., to periods with greater determinism [55]. They also reveal transitions between chaotic periods [56]. Transitions between chaotic periods are detected from measurements that quantify the vertical structures of the plot.

The measurements presented in this work are based on a calculation of the recurrent point density and frequency distribution of diagonal and vertical line lengths. In particular, the histograms of the lengths of these lines are fundamental pillars of recurrence quantification analysis. The measures associated with the diagonal lines of the recurrence plot are based on the following histogram of diagonal lines of length l.

(14)$P\left(l\right)={{\displaystyle \sum }}_{i=1,j=1}^N\left(1-R_{i-1,j-1}\right)\left(1-{\mathrm{R}}_{i+l,j+l}\right) {{\displaystyle \prod }}_{k=0}^{l-1}{R}_{i+k,j+k},$

Similarly, the histogram of the lengths of the vertical lines is used for measurements associated with vertical lines.

(15)$P\left(v\right)={\displaystyle \sum }_{i=1,j=1}^N\left(1-R_{i,j}\right)\left(1-R_{i,j+v}\right) {\displaystyle \prod }_{k=0}^{v-1}{R}_{i,j+k}$

From these histograms we can calculate most of the measures presented below, as presented in [51].

Recurrence rate (RR): This is a measure based on the density of recurrences and can be considered as the recurring probability of any state. Specifically, it measures the proportion of points considered as recurrences with respect to the total number of possible recurrences. The rate is calculated as follows:

(16)$RR=\frac{{{\displaystyle \sum }}_{l,j=1}^N{R}_{i,j}}{N^2}$

Determinism (DET): Determinism allows the quantification of how predictable the dynamics of a system is. Specifically, determinism is observed in the existence of diagonal lines parallel to the main diagonal in the recurrence plot [54]. The diagonal lines show that the system revisits the same sequence of states over time. Thus, a deterministic process will consist of long diagonal lines and a few isolated points, while a white noise process will show almost no diagonal lines and consist of many isolated points. As shown in the following formula, determinism is quantified as the proportion of points in the plot that form diagonal lines. To calculate it, we must first determine the minimum number of points necessary to consider a line as a diagonal line. Herein we take two points.

(17)$DET=\frac{{{\displaystyle \sum }}_{l=l_{min}}^{N}lP\left(l\right)}{{{\displaystyle \sum }}_{l=1}^{N}lP\left(l\right)} ,$

where l_min is the minimum length for a line to be considered a diagonal line. In our case l_min = 2. Note that for l_min = 1 the determinism is 1.

Maximum (L_max) and mean (L) length of diagonal lines: The maximum length (L_max) and mean length (L) of the diagonal lines of the plot are other measures that are also based on the diagonal lines and allow us to obtain information about how predictable the system is. As discussed above, a plot with long diagonal lines indicates that the system is moving through a sequence of states it has been in in the past, so the longer the diagonal lines, the more predictable the system is. These metrics are given by the following:

(18)$L=\frac{{{\displaystyle \sum }}_{l=l_{min}}^{N}lP\left(l\right)}{N_l} ; L_{max}=max\left({(l_i)}_{i=1}^{N_l}\right),$

where $N={\displaystyle \sum }_{l=l_{min}}^{N}lP\left(l\right)$ represents the total number of diagonal lines.

Laminarity (LAM): Similar measures can also be applied to vertical lines or bands. Vertical lines are observed when laminar states exist in intermittent regimes, i.e., when the system is trapped in one state for a certain period or when there are abrupt changes that affect the dynamics of the system. It is a measure of the graduality of the data variations over time. As for the diagonal lines, the measurements are based on the frequency distribution of the lengths of the vertical lines, used to obtain the proportion of recurrent points forming vertical structures (LAM). This proportion is expressed as follows:

(19)$LAM=\frac{{{\displaystyle \sum }}_{v=v_{min}}^{N}vP\left(v\right)}{{{\displaystyle \sum }}_{v=1}^{N}vP\left(v\right)}$

Trapping Time (TT): Just as the average length of diagonal lines is defined, the average length of vertical structures is defined as follows:

(20)$TT=\frac{{{\displaystyle \sum }}_{v=v_{min}}^{N}vP\left(v\right)}{{{\displaystyle \sum }}_{v=v_{min}}^{N}P\left(v\right)}$

Trapping Time estimates the average time the system remains trapped in a certain state. Note that to calculate these last two measures it is necessary to choose a value for vmin, i.e., the minimum length for a structure to be considered a vertical line. In this work we take vmin = 2, i.e., two points are enough to form a vertical structure.

Entropy (ENTR): Finally, the entropy measures the distribution of those line segments that are parallel to the main diagonal and reflects the complexity of the deterministic structure in the system. In an RP, the Shanon entropy [57], which refers to the probability of finding a diagonal line of length l in the graph, is often used. Systems with high entropy show an RP with diagonal lines of different lengths, indicating that there are several deterministic patterns taking place in the system over time. On the contrary, a system with low entropy will show diagonal lines of equal length. Marwan et al. [50] present the Shannon entropy as follows:

(21)$ENTR=-{\displaystyle \sum }_{l=l_{min}}^{N}p\left(l\right)\ln p\left(l\right),$

where p(l) is the probability of finding a diagonal line of length l in the recurrence graph, i.e., $p\left(l\right)=\frac{P\left(l\right)}{N_l}$.

3. Results

Utilizing the aforementioned methodology, we present the major results.

3.1. Data Analysis

We used the daily spot prices of the Iberian electricity market (for the Spanish area), known as SPEL, and determined the average of the hourly prices for each day. The dataset covers the period from 1 January 2013 to 31 March 2021, with 3012 observations being used overall. The dataset was sourced from the site http://www.omel.es (accessed on 1 February 2002). The series is shown in Figure 1.

Figure 1 illustrates that the series of wholesale electricity prices is highly volatile and is characterized by large fluctuations and spikes. It is worth noting that the main variations in wholesale electricity prices are due to climatic factors, changes in the prices of raw materials (mainly gas, oil, and coal), regulatory changes, and factors influencing economic activity. For example, a decrease in prices is evident at the end of the first quarter of 2021 and throughout the subsequent quarter due to the decrease in economic activity caused by COVID-19. It also displays an increase in prices driven by a rise in oil prices in the first quarter of 2017 and a decrease in prices caused by a regulatory change applicable as of 1 April 2013.

The main descriptive statistics of the series are displayed in Table 1. The SPEL series exhibits a negative skewness coefficient, suggesting a higher occurrence of larger values. Moreover, it has a leptokurtic distribution. According to the kurtosis, skewness, and Jarque–Bera test (p-value < 0.05), it can be inferred that data deviates from a normal distribution.

3.2. Returns

Unit root tests are performed to examine the stationarity of the SPEL series. We tested the null hypothesis that SPEL has unit roots, i.e., is non-stationary. In particular, we examine the stationary property by applying several procedures: the Augmented Dickey–Fuller test (ADF) [29], the Phillips–Perron test (PP) [30], and the Kwiatkowski–Phillips–Schmidt–Shin test (KPSS) [31]. Results obtained are shown in Table 2. It is observed that, for the original series of the Spanish electricity market, both ADF and KPSS tests do not give evidence to reject the null hypothesis of unit root existence at the 5% confidence level, however, the PP test shows evidence to reject the null hypothesis of being nonstationary. Therefore, we conclude that there is evidence of non-stationarity in the original series. Thus, we proceed to calculate the returns series by considering the logarithms of the ratio of two consecutives prices. Applying these same tests on the returns series, it is found that the data are already stationary under the three criteria (see Table 2).

From both Table 1 and Figure 2, which displays the time evolution of the returns, it can be concluded that the returns series exhibits a positive skewness coefficient. This suggests that the distribution is right skewed, indicating that larger upward jumps are more frequent than smaller downward jumps, and the tail to the right of the mean is longer than the one to the left. Moreover, it is worth noting the high value of the kurtosis coefficient, which implies that large jumps are more likely to happen more frequently in comparison with a normal distribution. Excess kurtosis suggests leptokurtic behavior in the returns series and is a common feature with financial series. In this line, the Jarque–Bera test rules out that the data come from a normal distribution. Furthermore, it is noteworthy that the kurtosis coefficient has a high value, indicating that large jumps are more likely to occur more frequently compared to a normal distribution. The excess kurtosis points toward leptokurtic behavior in the returns series and is a common characteristic in financial series. In line with this, the Jarque–Bera test rejects the hypothesis that the data follow a normal distribution.

3.3. Linear Dependence

The periodogram and the autocorrelation function of the returns series are obtained to investigate the existence of a seasonal component in the data, and to determine if it is necessary to apply a seasonal difference to eliminate it. As depicted in Figure A1, shown in Appendix A, the periodogram of this series shows a clear peak when the frequency is 0.143. Therefore, there is evidence of weekly seasonality. The autocorrelation function confirms the existence of weekly seasonality of order = 7, as already detected in the periodograms. Specifically, it is found that, on average, prices are reduced on the weekend, possibly attributed to a decline in the economic activity and, consequently, a decrease in the demand for electricity. On the other hand, there is no evidence of other types of seasonality.

As certain techniques, like the BDS test, are sensitive to the presence of linear structure, the linear dependence was removed by employing a SARMA filter after applying seasonal differencing to the returns series. The series with which we work is stationary, so the necessary hypotheses to apply the Box–Jenkins methodology [27] are fulfilled. The SARMA model that best fits the series is obtained. To choose the best SARMA model, the Schwarz information criterion (SIC) is used, selecting the model that minimizes this value. Following this methodology, the appropriate model to fit the series is a SARMA(5,1)(1,2)₇. Details of the model are given in Appendix A (Table A1). The residuals series (SARMA) is stationary (see Table 3).

3.4. Nonlinear Dependence

Nonlinearity is a necessary, although not sufficient, condition to detect the presence of chaos [51]. Nonlinearity is studied by means of the BDS test. This test was employed after removing linear dependence, as recommended by [37]. In this manner, it functioned as an indirect technique for scrutinizing nonlinearity: rejecting the null hypothesis implies the affirmation of nonlinearity in the series. There are two key parameters when running the test: the embedding dimension m and the threshold term ε. Following [37], the test is performed for embedding dimensions ranging from two to eight and for different values of ϵ. As recommended by [60], four different values are considered, namely 0.5σ, σ, 1.5σ, and 2σ, σ being the standard deviation of the series. Table A3 in Appendix A shows that for all dimensions m, and for all distances ϵ, the null hypothesis of independence is rejected, leading to the conclusion that the SARMA series is nonlinear.

3.5. Volatility Clustering

As the BDS test provides evidence that the SARMA series is nonlinear, we modelled the conditional variance by fitting an ARCH family model. The best model selected based on SIC criteria for the conditional variance is a GARCH(1,1). Next, standardized residuals were acquired, hereafter referred to as the GARCH series. Additional details regarding the model are presented in Appendix A (Table A2).

As nonlinearity is a necessary, but not sufficient, condition for chaotic behavior, its existence in the GARCH series was analyzed using the BDS test. Results are outlined in Appendix A (Table A4). In broad terms, the linear hypothesis is not rejected, but instances of nonlinearity were identified. The BDS test reveals some remaining dependence in the series, which is consistent with the potential presence of a chaotic component in the GARCH model residuals.

3.6. Chaotic Behavior Analysis

The following methods were used to test the presence of chaotic behavior.

3.6.1. Lyapunov Exponent

These were employed to assess the sensitivity to initial conditions, a characteristic of chaotic processes. The LE (λ) of a dynamical system is a quantity that describes the rate of the divergence of trajectories that are infinitesimally close. The largest Lyapunov exponent determines a notion of predictability for a dynamical system. Under certain conditions, a system is defined as chaotic when the maximum LE is positive, while a negative exponent suggests a mean reverting behavior. We calculated the maximum LE using an algorithm grounded in neural networks which approximates the nonlinear function of the underlying system dynamics. This technique can deal with the presence of noise in the data. The null hypothesis is H₀: λ ≥ 0, suggesting the presence of chaos, while the alternative hypothesis is H₁: λ < 0, indicating the absence of chaos.

The estimated largest LE λ is negative for returns, SARMA, and GARCH series (see Table 3). The results show that the null hypothesis of the presence of chaos is strongly rejected at a 5% significance level as indicated by the computed p-values. In conclusion, there is no evidence suggesting the presence of chaos in the Spanish electricity market.

3.6.2. Recurrence Analysis

Recurrence Plots

Creating a recurrence plot (RP) allows us to capture the recurrence property of states, a characteristic feature of chaotic time series. Thus, with this tool we analyzed the level of complexity and the presence of chaos. The analysis of the original series using an RP is justified, as this tool does not necessitate the data to be stationary. An RP is shown in Figure 3 for the original, returns, SARMA, and GARCH series.

For the graphical representation, the embedding dimensions and delays were determined, as well as the value of the threshold that determines whether a point is considered recurrent. The methodology described in Section 2.5.1 was used. Specifically, the optimal values of m, the embedding dimension, and τ, the considered delay time, were estimated using the mutual information criterion and the false neighbours algorithm, respectively, as previously described [52,53]. Finally, following [61], the threshold to determine when a point is considered recurrent was set by requiring the recurrence rate to below (0.1 to 5%). Once these values were calculated, the recurrence plots shown in Figure 3 were obtained.

In both of the plots of the returns and residuals of the SARMA model it is possible to observe clusters of recurrent points separated from each other by horizontal and vertical white bands. Moreover, no isolated diagonal lines parallel to the mean diagonal or line of identity (LOI) are identified. It should be noted that chaotic dynamics are detected in the recurrence plots when small diagonal lines parallel to the LOI are observed, a feature that is not present in any RP. Neither is there any indication of chaos, since the previously described characteristic is not observed. Hence, after the analysis by recurrence plot, it is possible to affirm that there is not a chaotic behavior in the series.

As for the plots of the standardized residuals of the GARCH models, although they are mainly formed by isolated points characteristic of a stochastic system, some patterns in the form of vertical bands can still be observed. Note that these bands are observed in the plots of the original series and are maintained in both the returns and residuals series. This could be due to the existence of some pattern intrinsic to the series, such as changes in regulation or market structure. For example, in all of the plots there is a vertical white band somewhat after t = 2500, corresponding to 14 November 2019. This could be related to the approval of Regulation (EU) 2019/944 concerning the internal electricity market. This regulatory framework, among other objectives, aims to promote the opening of the daily and intraday electricity markets. It is possible that this regulatory change, characteristic of the series itself, is shown throughout all of the plots.

Finally, the plots of the original series show quite different patterns from those of the rest of the series. This is because they still incorporate trends and seasonality. In any case, small diagonal lines parallel to the LOI, which would be characteristic of a chaotic system, are also not observed.

Recurrence Quantification Analysis

Although the RP is a visual tool enabling the exploration of the general structure and main characteristics of the series, its interpretation is not straightforward. Therefore, Recurrence Quantification Analysis (RQA)—a statistical quantification of the RP based on a set of measures extracted from the segments in an RP [55] is also used. The key measures of the quantification of the patterns present in the generated RP are detailed in Table 4. In addition, the values of these measures for the Henon chaotic attractor are shown to contrast the existence of chaos in the series.

We first consider the recurrence variable DET which gauges the proportion of recurrent points forming diagonal line structures and serves as a measure of determinism. For the Hénon attractor DET = 88.70, showing that most of the recurrent points present are found in deterministic structures. On the contrary, the value of DET in the four series analyzed is much lower than that of the Hénon attractor, with the GARCH series exhibiting the lowest value of 6.57. Another noteworthy recurrence variable is ENTR, representing the Shannon information entropy of all diagonal line lengths distributed over integer bins in a histogram. ENTR serves as a measure of signal complexity and reflects the complexity of the deterministic structure in the system. Across all series, its value is consistently lower than that associated with the Hénon attractor (2.56). Once again, the GARCH model residuals series records the lowest value at 0.14.

In summary, based on the RQA performed, it can be inferred that the electricity price series does not exhibit chaotic behavior.

In conclusion, the results obtained by applying all of the described methods show the presence of nonlinear dynamics, but they are not consistent with a chaotic structure. Therefore, it is concluded that traditional models that do not consider nonlinearity would not be suitable for the analysis and prediction of the electricity price time series. The results also reveal that GARCH/EGARCH models explain a significant part of the nonlinear structure. Therefore, it is necessary to investigate in the future whether other types of nonlinear models, both parametric and nonparametric, could explain the rest of the nonlinear dependence of the series.

4. Conclusions

This research explores the nature of the underlying dynamics for spot electricity prices in Spain. Using daily data spanning from 1 January 2013 to 29 March 2021 we investigated nonlinearity and the possible presence of chaos in the Spanish electricity market. Nonlinearity is scrutinized using the BDS method while chaotic dynamics are analyzed applying both metric (Lyapunov exponents) and topological (recurrence plots in combination with Recurrence Quantification Analysis) diagnostic tools. This approach relies on the fact that each method checks the different properties that characterize a chaotic behavior. The methods used to detect the presence of chaos encompass the best algorithms for analyzing noisy time series such as the one being investigated.

The findings show the presence of nonlinear dynamics but are not consistent with chaos. Thus, the empirical data propose that stochastic patterns, rather than deterministic principles, govern the system dynamics within the electricity spot market. Furthermore, GARCH models elucidate a substantial portion, although not the entirety, of the nonlinear structure apparent in the time series of electricity prices.

In a context where obtaining accurate forecasts is crucial for optimizing market performance, this investigation delivers fresh insights into the dynamics of electricity market prices. Regarding the benefits of our research for both theoretical comprehension and practical application, the outcomes substantially contribute to illustrating the behavior and underlying dynamics of the electricity price.

Similarly, the modelling and examination of the presence of nonlinearity and the chaotic dynamics of the spot electricity market are particularly useful for all stakeholders involved in the energy sector. This includes regulators, institutional agencies, producers, market forecasters, and policymakers.

Subsequent investigations should explore whether alternative nonlinear models could potentially capture the remaining nonlinear dependencies in the series. On the other hand, a study of the hourly price series, instead of the daily ones, may unveil patterns or behaviors that have not been observed in the daily series. It would also be of great interest to extend this research to other electricity markets and compare the outcomes obtained for each one. Additional work could be done by dividing the total sample into several periods according to the levels of volatility in the analyzed series.

Author Contributions

Conceptualization, L.I.-P. and S.G.y.G.; methodology, L.I.-P. and S.G.y.G.; software, L.I.-P. and S.G.y.G.; validation, L.I.-P. and S.G.y.G.; formal analysis, L.I.-P. and S.G.y.G.; investigation, L.I.-P. and S.G.y.G.; resources, L.I.-P. and S.G.y.G.; data curation, S.G.y.G.; writing—original draft preparation, L.I.-P. and S.G.y.G.; writing—review and editing, L.I.-P. and S.G.y.G.; visualization, L.I.-P. and S.G.y.G.; supervision, L.I.-P. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

The data presented in this study are available on request from the corresponding author.

Conflicts of Interest

The authors declare no conflicts of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1. Daily spot prices in the Spanish electricity market (SPEL) in EUR/MWh from 1 January 2013 to 31 March 2021. Time points (year) are on the x-axis and observations are on the y-axis.

Figure 2. Evolution of returns.

Figure 3. (A) Recurrence plot for the SPEL and returns series; (B) recurrence plot for the SARMA and GARCH series.

Table 1

Descriptive statistics.

	Original Series	Returns	SARMA	GARCH
Mean	45.9161	0.0002	0.0006	−0.0446
Median	47.3050	−0.0023	0.0140	0.0764
Maximum	94.9900	5.9376	4.1125	2.9347
Minimum	0.0100	−4.8675	−4.647	−9.7683
Std. Dev.	14.5141	0.3605	0.3099	0.9993
Skewness	−0.5139	0.1360	−2.6833	−2.2981
Kurtosis	3.7931	61.2181	61.1202	17.2129
Jarque–Bera	211.6 *	425,229.6 *	426,413.7 *	27,817.4 *
N	3012	3011	3004	2992

Note: Statistics for the considered time series: original series, return series, and residuals of the SARMA process, as well as standardized residuals from the GARCH model process. Skewness corresponds to Fisher’s skewness coefficient. Asterisk denotes that the Jarque-Bera test [58] result is significant at the 5 % confidence level. N is the number of observations. Source: Our own elaboration based on OMIP.

Table 2

Stationarity analysis.

Test	Original Series	Returns Series
Augmented Dickey Fuller Test	−3.197 (0.085)	−17.014 (0.000)
Phillips–Perron Test	−19.839 (0.000)	−82.278 (0.000)
Kwiatkowski–Phillips–Schmidt–Shin Test	0.389	0.004

Note: Constant and linear trend Augmented Dickey–Fuller p-values correspond to [59] one-sided p-values. For the Phillips–Perron test, lags were based on bandwidth Newey–West using Bartlett kernel. Critical value for the Kwiatkowski, Phillips, Schmidt, and Shin test is 0.146 for constant plus linear trend model. In bold are significant values at the 5% confidence level.

Table 3

Estimates of the maximum Lyapunov exponent.

	Lyapunov Exponent	p-Value *	Conclusion
Returns	−0.2739	0.000	No chaos
SARMA	−0.4283	0.000	No chaos
GARCH	−2.6099	0.000	No chaos

Note: *: At 5% significance level, the null hypothesis of the existence of chaos is rejected for those p-values below 0.05.

Table 4

Recurrence quantification measures.

Recurrence Measures	Original Series	Returns	SARMA Residuals	GARCH Residuals	Hénon Attractor (Chaotic System)
RR	3.76	4.58	2.06	2.79	1.21
DET	52.96	60.57	26.54	6.57	88.70
Lmax	2992	2976	2989	2977	N.A.
L	3.06	3.67	2.69	2.49	N.A.
ENTR	1.17	1.26	0.54	0.14	2.56
LAM	68.54	65.06	27.89	9.55	2.51
TT	3.53	3.88	2.29	2.07	2.00

Note: RR is the recurrence rate, DET stands for a measure of determinism, Lmax is the number of points forming the longest diagonal line, L is the mean value of the length of the diagonal lines, ENTR is the Shannon entropy, LAM is a measure of laminarity, and TT is the time the system spends trapped in a state. The values of the measures corresponding to the Hénon attractor have been taken from [57]. N.A.: Not available.

Appendix A

Figure A1. Periodogram (left) and autocorrelation function (right) for returns.

Table A1

SARMA(5,1)(1,2)₇ Model.

	Coefficient	Standard Deviation	t-Statistic	p-Value
C	3.72 × 10−6	4.73 × 10−5	0.0786	0.9374
AR(1)	0.6656	0.0101	66.2286	0.0000
AR(2)	−0.2171	0.0088	−24.5921	0.0000
AR(3)	0.1443	0.0094	15.4266	0.0000
AR(4)	0.0149	0.0079	1.8864	0.0593
AR(5)	−0.0779	0.0097	−8.0747	0.0000
SAR(7)	−0.7820	0.0384	−20.3332	0.0000
MA(1)	−0.8946	0.0097	−92.6031	0.0000
SMA(7)	−0.1267	0.0364	−3.4810	0.0005
SMA(14)	−0.8377	0.0353	−23.7352	0.0000

Table A2

GARCH(1,1) Model.

	Coefficient	Standard Deviation	t-Statistic	p-Value
C	0.0002	2.33 × 10−5	8.1965	0.0000
RESID(−1)2	0.2653	0.0093	28.5904	0.0000
GARCH(−1)	0.8121	0.0055	147.9187	0.0000

The linear model has been presented before.

Table A3

BDS test results for SARMA series.

ϵ	m
	0.5σ		σ		1.5σ		2σ
	BDS Statistic	p-Value	BDS Statistic	p-Value	BDS Statistic	p-Value	BDS Statistic	p-Value
2	0.051997	0.0000	0.054230	0.0000	0.039411	0.0000	0.028007	0.0000
3	0.080180	0.0000	0.105672	0.0000	0.083443	0.0000	0.061955	0.0000
4	0.089420	0.0000	0.148452	0.0000	0.126499	0.0000	0.097917	0.0000
5	0.087115	0.0000	0.180151	0.0000	0.164588	0.0000	0.132529	0.0000
6	0.079014	0.0000	0.203206	0.0000	0.197497	0.0000	0.164326	0.0000
7	0.069394	0.0000	0.220141	0.0000	0.225774	0.0000	0.193713	0.0000
8	0.060912	0.0000	0.232535	0.0000	0.250296	0.0000	0.220545	0.0000

Note: p-value below 0.05 is shown in bold. m stands for the embedding dimension and ϵ is the distance expressed by the number of standard deviations of the data: (0.5–2) · σ.

Table A4

BDS test results for GARCH series.

ϵ	m
	0.5σ		σ		1.5σ		2σ
	BDS Statistic	p-Value	BDS Statistic	p-Value	BDS Statistic	p-Value	BDS Statistic	p-Value
2	0.004098	0.0000	0.007652	0.0000	0.005329	0.0002	0.002470	0. 0134
3	0.003245	0.0000	0.009711	0.0000	0.008367	0.0012	0.004287	0.0236
4	0.001584	0.0014	0.007708	0.0022	0.008536	0.0098	0.005096	0.0672
5	0.000690	0.0078	0.004954	0.0252	0.006699	0.0846	0.004508	0.2014
6	0.000232	0.0500	0.002500	0.1608	0.003913	0.3392	0.003015	0.4610
7	6.85E-05	0.1892	0.001092	0.4202	0.001870	0.6316	0.002040	0.6500
8	3.46E-05	0.1626	0.000845	0.4262	0.001333	0.7146	0.001701	0.7138

Note: p-value below 0.05 is shown in bold. m stands for the embedding dimension and ϵ is the distance expressed by the number of standard deviations of the data: (0.5–2) · σ.

References

1. Weron, R. Electricity price forecasting: A review of the state-of-the-art with a look into the future. Int. J. Forecast.; 2014; 30, pp. 1030-1081. [DOI: https://dx.doi.org/10.1016/j.ijforecast.2014.08.008]

2. Tashpulatov, S.N. Modeling Electricity Price Dynamics Using Flexible Distributions. Mathematics; 2022; 10, 1757. [DOI: https://dx.doi.org/10.3390/math10101757]

3. Nowotarski, R.; Weron, R. Recent advances in electricity price forecasting: A review of probabilistic forecasting. Renew. Sust. Energy Rev.; 2018; 81, pp. 1548-1568. [DOI: https://dx.doi.org/10.1016/j.rser.2017.05.234]

4. Lago, J.; Marcjasz, G.; De Schutter, B.; Weron, R. Forecasting day-ahead electricity prices: A review of state-of-the-art algorithms, best practices and an open-access benchmark. Appl. Energy; 2021; 293, 116983. [DOI: https://dx.doi.org/10.1016/j.apenergy.2021.116983]

5. Lu, H.; Ma, X.; Ma, M.; Zhu, S. Energy price prediction using data-driven models: A decade review. Comput. Sci. Rev.; 2021; 39, 100356. [DOI: https://dx.doi.org/10.1016/j.cosrev.2020.100356]

6. Zema, T.; Sulich, A. Models of Electricity Price Forecasting: Bibliometric Research. Energies; 2022; 15, 5642. [DOI: https://dx.doi.org/10.3390/en15155642]

7. Chai, S.; Li, Q.; Abedin, M.Z.; Lucey, B.M. Forecasting electricity prices from the state-of-the-art modeling technology and the price determinant perspectives. Res. Int. Bus. Financ.; 2024; 67, 102132. [DOI: https://dx.doi.org/10.1016/j.ribaf.2023.102132]

8. Castro, A.L.; Marcato, A.L.M.; De Aguiar, E.P. Multifractal Analysis of the Brazilian Electricity Market. IEEE Access; 2023; 11, pp. 98939-98957. [DOI: https://dx.doi.org/10.1109/ACCESS.2023.3313099]

9. Facchini, A.; Rubino, A.; Caldarelli, G.; Di Liddo, G. Changes to Gate Closure and its impact on wholesale electricity prices: The case of the UK. Energy Policy; 2019; 125, pp. 110-121. [DOI: https://dx.doi.org/10.1016/j.enpol.2018.10.047]

10. Rassi, S.; Kanamura, T. Electricity price spike formation and LNG prices effect under gross bidding scheme in JEPX. Energy Policy; 2023; 177, 113552. [DOI: https://dx.doi.org/10.1016/j.enpol.2023.113552]

11. Avci-Surucu, E.; Aydogan, A.K.; Akgul, D. Bidding structure, market efficiency and persistence in a multi-time tariff setting. Energy Econ.; 2016; 54, pp. 77-87. [DOI: https://dx.doi.org/10.1016/j.eneco.2015.10.017]

12. Rypdal, M.; Løvsletten, O. Modeling electricity spot prices using mean-reverting multifractal processes. Phys. A; 2013; 392, pp. 94-207. [DOI: https://dx.doi.org/10.1016/j.physa.2012.08.004]

13. Erzgräber, H.; Strozzi, F.; Zaldívar, J.M.; Touchette, H.; Gutiérrez, E.; Arrowsmith, D.K. Time series analysis and long range correlations of Nordic spot electricity market data. Phys. A; 2008; 387, pp. 6567-6574. [DOI: https://dx.doi.org/10.1016/j.physa.2008.07.030]

14. Strozzi, F.; Zaldívar, J.M.; Zbilut, J.P. Application of nonlinear time series analysis techniques to high-frequency currency exchange data. Phys. A; 2002; 312, pp. 520-538. [DOI: https://dx.doi.org/10.1016/S0378-4371(02)00846-4]

15. Papaioannou, G.P.; Dikaiakos, C.; Stratigakos, A.C.; Papageorgiou, P.C.; Krommydas, K.F. Testing the Efficiency of Electricity Markets Using a New Composite Measure Based on Nonlinear TS Tools. Energies; 2019; 12, 618. [DOI: https://dx.doi.org/10.3390/en12040618]

16. Yousefpoor, P.; Esfahani, M.S.; Nojumi, H. Looking for systematic approach to select chaos tests. Appl. Math. Comput.; 2008; 198, pp. 73-91. [DOI: https://dx.doi.org/10.1016/j.amc.2007.08.070]

17. Lahmiri, S. Investigating existence of chaos in short and long term dynamics of Moroccan exchange rates. Phys. A; 2017; 465, pp. 655-661. [DOI: https://dx.doi.org/10.1016/j.physa.2016.08.024]

18. Hsieh, D.A. Chaos and nonlinear dynamics: Application to financial markets. J. Financ.; 1991; 46, pp. 1839-1877. [DOI: https://dx.doi.org/10.1111/j.1540-6261.1991.tb04646.x]

19. Faggini, M. Chaotic time series analysis in economics: Balance and perspectives. Chaos; 2014; 24, 042101. [DOI: https://dx.doi.org/10.1063/1.4903797]

20. Bask, M.; Tung, L.; Widerberg, A. The stability of electricity prices: Estimation and inference of the Lyapunov exponents. Phys. A; 2007; 376, pp. 565-572. [DOI: https://dx.doi.org/10.1016/j.physa.2006.10.016]

21. Bigdeli, N.; Afshar, K. Chaotic behavior of price in the power markets with pay-as-bid payment mechanism. Chaos Solitons Fract.; 2009; 42, pp. 2560-2569. [DOI: https://dx.doi.org/10.1016/j.chaos.2009.03.193]

22. Brock, W.A.; Dechert, W.D.; Sheinkman, J.A. A Test of Independence Based on the Correlation Dimension; SSRI no. 8702 Department of Economics, University of Wisconsin: Madison, WI, USA, 1987.

23. Brock, W.A.; Scheinkman, J.A.; Dechert, W.D.; LeBaron, B. A test for independence based on the correlation dimension. Econ. Rev.; 1996; 15, pp. 197-235. [DOI: https://dx.doi.org/10.1080/07474939608800353]

24. BenSaïda, A.; Litimi, H. High level chaos in the exchange and index markets. Chaos Solitons Fract.; 2013; 54, pp. 90-95. [DOI: https://dx.doi.org/10.1016/j.chaos.2013.06.004]

25. Inglada-Pérez, L.; Coto-Millán, P. A Chaos Analysis of the Dry Bulk Shipping Market. Mathematics; 2021; 9, 2065. [DOI: https://dx.doi.org/10.3390/math9172065]

26. Inglada-Perez, L. A Comprehensive Framework for Uncovering Non-Linearity and Chaos in Financial Markets: Empirical Evidence for Four Major Stock Market Indices. Entropy; 2020; 22, 1435. [DOI: https://dx.doi.org/10.3390/e22121435] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/33353243]

27. Box, G.E.P.; Jenkins, G.M. Time Series Analysis: Forecasting and Control; Holden-Day: San Francisco, CA, USA, 1970.

28. Engle, R.F.; Bollerslev, T. Modelling the persistence of conditional variances. Econom. Rev.; 1986; 5, pp. 1-50. [DOI: https://dx.doi.org/10.1080/07474938608800095]

29. Dickey, D.; Fuller, W. Distribution of the estimators for autoregressive time series with a unit root. J. Am. Stat. Assoc.; 1979; 74, pp. 426-431.

30. Phillips, P.C.B.; Perron, P. Testing for Unit Roots in Time Series Regression. Biometrika; 1988; 75, pp. 335-346. [DOI: https://dx.doi.org/10.1093/biomet/75.2.335]

31. Kwiatkowski, D.; Phillips, P.C.B.; Schmidt, P.; Shin, Y. Testing the Null Hypothesis of Stationarity against the Alternative of a Unit Root. J. Econom.; 1992; 54, pp. 159-178. [DOI: https://dx.doi.org/10.1016/0304-4076(92)90104-Y]

32. Schwarz, G. Estimating the dimension of a model. Ann. Stat.; 1978; 6, pp. 461-464. [DOI: https://dx.doi.org/10.1214/aos/1176344136]

33. Engle, R.F. Autoregressive conditional heteroscedasticity with estimates of the variance of United Kingdom inflation. Econometrica; 1982; 50, pp. 987-1007. [DOI: https://dx.doi.org/10.2307/1912773]

34. Eckmann, J.; Ruelle, D. Ergodic theory of chaos and strange attractors. Rev. Mod. Phys.; 1985; 57, pp. 617-656. [DOI: https://dx.doi.org/10.1103/RevModPhys.57.617]

35. Takens, F. Detecting strange attractors in turbulence. Dynamical Systems and Turbulence; Rand, D.; Young, L. Springer: Berlin/Heidelberg, Germany, 1981; pp. 366-381.

36. Brock, W.A.; Hsieh, D.A.; LeBaron, B.D. Nonlinear Dynamics, Chaos, and Instability: Statistical Theory and Economic Evidence; MIT Press: Cambridge, MA, USA, 1991.

37. Barnett, W.A.; Gallant, A.R.; Hinich, M.J.; Jungeilges, J.A.; Kaplan, D.T.; Jensen, M.J. A single-blind controlled competition among tests for nonlinearity and chaos. J. Econom.; 1997; 82, pp. 157-192. [DOI: https://dx.doi.org/10.1016/S0304-4076(97)00081-X]

38. Kim, H.; Eykholt, R.; Salas, J. Nonlinear dynamics, delay times, and embedding windows. Phys. D; 1999; 127, pp. 48-60. [DOI: https://dx.doi.org/10.1016/S0167-2789(98)00240-1]

39. Grassberger, P.; Procaccia, I. Dimensions and entropies of strange attractors from a fluctuating dynamics approach. Physica; 1984; 13, pp. 34-54. [DOI: https://dx.doi.org/10.1016/0167-2789(84)90269-0]

40. Lahmiri, S. On fractality and chaos in Moroccan family business stock returns and volatility. Phys. A; 2017; 473, pp. 29-39. [DOI: https://dx.doi.org/10.1016/j.physa.2017.01.033]

41. Nychka, D.; Ellner, S.; Gallant, A.R.; McCaffrey, D. Finding chaos in noisy systems. J. R. Stat. Soc. Series B Stat. Methodol.; 1992; 54, pp. 399-426. [DOI: https://dx.doi.org/10.1111/j.2517-6161.1992.tb01889.x]

42. Hernández-Gómez, J.J.; Couder-Castañeda, C.; Herrera-Díaz, I.E.; Flores-Guzmán, N.; Gómez-Cruz, E. Weakly Coupled Distributed Calculation of Lyapunov Exponents for Non-Linear Dynamical Systems. Algorithms; 2017; 10, 137. [DOI: https://dx.doi.org/10.3390/a10040137]

43. Wolf, A.; Swift, J.B.; Swinney, H.L.; Vastano, J.A. Determining Lyapunov exponents from a time series. Phys. D; 1985; 16, pp. 285-317. [DOI: https://dx.doi.org/10.1016/0167-2789(85)90011-9]

44. Shintani, M.; Linton, O. Nonparametric neural network estimation of Lyapunov exponents and direct test for chaos. J. Econom.; 2004; 120, pp. 1-33. [DOI: https://dx.doi.org/10.1016/S0304-4076(03)00205-7]

45. McCaffrey, D.; Ellner, S.; Gallant, A.; Nychka, D. Estimating the Lyapunov exponent of a chaotic system with nonparametric regression. J. Am. Stat. Assoc.; 1992; 87, pp. 682-695. [DOI: https://dx.doi.org/10.1080/01621459.1992.10475270]

46. Bailey, B.A.; Ellner, S.; Nychka, D.W. Chaos with confidence: Asymptotics and applications of local Lyapunov exponents. Nonlinear Dynamics and Time Series: Building a Bridge Between the Natural and Statistical Sciences; Cutler, C.D.; Kaplan, D.T. Fields Institute Communications, American Mathematical Society: Providence, RI, USA, 1997; pp. 115-133.

47. Hornik, K.; Stinchcombe, M.; White, H. Multilayer feedforward networks are universal approximators. Neural Netw.; 1989; 2, pp. 359-366. [DOI: https://dx.doi.org/10.1016/0893-6080(89)90020-8]

48. Sandubete, J.E.; Escot, L. DChaos: An R Package for Chaotic Time Series Analysis. R J.; 2021; 13, pp. 232-252. [DOI: https://dx.doi.org/10.32614/RJ-2021-036]

49. Eckmann, J.P.; Kamphorst, S.O.; Ruelle, D. Recurrence plots of dynamical systems. Europhys. Lett.; 1987; 4, pp. 973-977. [DOI: https://dx.doi.org/10.1209/0295-5075/4/9/004]

50. Marwan, N.; Romano, M.C.; Thiel, M.; Kurths, J. Recurrence plots for the analysis of complex systems. Phys. Rep.; 2007; 438, pp. 237-329. [DOI: https://dx.doi.org/10.1016/j.physrep.2006.11.001]

51. Inglada-Perez, L. Uncovering nonlinear dynamics in air transport demand. Int. J. Transp. Econ.; 2016; XLIII, pp. 33-66.

52. Fraser, A.M.; Swinney, H.L. Independent Coordinates for Strange Attractors from Mutual Information. Phys. Rev. A; 1986; 33, pp. 1134-1140. [DOI: https://dx.doi.org/10.1103/PhysRevA.33.1134] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/9896728]

53. Kennel, M.B.; Brown, R.; Mand Abarbanel, H.D. Determining embedding dimension for phase-space reconstruction using a geometrical construction. Phys. Rev. A; 1992; 45, 3403. [DOI: https://dx.doi.org/10.1103/PhysRevA.45.3403] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/9907388]

54. Webber, C.L., Jr.; Zbilut, J.P. Dynamical assessment of physiological systems and states using recurrence plot strategies. J. Appl. Physiol.; 1994; 76, pp. 965-973. [DOI: https://dx.doi.org/10.1152/jappl.1994.76.2.965] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/8175612]

55. Trulla, L.; Giuliani, A.; Zbilut, J.; Webber, C.L., Jr. Recurrence quantification analysis of the logistic equation with transients. Phys. Lett. A; 1996; 223, pp. 255-260. [DOI: https://dx.doi.org/10.1016/S0375-9601(96)00741-4]

56. Marwan, N.; Wessel, N.; Meyerfeldt, U.; Schirdewan, A.; Kurths, J. Recurrence- plot-based measures of complexity and their application to heart-rate-variability data. Phys. Rev. E; 2002; 66, 026702. [DOI: https://dx.doi.org/10.1103/PhysRevE.66.026702] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/12241313]

57. Shannon, C.E. A mathematical theory of communication. Bell Syst. Tech. J.; 1948; 27, pp. 379-423. [DOI: https://dx.doi.org/10.1002/j.1538-7305.1948.tb01338.x]

58. Jarque, C.M.; Bera, A.K. A test for normality of observations and regression residuals. Int. Stat. Rev.; 1987; 65, pp. 163-172. [DOI: https://dx.doi.org/10.2307/1403192]

59. MacKinnon, J.G. Numerical distribution functions for unit root and cointegration tests. J. Appl. Econom.; 1996; 11, pp. 601-618. [DOI: https://dx.doi.org/10.1002/(SICI)1099-1255(199611)11:6<601::AID-JAE417>3.0.CO;2-T]

60. Brock, W.A.; Potter, S. Nonlinear Time Series and Macroeconomics. Handbook of Statistics; Maddala, G.S.; Rao, C.R.; Vinod, H.D. Elsevier: Amsterdam, The Netherlands, 1993.

61. Zbilut, P.; Zaldivar-Comenges, J.M.; Strozzi, F. Recurrence quantification based Lyapunov exponents for monitoring divergence in experimental data. Phys. Lett. A; 2002; 297, pp. 173-181. [DOI: https://dx.doi.org/10.1016/S0375-9601(02)00436-X]