Short-Term Electricity Prices Forecasting Using Functional Time Series Analysis

Abstract

In recent years, efficient modeling and forecasting of electricity prices became highly important for all the market participants for developing bidding strategies and making investment decisions. However, as electricity prices exhibit specific features, such as periods of high volatility, seasonal patterns, calendar effects, nonlinearity, etc., their accurate forecasting is challenging. This study proposes a functional forecasting method for the accurate forecasting of electricity prices. A functional autoregressive model of order $\mathbb{P}$ is suggested for short-term price forecasting in the electricity markets. The applicability of the model is improved with the help of functional final prediction error (FFPE), through which the model dimensionality and lag structure were selected automatically. An application of the suggested algorithm was evaluated on the Italian electricity market (IPEX). The out-of-sample forecasted results indicate that the proposed method performs relatively better than the nonfunctional forecasting techniques such as autoregressive (AR) and naïve models.

1. Introduction

In the late 1980s, the worldwide electricity industry had undergone numerous fundamental changes when the state-owned monopolistic structure was restructured into the deregulated and competitive electricity market. The main driving force behind the restructuring of the electricity market was to promote competition among producers, retailers, and consumers by boosting private investments in production, supply, and retail sectors. Liberalization of this sector brought many benefits to the stakeholders in terms of reliable, secure, and economical electricity trading. However, due to electricity’s inherent physical characteristic of non-storability in large volumes, the uncertainty related to electricity prices and demand forecasting increased. In addition, electricity prices and demand series generally exhibit specific features, such as multiple periodicities, long-trend, bank holiday effect, spikes, jumps, etc. In the presence of these features, the forecasting problem is challenging in all three forecasting horizons, i.e., short term, medium term, and long term [1].

In electricity markets, short-term forecasting refers to forecasting electricity prices from a few minutes to a week ahead. Apart from the power scheduling, management, and risk assessment, a short-term forecast is essential for market participants to optimize their bidding strategies. Medium-term forecast generally refers to the forecast made for a few weeks to a few months ahead. It is usually vital for expanding generation plants, scheduling maintenance, developing investment, fuel contracting, bilateral contracting, and hedging strategies. Forecasts ranging from a few months ahead to a few years ahead are commonly referred to as long-term-ahead forecasts. They are used for planning and investment profitability analysis, i.e., making decisions for future investments in power plants, inducing sites, and fuel sources [2,3]. In the literature, short-term forecasting has received greater research attention as the maximum electricity trade takes place in this market.

The literature concerning electricity price forecasting reported several statistical, machine learning, econometric, and hybrid models used to forecast short-term electricity prices [4,5,6,7]. Different linear time series models, including AR, ARMA, ARIMA, SARIMA, and ARIMAX [8,9,10,11,12], and nonlinear time series models, such as NPAR, ARCH, GARCH, and their extensions [13,14,15], are extensively used for forecasting electricity prices. Parametric and nonparametric regression-type models considering multiple, local polynomial, kernel, smoothing spline, and quantile regression are easy to implement and are widely studied in the case of electricity price forecasting [2,16,17,18,19,20,21]. In addition, models based on exponential smoothing including simple, double, and triple Holt’s winters models that account for various periodicities [22,23,24,25,26,27] are often used for forecasting purposes. Artificial intelligence models have also been used to predict day-ahead electricity prices [28,29,30,31,32], as well as state-space models [33,34]. Various researchers combined the characteristics of two or more models to build a new model generally referred to as a hybrid model [3,35,36,37,38,39]. Generally, the above-stated models have their own functional and structural form, and the forecasting performance varies from market to market [40].

In the last three decades, technological developments simplified and decreased the cost of data collection and storage processes. Such advancements helped us to examine and record practical life activities in great detail. Examples include curves, images, surfaces, or anything else varying over a continuum. Consequently, classical statistical analysis techniques are inadequate and inefficient due to the large dimensions of data. To analyze such datasets, some suitable statistical methods are required, and functional data analysis (FDA) is one of the prominent methods to tackle such data in an efficient way. The FDA presents the essential statistical background for the analysis of functional variables, where every observation is a continuous function. The application of the FDA exists in almost every field of science, including economics, environment, engineering, energy, etc. [41,42]. In this research work, the application of the FDA is proposed for the electricity market, which is of primary interest for many researchers working in this field, especially after the liberalization of this market.

Given the temporal dependence, the FAR models have been suggested for the time series of trajectories. The autoregressive Hilbertian (ARH) process proposed by [43], also called the FAR model under Hilbert space, is likely the most popular pioneering work that plays an important role in the FDA context.The FAR is an extension of the AR process to infinite-dimensional space and is also used in electricity price forecasting. For example, using functional analysis of variance (FANOVA) and FAR model, Ref. [44] studied the seasonal patterns and improved prediction accuracy for electricity demand time series used from the Nord Pool electricity market. The application of a local linear method with functional explanatory variables was studied by [45]. They compared their proposed approach with the functional Nadaraya–Watson (NW) method and other finite-dimensional nonparametric techniques. For empirical analysis, monthly electricity consumption data of the United States of America (USA) were used, and the results suggest the superior performance of their proposed methods. The forecasting performances of different parametric and nonparametric functional models for electricity demand were studied by [46]. The authors used data from the Italian and British electricity markets and concluded that the nonparametric functional models give superior performance to their parametric counterparts. In another study, Ref. [47] used different functional models and compared their results with the finite univariate dimension (univariate and multivariate) models. Data from four different electricity markets, namely, the Nord Pool electricity market (NP), Pennsylvania–New Jersey–Maryland electricity market (PJM), the Italian electricity market (IPEX), and the British electricity market (APX Power UK), were used, and the results were summarized using different descriptive measures. The results suggested that the functional approach produces better results than the rest. Ref. [48] used the electricity demand curves data from Southern Australia. The author sliced the univariate time series into curves and reduced their dimensionality by applying the functional principal components technique. Finally, the author used univariate time series models to predict short-term electricity demand.

The main aim of this research work is to propose a functional model that can efficiently predict electricity prices. To this end, a method based on a two-components estimation procedure is proposed. The first component, known as the deterministic component, is computed using the additive modeling technique. The stochastic component, on the other hand, is modeled using an FAR($\mathbb{P}$) model where the selection of the dimension and lags is automatic. Finally, the model is tested for a whole year to see its forecasting performance. The rest of this paper is organized as follows. Section 2 provides an overview of the preliminaries. Section 3 describes a comprehensive review of the FAR($\mathbb{P}$) and functional final prediction error (FFPE). Section 4 provides the application of the proposed method, while Section 5 concludes the study.

2. Functional Modeling

2.1. Preliminaries

Let $\{{\mathbb{Z}}_i\left(t\right):i\in \mathbb{N},t\in \mathcal{J}\}$ be an arbitrary stationary $\mathbb{N}$-dimensional time series where $\mathcal{J}$ represents a continuum bounded within a finite interval. For each i, the functional observation ${\mathbb{Z}}_i$ belongs to a Hilbert space $\mathbb{H}={\mathcal{L}}^2([0,1],\parallel ·\parallel )$ of square integrable functions which is equipped with a norm $\parallel ·\parallel$ induced by the inner product $<g,h>=\int g\left(t\right)h\left(t\right)dt$. The object $\left\{{\mathbb{Z}}_i\left(t\right)\right\}$ is referred to as FTS with i as the time index [49,50]. Furthermore, all stochastic functions are defined on a common probability space $(\mathsf{\Omega },\mathcal{A},P)$. The notation $\mathbb{Z}\in {\mathcal{L}}_{\mathbb{H}}^p(\mathsf{\Omega },\mathcal{A},P)$ is used to indicate $\mathbb{E}(\parallel \mathbb{Z}{\parallel }^p)<\infty$ for some $p>0$. When $p=1$, $\mathbb{Z}\left(t\right)$ has the mean curve $\mu \left(t\right)$; when $p=2$, the covariance operators $C(t,s)$ are defined as in Equations (1) and (2) as under

$$\mu \left(t\right)=\mathbb{E}\left[\mathbb{Z}\right(t\left)\right]$$

(1)

$$C(t,s)=\mathbb{E}\left[\right(\mathbb{Z}\left(t\right)-\mu \left(t\right)\left)\right(\mathbb{Z}\left(s\right)-\mu \left(s\right)\left)\right]$$

(2)

Mercer’s theorem [51] provides the following convenient spectral decomposition of Equation (2):

$$C(t,s)=\sum _{j=1}^{\infty }{\kappa }_j{\phi }_j\left(t\right){\phi }_j\left(s\right)$$

(3)

where ${\phi }_j$ denotes the jth orthonormal principal component, and ${\kappa }_j$ denotes the jth eigenvalue. The principal component scores (PCSs) ${\gamma }_{i,j}$ are given by the projection of $[{\mathbb{Z}}_i\left(t\right)-\mu \left(t\right)]$ in the direction of the jth eigenfunction ${\phi }_j$, i.e., ${\gamma }_{i,j}=〈{\mathbb{Z}}_i-\mu ,{\phi }_j〉$. Based on the separability of the Hilbert space, the Karhunen–Loève (KL) expansion [52,53] of the random function $\mathbb{Z}\left(t\right)$ can be expressed as

$$\mathbb{Z}\left(t\right)=\mu \left(t\right)+\sum _{j=1}^{\infty }{\gamma }_{i,j}{\phi }_j\left(t\right)$$

(4)

The KL expansion provides the theoretical background for FPCA; see [54,55] for more details about FPCA and its practical demonstration.

Expansion (4) facilitates dimension reduction as the first $\mathbb{D}$ terms often provide a good approximation to the infinite sums, and, thus, the information contained in $\mathbb{Z}\left(t\right)$ can be adequately summarized by the jth-dimensional vector (${\gamma }_1,\dots ,{\gamma }_j$). The approximated processes can be defined as

$$\mathbb{Z}\left(t\right)=\mu \left(t\right)+\sum _{j=1}^{\mathbb{D}}{\gamma }_j{\phi }_j\left(t\right)+ϵ\left(t\right)$$

(5)

where $ϵ\left(t\right)$ denotes the zero-mean white noise function that captures the variation excluded from the first $\mathbb{D}$ leading functional principal components (FPCs). There are different methods available in the literature for choosing the value of $\mathbb{D}$: (i) scree plots or the fraction of variation explained by first few PCSs [56], (ii) using the Akaike information and Bayesian information criteria [57], (iii) cross-validation with one-curve-leave-out or k-fold method [58], or (iv) bootstrap techniques [59].

Once the sample functional data are available, the sample mean can be obtained as

$$\widehat{\mu }\left(t\right)=\frac{1}{\mathbb{N}}\sum _{i=1}^{\mathbb{N}}{\mathbb{Z}}_i\left(t\right),t\in [0,1],$$

(6)

and the sample covariance function is defined as

$$\widehat{C}(t,s)=\frac{1}{\mathbb{N}-1}\sum _{i=1}^{\mathbb{N}}\left({\mathbb{Z}}_i\left(t\right)-\widehat{\mu }\left(t\right)\right)\left({\mathbb{Z}}_i\left(s\right)-\widehat{\mu }\left(s\right)\right)$$

(7)

Ref. [60] proved that the estimators are consistent for weakly dependent process.

2.2. Functional Autoregressive Model

Autoregressive (AR) models are one of the most popular forecasting models used in time series analysis. In the AR modeling framework, the response variable is linearly dependent on it past p lags with an error term. The theory of AR and more general linear processes in Hilbert spaces is developed in the monograph of [50], containing sufficient technical details. In addition, more relevant information can also be found in [49,61].

Recall a sequence of stationary random curves (${\mathbb{Z}}_i\left(t\right),i\in \mathcal{N}$) in ${\mathcal{L}}^2\left([0,1]\right)$ defined in Section 2.1. The functional AR model of order $\mathbb{P}$ (FAR $\left(\mathbb{P}\right)$) can be written as [50]:

$${\mathbb{Z}}_i\left(t\right)-\mu \left(t\right)=\sum _{k=1}^{\mathbb{P}}{\mathsf{\Psi }}_k\left({\mathbb{Z}}_{i-k}\left(t\right)-\mu \left(t\right)\right)+{\xi }_i\left(t\right)$$

(8)

where ${\mathsf{\Psi }}_k(k=1,\dots ,\mathbb{P})$ are the FAR operators (functional parameters), $\mu \left(t\right)$ is the mean function of ${\mathbb{Z}}_i\left(t\right)$, ${\mathbb{Z}}_{i-k}\left(t\right)$ denotes $k\mathrm{th}$ lag of curve ${\mathbb{Z}}_i$, and ${\xi }_i\left(t\right)$ is a strong $\mathbb{H}$-white noise with zero mean and finite second moment ($\mathbb{E}{∥{\xi }_i\left(t\right)∥}^2<\infty$). For the prediction and forecasting of the model given in Equation (8), the following forecasting algorithm is used, which is based on Equations (5)–(7) [62].

• First, the dimension which is denoted by $\mathbb{D}$ is fixed by using the method described in Section 2.3, and the estimated FPC scores are obtained as ${\widehat{\gamma }}_{i,j}=\int {\widehat{\mathbb{Z}}}_i\left(t\right){\widehat{\phi }}_j\left(t\right)dt$ for each observation ${\widehat{\mathbb{Z}}}_i\left(t\right),i=1,\dots ,\mathbb{N},j=1,\dots ,\mathbb{D}$, and the estimated j-variate FPC scores vectors ${\widehat{\mathit{\gamma }}}_i={({\widehat{\gamma }}_{1,i},\dots ,{\widehat{\gamma }}_{\mathbb{D},i})}^t,i=1,\dots ,\mathbb{N}$.
• Next, the order $\mathbb{P}$ is fixed using the technique described in Section 2.3 and we fit the vector AR model, VAR($\mathbb{P}$), as ${\mathit{\gamma }}_i={\sum }_{k=1}^{\mathbb{P}}{\mathsf{\Psi }}_k{\mathit{\gamma }}_{i-k}+{ϵ}_i$ for eigenscores vectors to produce forecasting ${\widehat{\mathit{\gamma }}}_{\mathbb{N}+1}={({\widehat{\gamma }}_{\mathbb{N}+1,1},\dots ,{\widehat{\gamma }}_{\mathbb{N}+1,\mathbb{D}})}^t$. Durbin–Levinson and innovations algorithm can be readily applied here, given the vectors ${\widehat{\mathit{\gamma }}}_1,\dots ,{\widehat{\mathit{\gamma }}}_{\mathbb{N}}$.
• In the last step, the multivariate time series are converted back to functional version using the KL theorem ${\widehat{\mathbb{Z}}}_{\mathbb{N}+1}\left(t\right)=\widehat{\mu }\left(t\right)+{\widehat{\gamma }}_{\mathbb{N}+1,1}{\widehat{\phi }}_1\left(t\right)+\cdots +{\widehat{\gamma }}_{\mathbb{N}+1,\mathbb{D}}{\widehat{\phi }}_{\mathbb{D}}\left(t\right)$. The FPC scores and sample eigenfunctions result in ${\widehat{\mathbb{Z}}}_{\mathbb{N}+1}\left(t\right)$, which is then used as a one-step-ahead forecast of ${\mathbb{Z}}_{\mathbb{N}+1}\left(t\right)$.

As can be seen, the selection of the dimension $\mathbb{D}$ and lags $\mathbb{P}$ is an important step in the above algorithm. The following section illustrates how to select the optimal values for these variables.

2.3. Selection of Order and Dimension of FAR($\mathbb{P}$)

The main goal of the current article is the accurate forecasting through FAR($\mathbb{P}$), which requires the appropriate order $\mathbb{P}$ selection as well as the dimension $\mathbb{D}$, in such a way that the mean square error (MSE) is minimized.

As the eigenfunctions ${\phi }_j$ and the PCS’s ${\gamma }_{\mathbb{N},j}$ are uncorrelated, the MSE can be decomposed as

$$\begin{array}{ccc}\hfill \mathbb{E}\left\{{∥{\mathbb{Z}}_{\mathbb{N}+1}-{\widehat{\mathbb{Z}}}_{\mathbb{N}+1}∥}^2\right\}& =& \mathbb{E}\left\{{∥\sum _{j=1}^{\infty }{\gamma }_{\mathbb{N}+1,j}{\phi }_j-\sum _{j=1}^{\mathbb{D}}{\widehat{\gamma }}_{\mathbb{N}+1,j}{\phi }_j∥}^2\right\}\hfill \\ & =& \mathbb{E}\left\{{∥{\mathbf{Z}}_{\mathbb{N}+\mathbb{1}}-{\widehat{\mathbf{Z}}}_{\mathbb{N}+1}∥}^2\right\}+\sum _{j=\mathbb{D}+1}^{\infty }{\kappa }_j\hfill \end{array}$$

where ${∥.∥}^2$ denotes the usual l-2 Euclidean norm of vectors. We suppose that the vector ${\mathbf{Z}}_{\mathbb{N}}$ is stationary and follows a $\mathbb{D}$-variables vector AR of order $\mathbb{P}$, VAR($\mathbb{P}$), that can be written as

$${\mathit{Z}}_{\mathbb{N}+1}={\mathsf{\Phi }}_1{\mathit{Z}}_{\mathbb{N}}+{\mathsf{\Phi }}_2{\mathit{Z}}_{\mathbb{N}-1}+,\dots ,+{\mathsf{\Phi }}_{\mathbb{P}}{\mathit{Z}}_{\mathbb{N}-\mathbb{P}+1}+{\mathit{Y}}_{\mathbb{N}+1}.$$

(9)

Ref. [63] showed that $\left({\mathit{Y}}_{\mathbb{N}}\right)$ is a white noise process such that

$$\sqrt{\mathbb{N}}(\widehat{\rho }-\rho )\stackrel{D}{\to }\mathbf{\mathbb{N}}(\mathbf{0},{\Sigma }_{\mathit{Y}}\otimes {\mathsf{\Delta }}_{\mathbb{P}}^{-1})$$

(10)

where $\rho$ = vec ${[{\mathsf{\Phi }}_1,\dots ,{\mathsf{\Phi }}_{\mathbb{P}}]}^t$ and $\widehat{\rho }$ = vec ${[{\widehat{\mathsf{\Phi }}}_1,\dots ,{\widehat{\mathsf{\Phi }}}_{\mathbb{P}}]}^t$ is the least squares estimator in vector form, and ${\mathsf{\Delta }}_{\mathbb{P}}$ = var[vec $({\mathit{Z}}_{\mathbb{P}},\dots ,{\mathit{Z}}_1)]$ and ${\Sigma }_{\mathit{Y}}=\mathbb{E}[{\mathit{Y}}_1,{\mathit{Y}}_1^t]$. Assume that the $\widehat{\rho }$ are estimated from independent training sample $({\mathit{X}}_1,\dots ,{\mathit{X}}_{\mathbb{N}})\stackrel{\mathbb{D}}{=}({\mathit{Z}}_1,\dots ,{\mathit{Z}}_{\mathbb{N}})$. It follows then that

$$\begin{array}{ccc}\hfill \mathbb{E}\left\{{∥{\mathit{Z}}_{\mathbb{N}+1}-{\widehat{\mathit{Z}}}_{\mathbb{N}+1}∥}^2\right\}& =& \mathbb{E}\left\{{∥{\mathit{Z}}_{\mathbb{N}+1}-({\widehat{\mathsf{\Phi }}}_1{\mathit{Z}}_{\mathbb{N}}+\cdots +{\widehat{\mathsf{\Phi }}}_{\mathbb{P}}{\mathit{Z}}_{\mathbb{N}-\mathbb{P}+1})∥}^2\right\}\hfill \\ & =& \mathbb{E}\left\{{∥{\mathit{Y}}_{\mathbb{N}+1}∥}^2\right\}+\mathbb{E}\left\{{∥({\mathsf{\Phi }}_1-{\widehat{\mathsf{\Phi }}}_1){\mathit{Z}}_{\mathbb{N}}+\cdots +({\mathsf{\Phi }}_{\mathbb{P}}-{\widehat{\mathsf{\Phi }}}_{\mathbb{P}}){\mathit{Z}}_{\mathbb{N}-\mathbb{P}+1}∥}^2\right\}\hfill \\ & =& \mathrm{trace}\left\{{\Sigma }_{\mathit{Y}}\right\}+\mathbb{E}\left\{{∥I_{\mathbb{P}}\otimes ({\mathit{Z}}_{\mathbb{N}}^t,\dots ,{\mathit{Z}}_{\mathbb{N}-\mathbb{P}+1}^t)(\rho -\widehat{\rho })∥}^2\right\}\hfill \end{array}$$

(11)

For some further derivation by using Equation (10), Ref. [64] showed that Equation (11) can be approximated as

$\mathbb{E}{∥{\mathit{Z}}_{\mathbb{N}+1}-{\widehat{\mathit{Z}}}_{\mathbb{N}+1}∥}^2\approx \frac{\mathbb{N}+\mathbb{P}\ast \mathbb{D}}{\mathbb{N}-\mathbb{P}\ast \mathbb{D}}$ trace $\left({\widehat{\Sigma }}_{\mathit{Y}}\right)+{\sum }_{j>\mathbb{D}}{\kappa }_j$. The suggested functional final prediction error selects order $\mathbb{P}$ and dimension $\mathbb{D}$ simultaneously by minimizing error term.

$$fFPE(\mathbb{P},\mathbb{D})=\frac{\mathbb{N}+\mathbb{P}\ast \mathbb{D}}{\mathbb{N}-\mathbb{P}\ast \mathbb{D}}\mathrm{trace}\left({\widehat{\Sigma }}_{\mathit{Y}}\right)+\sum _{j>\mathbb{D}}{\kappa }_j$$

(12)

Using the fFPE method, the suggested forecasting procedure works in a completely data-driven-based way and does not require any subjective specification of parameters. It is specifically important that the choice of $\mathbb{D}$ depends upon the sample size $\mathbb{N}$. For more technical details, the interested readers are referred to [64] and the references cited therein.

3. Modeling Framework

This section provides the general modeling framework used to model and forecast electricity prices. As described in Section 1, electricity prices exhibit specific features, e.g., extreme values (outliers), multiple periodicities, bank holidays effect, etc. Incorporating these specific features in the model greatly improves the forecasting accuracy [47]. To this end, the price time series is first filtered using the moving window filter on prices discussed in the following section.

3.1. Moving Window Filter on Prices

The identification of outliers, also known as the extreme values, in the data is one of the growing research areas. Various methods and ideas have been used in the literature to detect and impute outliers in the data. The significant developments in terms of outliers detection techniques in time series are suggested by [65,66,67,68]. Generally, the presence of outliers in the original electricity price data can substantially influence most forecasting models, which can result in poor forecasting performance. Therefore, identifying and analyzing outliers in the data is an essential step in constructing a forecasting model.

The moving window filter on price (MFP) [69] is an extension of the standard deviation filter on prices (SFP) technique. The SFP technique is based on the idea that the prices whose absolute deviation is taken from the mean $\widehat{\mu }$ and are greater than some multiple of the sample standard deviation $\widehat{\sigma }$ are referred to as outliers. However, the MFP technique differs from the SFP in the sense that it works out with the rolling window having fixed width of intervals. Using the MFP technique, the original price series is divided into $\mathbb{N}=\mathbb{T}/\mathbb{M}$ parts, where $\mathbb{M}$ is the width of the windows. Then, the SFP technique is applied to the first window of the given time series. Next, the window is shifted into the next fixed interval of $\mathbb{M}$ width, and the SFP is applied. Finally, the process is repeated until the last window is treated. Our work considers the same predictive interval used in [69], with the width of the window being equal to ten weeks. Thus, the subset of outliers ${\mathbb{Z}}^{\ast }$, obtained by the MFP with a moving window of width $\mathbb{M}$, is obtained as

$$\begin{array}{cc}\hfill {\mathbb{Z}}_i^o=& \bigcup _{i=1,\dots ,\mathbb{N}}\left\{{\mathbb{Z}}_{{\tau }_i}:\left|{\mathbb{Z}}_{{\tau }_i}-{\widehat{\mu }}_i\right|⩾1.64·{\widehat{\sigma }}_i\right.\hfill \\ & \left.{\tau }_i\in ((i-1)·\mathbb{M}+1,i·\mathbb{M})\right\}\hfill \end{array}$$

(13)

Once the outliers are identified, they are replaced by normal values [70]. In this work, they are replaced by the median value price of the specific window period.

3.2. The Model

Once the filtered price series is obtained, it is modeled using the following model.

$${\mathbb{Y}}_i=D_i+{\mathbb{Z}}_{i}i=1,\cdots ,N$$

(14)

where ${\mathbb{Y}}_i$ is the filtered time series and ${\mathbb{Z}}_i$ is a stochastic term. The deterministic component captures the long trend, the yearly and weekly periodicity, and the bank holidays effect. Mathematically, it is defined as

$$D_i=l_i+y_i+w_i+b_i$$

where the terms $l_i,y_i,w_i,\mathrm{and}{b}_i$ represent the long-term trend, yearly periodicity, weekly periodicity, and bank holidays effect, respectively. In this work, the estimation procedure for the deterministic component described in [71] is used.

Once the deterministic component is estimated, the stochastic component ${\mathbb{Z}}_i$ is obtained as

$${\mathbb{Z}}_i={\mathbb{Y}}_i-D_i$$

(15)

which is modeled using the aforementioned FAR($\mathbb{P}$) and two alternate competing models. The alternate competing models used in this work are the univariate AR(P) model and a naïve benchmark model. The details of the competing models are as below.

3.2.1. Autoregressive (AR) Model

The univariate AR is one of the popular forecasting models used in time series analysis. It is similar to a regression model where the response variable is regressed over its lagged values. More specifically, in the AR modeling, a response variable is linearly dependent on its P lagged (past) values and an error term. Denoted by AR(P), mathematically, it can be written as

$$Z_i=\beta +\sum _{k=1}^P{\alpha }_k{Z}_{i-k}+{\epsilon }_i$$

(16)

where $Z_i$ is a univariate stationary time series, $\beta$ is a constant, ${\alpha }_k(k=1,\dots ,P)$ are the autoregressive parameters, and ${\epsilon }_i$ is a white noise process having zero mean and a constant variance. The choice of appropriate lag order selection is one of the most important steps in AR modeling. Different methods, including the Akaike information criterion (AIC) or Bayesian information criterion (BIC), or residual plots, e.g., autocorrelation function (ACF) and partial autocorrelation function (PACF), can be used to determine the lag order to be used in the model. In our work, the ACF and PACF are used, which indicate to use a restricted AR(7) model with ${\alpha }_k=0\mathrm{for}k=3,4,5,6$. The maximum likelihood estimation (MLE) method is used to estimate the parameters of the above model.

Once both the deterministic and stochastic components are modeled and forecasted, the final forecast is obtained as

$${\widehat{\mathbb{Y}}}_{i+1}={\widehat{D}}_{i+1}+{\widehat{\mathbb{Z}}}_{i+1}i=1,\cdots ,N.$$

(17)

The flowchart of the proposed general modeling framework is given in Figure 1.

3.2.2. The Naïve Benchmark

This section provides details about a naïve forecasting method that belongs to a similar day technique and has reported greater accuracy than other naïve methods [2]. This method works as follows.

• To forecast a given day, for example, Thursday, select the day before Thursday, which is Wednesday, and denote it by $x^{\ast }$.
• From the validation dataset, select all the Wednesdays (except $x^{\ast }$) and compare them with $x^{\ast }$ using the mean absolute error (MAE).
• Obtain a value of the MAE for each comparison that will result in a vector of the MAE values.
• Find and locate the smallest value of the MAE in the vector. Once the Wednesday having the lowest MAE is located, use its next day, i.e., Thursday, as the forecast for the concerned Thursday. This process is repeated for all the remaining days of the week.

4. Out-of-Sample Forecast

The dataset used in this empirical study includes electricity prices data called “Prezzo Unico Nazionale (PUN)” from the Italian Electricity Market (IPEX), collected from 1 January 2012 until 31 December 2017. Each day consists of 24 observations, where each observation corresponds to a load period. For modeling and forecasting purposes, we split the data into two periods. The period from 1 January 2012 to 31 December 2016 (1827 days) is used for model estimation. This period is used to optimize the parameters of the models. The out-of-sample period ranges from 1 January 2017 to 31 December 2017 (365 days). This period is used for forecasting the performance of the models. The one-day-ahead out-of-sample forecast is obtained through the window expending technique. In Figure 2, the spot electricity prices series is depicted for six years with a sample of functional (smoothed) curves for a week plotted on the right-hand side. The weekly periodicity is evident in the price time series as the prices profile for working days is relatively different from the non-working days.

The forecasting performance of the proposed and alternative models is compared using three standard descriptive forecast error measures. The point forecast accuracy is evaluated using three standard accuracy measures, namely, mean absolute percentage error (MAPE), MAE, and root mean square error (RMSE). Mathematically, the MAPE, MAE, and RMSE are given as

$$\begin{array}{ccc}\hfill MAPE& =& \frac{1}{N}\sum _{i=1}^N\frac{\left|{\mathbb{Z}}_{i,j}-{\widehat{\mathbb{Z}}}_{i,j}\right|}{{\mathbb{Z}}_{i,j}}\times 100\hfill \\ \hfill MAE& =& \frac{1}{N}\sum _{i=1}^N\left|{\mathbb{Z}}_{i,j}-{\widehat{\mathbb{Z}}}_{i,j}\right|\hfill \\ \hfill \mathrm{RMSE}& =& \sqrt{\frac{1}{N}\sum _{i=1}^N{[{\mathbb{Z}}_{i,j}-{\widehat{\mathbb{Z}}}_{i,j}]}^2}\hfill \end{array}$$

where N represents the number of observations in the out-of-sample forecasting period, ${\mathbb{Z}}_{i,j}$ denotes the original observed prices of the ith day and jth hour, and ${\widehat{\mathbb{Z}}}_{i,j}$ denotes the forecasted price of the aforementioned day and hour with $j=1,2,\dots ,24$.

In addition, directional forecast statistics can be very beneficial for traders in the electricity market in making investment decisions. These direction moments or turning points can be measured using directional statistic defined as [72]

$$\begin{array}{c}\hfill {\mathrm{D}}_{stat}=\frac{1}{N}\sum _{i=1}^N{\alpha }_i\ast 100\end{array}$$

where

$${\alpha }_i=\left\{\begin{array}{cc}1,\hfill & \mathrm{if}({\mathbb{Z}}_{i+1,j}-{\mathbb{Z}}_{i,j})({\widehat{\mathbb{Z}}}_{i+1,j}-{\mathbb{Z}}_{i,j})\ge 1\hfill \\ 0,\hfill & \mathrm{otherwise}\hfill \end{array}\right.$$

The electricity prices forecast through the FAR($\mathbb{P}$) model have the following steps. In the first step, the moving window filter method was used for the identification and accommodation of outliers. In the second step, a logarithm (log) transformation was performed to stabilize the variance of the series. In the third step, model (14) is applied to the data and the series ${\mathbb{Z}}_i$ is obtained using Equation (15). In the fourth step, the Fourier basis functions are used to transform the discrete data into functional data to obtain 2192 daily functional trajectories, say, ${\mathbb{Z}}_1\left(t\right),\dots ,{\mathbb{Z}}_{2192}\left(t\right),t\in J$. Once the functional data are obtained, the FAR($\mathbb{P}$) model described in Section 2.1 is applied, and one-day-ahead forecasts are obtained for the whole year. In the case of the competing models, the univariate AR model and the naïve benchmark are applied directly to ${\mathbb{Z}}_i$, and the one-day-ahead forecasts are obtained for the whole out-of-sample period.

Figure 3 highlights the population mean function $\mu \left(t\right)$ and the functions obtained by adding and subtracting a suitable multiple of the eigenfunctions to the mean. Such plots are helpful to understand the variability in the direction of certain eigenfunctions. The first eigenfunction is positive, indicating that subjects with positive scores on this component will contribute to obtaining a consistently larger proportion ($77.1\%$) of the total variation of the data. The second eigenfunction displays an oscillatory behavior, suggesting that subjects with positive scores will have lower electricity prices from midnight till early morning and then slightly more between hours 7 a.m. and 10 a.m., and explain $10.5\%$ of the total variation of the data. Similarly, the third and fourth eigenfunctions explain $4.8\%$ and $2.7\%$ of the total variation of the data, respectively. The first four eigenfunctions collectively explain more than $95\%$ of the total variability in the electricity prices data.

Concerning the forecasting results for the proposed and alternative models,

Table 1. IPEX electricity prices: out-of-sample forecasting errors MAE, MAPE, and RMSE for FAR($\mathbb{P}$), AR(7), naïve models, and the directional statistics D${}_{stat}$ with number of forecasting directions (same, up, down).

Model	MAE	MAPE	RMSE	Dstat (%)	SAME	UP	DOWN
FAR($\mathbb{P}$)	5.16485	8.99009	8.65032	88.34342	1525	3743	3492
AR(7)	5.65833	10.09469	9.20305	82.95525	1272	4084	3404
Naive	6.86278	12.63467	10.09929	53.63626	385	4137	4238

compares the overall forecasting ability of the FAR($\mathbb{P}$), AR(7), and naïve models through out-of-sample forecasting errors computed by MAE, MAPE, and RMSE. The table also provides the directional forecasting performance for these models. From the results, it is evident that our proposed functional model performs significantly better than the other competing models. The proposed FAR($\mathbb{P}$) models produce MAE, MAPE, and RMSE of 5.16, 8.99, and 8.65, respectively. Although the univariate AR model produces better results than the naïve model, it produces considerably higher forecasting errors compared to the proposed functional model. Looking at the directional forecasting results, note that the value of D ${}_{stat}$ for FAR($\mathbb{P}$) is 88.34%, whereas values of 82.96% and 53.64% are obtained in the case of AR and naïve models, respectively. Hence, our proposed functional model performs relatively well compared to the competing models. From the number of forecast direction moments, it can be seen that the FAR($\mathbb{P}$) forecast 1525 out of the total 8760 load periods accurately (the“SAME” in Table 1 refers to the absolute difference of the forecasted value minus the actual value to be less than EUR 1), whereas this value for the AR and naïve models is 1272 and 385, respectively. The number of over-forecasted values for FAR($\mathbb{P}$) and AR(7) are 3743 and 4084, respectively. Again, the poor performance of the naïve model is evident from the results of the directional forecast.

Table 2. IPEX electricity prices: daily forecast errors for FAR($\mathbb{P}$), AR(7), and naïve models.

Model	Error	Days of a Week
		Monday	Tuesday	Wednesday	Thursday	Friday	Saturday	Sunday
FAR($\mathbb{P}$)	MAE	5.63112	5.87795	6.53440	5.00897	4.88412	4.06196	4.17448
AR(7)		6.16619	6.34358	6.95023	5.86286	5.45145	4.54586	4.31402
Naive		7.58172	7.38245	6.39990	6.93643	6.40897	6.86697	6.47056
FAR($\mathbb{P}$)	MAPE	9.87044	9.23202	9.30545	7.72602	7.85736	8.58661	10.32703
AR(7)		11.14407	10.06261	10.35827	9.02434	8.86721	9.91458	11.26919
Naive		13.82605	12.47558	11.01666	10.93947	11.04299	13.52881	18.26894
FAR($\mathbb{P}$)	RMSE	8.72362	9.56184	11.93833	9.59608	7.97935	5.43820	5.36185
AR(7)		9.34377	10.08267	11.97317	10.87680	8.79962	5.86412	5.60787
Naive		11.60842	10.63521	10.37988	10.76548	8.88560	9.54221	8.54465

reports the daily forecast accuracy for the electricity prices using different models. From the table, one can see that the FAR($\mathbb{P}$) model produces lower forecasting errors compared to the univariate AR(7) and naïve models. Although the forecast errors vary from day to day, they are lower on Thursday and Friday when considering MAPE. The poor performance of the naïve model is evident from this table. The hourly forecast errors for different models are listed in

Table 3. IPEX electricity prices: hourly forecast errors for FAR($\mathbb{P}$), AR(7), and naïve models.

Model	Hour	MAE	MAPE	RMSE	Hour	MAE	MAPE	RMSE
FAR($\mathbb{P}$)	1	3.17448	10.32703	5.36185	13	5.03299	9.77343	7.50008
AR(7)		4.09499	8.20630	5.48352		5.50527	10.81530	8.32688
Naive		2.99810	5.71016	4.35740		5.89142	10.72254	8.54265
FAR($\mathbb{P}$)	2	5.63112	9.87044	8.72362	14	5.13495	10.97364	7.45864
AR(7)		3.90162	8.62680	5.04521		5.69565	12.52584	8.52711
Naive		3.99566	8.14583	5.04110		5.15640	10.87232	7.39582
FAR($\mathbb{P}$)	3	5.87795	9.23202	6.46184	15	4.91392	10.19534	9.29372
AR(7)		3.82670	9.05766	4.96036		6.52555	13.92227	10.37498
Naive		3.66855	8.12195	4.42199		5.40823	11.21776	7.62226
FAR($\mathbb{P}$)	4	3.39714	8.41575	4.46110	16	6.20807	11.99534	9.73345
AR(7)		3.97490	10.11668	5.28032		6.91072	13.74287	10.63356
Naive		3.52617	8.48418	4.64024		6.07377	11.75996	9.43945
FAR($\mathbb{P}$)	5	3.39652	8.35320	4.46989	17	6.63385	10.73122	10.96537
AR(7)		3.89939	9.77731	5.18900		7.02079	11.61225	11.34363
Naive		6.33318	12.83398	9.91261		6.47302	10.62784	10.48912
FAR($\mathbb{P}$)	6	4.88413	7.85736	7.97935	18	6.95788	9.73444	12.11111
AR(7)		3.73582	8.64474	5.02575		7.32959	10.47554	12.35075
Naive		7.59001	8.59831	4.74848		7.04453	13.45730	10.9536
FAR($\mathbb{P}$)	7	4.06196	8.58661	5.43820	19	7.83811	9.91660	14.00040
AR(7)		4.27194	8.50246	5.79461		7.99731	10.31619	13.91400
Naive		4.98340	11.28051	6.28747		7.35881	11.86122	11.68502
FAR($\mathbb{P}$)	8	4.94668	8.16684	7.92432	20	7.35007	9.7326	11.80230
AR(7)		5.60254	9.18891	9.23695		7.40072	9.98525	11.66830
Naive		6.25091	11.78682	8.76700		8.46851	12.49491	13.42239
FAR($\mathbb{P}$)	9	7.28609	10.34856	12.47580	21	6.28777	9.03979	9.63504
AR(7)		8.17881	11.92595	13.42361		6.26023	9.14635	9.46692
Naive		7.14544	11.18423	11.62178		8.05637	12.35241	11.55807
FAR($\mathbb{P}$)	10	6.73647	9.88567	11.63280	22	5.11291	7.86378	8.41742
AR(7)		7.67803	11.57003	12.60697		5.18257	8.05526	8.48756
Naive		6.95444	10.27270	11.25744		7.08351	10.66474	10.11036
FAR($\mathbb{P}$)	11	6.09630	9.80792	9.88950	23	3.86946	6.64648	6.09351
AR(7)		6.80914	11.15632	10.81531		3.93445	9.78030	6.12057
Naive		6.40122	10.11695	9.80316		6.14502	9.13875	9.47691
FAR($\mathbb{P}$)	12	5.79828	10.13717	8.84841	24	3.55292	6.70792	5.27639
AR(7)		6.45438	11.43950	9.81589		3.55885	6.69256	5.32193
Naive		5.94363	10.01007	8.73133		4.49332	7.98897	6.5533

, which shows that the forecast errors vary throughout the day. Although the FAR($\mathbb{P}$) model produces better results on most hours, the AR(7) has better results on two hours when considering the MAPE. It is worth mentioning that the proposed FAR($\mathbb{P}$) model performs significantly well during peak hours compared to the competing models. Again, the poor performance of the naïve model is evident from the results.

Finally, the results obtained by our proposed functional model in this study are compared with the results listed in the literature. Here, it is worth mentioning that such a comparison is only to evaluate the performance of our model, as different authors considered different forecasting horizons, different periods, and different error summary measures. Using the Italian electricity market and considering a one-day-ahead forecast, Ref. [2] obtained an MAPE value of 9.74 using the NPAR model, which is significantly higher than our proposed model MAPE value of 8.99. The research work of [73] used the Italian electricity market data, and their proposed model produced an MAE of 8.58, whereas our proposal reported an MAE value of 5.16, 60% lower. For a one-day-ahead forecast, Ref. [70] reported an MAPE value of 9.05, which is slightly higher than our obtained MAPE value. Using an ARX-EGARCH model for the Italian electricity prices time series, Ref. [74] obtained an RMSE of 11.58, whereas our proposed model produced an RMSE value of 8.65. The work of [75] reported RMSE values of 16.72 and 15.79 using ARMA and GARCH models, respectively, significantly higher than our value of 8.65.

5. Conclusions and Future Direction

In today’s competitive electricity market, modeling and forecasting electricity prices are critical for market participants to optimize their strategies. However, electricity prices exhibit specific features, including long-trend, periodicities, spikes or jumps, bank holidays, etc. In the presence of these features, the forecasting problem is a great challenge for researchers. This paper proposes a functional model for modeling and forecasting electricity prices. To this end, the price time series is first treated for the extreme values. The filtered series is then divided into deterministic and stochastic parts. The deterministic part modeled the effects of long-trend, annual, and weekly periodicities, and bank holidays. For the stochastic component, a functional AR model (FAR) is proposed that is capable of automatic selection of lags and dimensions. To evaluate the performance of our proposed model, two alternate models, namely, the univariate AR and a naïve benchmark, are also used in this study. For empirical comparison, data from the Italian electricity market are used and the out-of-sample one-day-ahead forecast errors measured through MAPE, MAE, and RMSE are calculated for a complete year.

The empirical results suggest that the proposed FAR($\mathbb{P}$) model is significantly better than the competing model, as it produced considerably lower forecasting errors. Furthermore, the component estimation procedure is highly effective in forecasting electricity prices. Moreover, the directional forecast results suggest that this approach can significantly increase the number of accurate forecasts. Accurate forecasting can be very helpful for the traders (buyers and suppliers) to optimize their bidding strategies to maximize their gains and to use the resources required for electricity generation more effectively. Consequently, this will also benefit the end-user in terms of reliable and economical electricity facilities.

As the current study does not consider any exogenous variable effect in the model, this effect can be investigated in the future. Furthermore, as the current study only considers linear models, nonlinear models can also be compared with the proposed functional model.