# Functional Data Approach for Short-Term Electricity Demand Forecasting

Authors: Ismail Shah; Faheem Jan; Sajid Ali
Journal: Mathematical Problems in Engineering
 (2022)

---

## Abstract

In today’s liberalized electricity markets, modeling and forecasting electricity demand data are highly important for the effective management of the power system. However, electricity demand forecasting is a challenging task due to the specific features it exhibits. These features include the presence of extreme values, spikes or jumps, multiple periodicities, long trend, and bank holiday effect. In addition, the forecasts are required for a complete day as electricity demand is decided a day before the physical delivery. Therefore, this study aimed to investigate the forecasting performance of models based on functional data analysis, a relatively less explored area in energy research. To this end, the demand time series is first treated for the extreme values. The filtered series is then divided into deterministic and stochastic components. The generalized additive modeling technique is used to model the deterministic component, whereas functional autoregressive (FAR), FAR with exogenous variable (FARX), and classical univariate AR models are used to model and forecast the stochastic component. Data from the Nord Pool electricity market are used, and the one-day-ahead out-of-sample forecast obtained for a whole year is evaluated using different forecasting accuracy measures. The results indicate that the functional modeling approach produces superior forecasting results, while FARX outperforms FAR and classical AR models. More specifically, for the NP electricity demand, FARX produces a MAPE value of 2.74, whereas 6.27 and 9.73 values of MAPE are obtained for FAR and AR models, respectively.


---

## Body

## 1. Introduction

Electricity is one of the essential needs in today’s world. It is a basic necessity for every person to complete day-to-day activities. However, it poses a great challenge to every country of ensuring the accessibility of a reliable and cost-efficient availability of electricity. Electricity demand forecasting is an important issue for different percipients in today’s liberalized electricity market. Efficient and reliable demand forecasting plays a crucial role in electricity generation capacity, transmission planning, and pricing. Generally, electricity forecasting is required from three different forecasting perspectives. Long-term electricity demand forecasting is needed for capacity planning that is generally based on economic or demographic variables. Medium-term forecasts are necessary for scheduling, maintenance, and acquiring resources to ensure an efficient electricity demand. On the other hand, a short-term electricity forecast, which is required for a day-ahead market, is helpful in day-to-day market activity management, risk assessment, optimizing bidding strategies, and increasing the profit margins. In the literature, short-term electricity demand forecasting received greater attention as maximum electricity is traded in this market [1].The liberalization brings many benefits to the end user; however, the electricity demand forecasting became challenging due to the specific features that the demand series exhibit. In general, the electricity demand time series contains a long trend, multiple periodicities, the bank holiday effect, and the presence of extreme values. Over time, electricity demand increases, introducing a long trend in the demand time series. The yearly periodicity refers to the cyclical demand behavior from one year to another that is evident in the demand time series. The yearly seasonality is evident from the graph as the electricity demand pattern over a year is similar to the adjacent year. The demand is higher in winter compared with other seasons. The load profile on working days shows different behaviors than the nonworking days resulting in weekly periodicity. Moreover, the demand series shows high volatility and unpredicted spikes/jumps caused by sudden shutdowns of power plants, unpredictable weather fluctuations, etc. In the presence of these specific features, forecasting electricity demand is a challenging task, and continuous efforts are made to predict it accurately [2–4].Many statistical, econometric, machine learning, and hybrid models are extensively used in the literature for electricity demand forecasting. These methods generally vary in methodology, complexity, and performance [5–9]. The most common methods used for forecasting electricity demand are based on linear time-series models such as autoregressive (AR) models, AR moving average (ARMA), AR integrated MA (ARIMA), seasonal ARIMA (SARIMA), and their different extensions [10–12]. Nonlinear time-series models, for example, nonparametric AR (NPAR), AR conditional heteroscedasticity (ARCH), generalized ARCH (GARCH), and various extensions of GARCH, are widely used for forecasting electricity demand and prices [13, 14]. In addition, exponential smoothing models, i.e., simple, double, and triple Holt-Winters models, are easy to implement and can accommodate multiple periodicities. These types of models are widely used in different fields of study, including energy markets. For example, using data from Pakistan, Rehman et al. [15] used three forecasting models, namely ARIMA, Holt-Winters, and long-range energy alternative planning (LEAP) methods to model and forecast the energy consumption of five fuels, i.e., electricity, oil, coal, natural gas, and liquefied petroleum gas (LPG) in six different fields including domestic, industrial, commercial, transportation, agriculture, and other governmental sectors. Data from 1992 to 2014 are used in the study, and the energy consumption forecasts are made for the coming 21 years. The study results suggested that the ARIMA model was the most appropriate model to predict the energy demand compared with the rest. Bin [16] used the ARIMA model for the short-term prediction of monthly load data in Dayton. The results indicated that ARIMA produced better results than other models used in the study. The exponential smoothing method for medium-term and long-term demand forecasting is used by Lv et al. [17], which is suitable for large areas with relatively little uncertainty, and some seasonal factors can also be included in the model. The results indicated that their proposed model produces better results compared with other predictive models. On the other hand, parametric and nonparametric regression-type models, including multiple regression, local polynomial regression, kernel regression, smoothing splines, and quantile regression models, are also used for electricity demand and price forecasting [18–20].Recently, the artificial intelligence models are also used to predict day-ahead electricity demand [21, 22]. For instance, using monthly data from Canada, Bouktif et al. [23] forecasted the short- and medium-term electric load using the long short-term memory (LSTM)-based recurrent neural networks (LSTM-RNNs). The results of the proposed models are compared with other machine learning models, including linear regression (LR), ridge, regression K-nearest neighbors (KNN), random forest (RF), gradient boosting (GB), ANN, and extra trees regressors. Besides the electricity demand time series, other predictors such as temperature, humidity, wind speed, and schedule-related variables (month number, weekends, and weekdays) are also used in their study. The forecasting performance of different models is evaluated using different descriptive statistic measures. The results suggested the superior performance of their proposal compared with the rest. To forecast electricity demand or prices, various researchers combined the characteristics of two or more models resulting in a hybrid model [24–26]. In general, every model has its own functional and structural form, and the forecasting performance varies from market to market.The main aim of this research was to propose forecasting models for electricity demand based on the functional data analysis, a relatively less explored field of statistics in the energy sector. To this end, the main contribution of this study is the use of functional models within the approach based on component estimation. The functional model is capable of selecting the lags and dimensions automatically. The forecasting results are obtained for a complete year to assess the performance of the proposed approach. Furthermore, different descriptive measures are used to compare the forecasting performance of the proposed approach with those reported in the literature.The rest of the article is structured as follows. Introduction to functional data analysis (FDA) and literature review related to FDA is provided in Section2. Section 3 presents the general modeling framework, including a detailed description of extreme value treatment and different component estimation procedures. In Section 4, an analysis of the Nord Pool electricity market is provided, whereas empirical results for different models used in this study are discussed in Section 5. Finally, the conclusion is provided in Section 6.


## 2. Functional Data Analysis

Modern technological developments have simplified and decreased the cost of the data collection and storage process, which allows us to get an increasingly common sketch of several real-life phenomena. This novel idea has made data scientists to confront big data sets with complicated structures. Analyzing big data gives both benefits and challenges in different research areas, including engineering, medicine, public health, environmental science, and finance. Big data provide substantial facts for statistical inference, while on the other hand, complicated frameworks, large sample sizes, and high dimensions make the analysis of such data a challenging task. Functional data analysis (FDA) has established itself as an important and dynamic area of statistics to deal with such complex data and has received greater attention in recent years. It also provides new tools and has simulated new methodological and theoretical developments. The FDA generally arises in situations to exploit the information recorded over a continuum, such as time or space. Therefore, this area has become very broad, with multiple specialized directions of investigation. The monographs by Ramsay and Silverman [27, 28]; Locantore et al. [29]; Ferraty and Vieu [30]; and Ferraty and Romain [31] provide a comprehensive sketch of the methodological and theoretical improvements. Functional principal component analysis (FPCA) plays a crucial job in the development of the FDA that is widely used in functional forecasting [32].In the FDA, the observed time series is considered discrete observations of a continuous function. With a real and parsimonious functional representation, the large-dimensional data are transformed into a sequence of functional trajectories that can be evaluated with enhanced efficiency and accuracy. Functional time series (FTS) usually consists of stochastic functions observed at regular time intervals. FTS also arises when the observations in a period are considered definite realizations of an underlying continuous function. While there are a variety of methods dealing with discrete time series, there are fewer contributions paying attention to the FTS [33, 34]. The FTS analysis is so far a fast-developing area of research. Reference [35] derived a functional Yule–Walker estimator for the serially dependent functional data that are likely the most popular pioneering study that plays a prominent role in FTS. Reference [36] studied the theoretical structure to explore the serial dependence in FTS modeling. Reference [37] proposed a dimension reduction technique to model the functional ARMA (FARMA) model with an application to traffic data. Among the existing literature, the FTS models have been employed to forecast intraday trajectories; see, for example, Shang [38, 39] and Horváth et al. [40].Due to the popularity of the FDA for high-dimensional data, the FTS methods are also used in energy markets [41–44]. Within the framework of the FTS, several approaches have been proposed. For instance, Goia et al. [45] studied the functional linear model for the peak load demand forecast. Antoniadis et al. [46] developed a nonparametric functional method based on the functional kernel regression estimator for forecasting half-hourly electricity consumption in France. The research of Antoch et al. [47] used a parametric functional model to forecast electricity consumption curves, whereas Cho et al. [48] studied a hybrid method that was applied to the French demand data. Chen and Li [49] suggested an adaptive functional AR (AFAR) prediction model to predict electricity price trajectories. Li et al. [50] suggested FAR fractionally integrated moving average (FARFIMA). Using monthly sea surface temperature data, Li et al. [51] study stationary functional time series with long-range dependence and estimate the memory parameter involved. Forecasting electricity price time series for the Italian electricity market, Jan et al. [52] used the functional autoregressive model and compared it with classical and naive models. The results suggested that the functional model produces significantly better results than the competing models.


## 3. Methods

The general framework for modeling and forecasting electricity demand is provided in this section. As described in Section1, electricity demand exhibits specific features such as extreme values, long trend, multiple periodicities, and the bank holiday effect. Incorporating these features into the model greatly improves forecasting accuracy. Thus, the series is preprocessed for extreme values before modeling the demand series.


### 3.1. Moving Filter on Demand

The presence of outliers may significantly affect the estimation and the forecasting. The outliers are mostly present in the time-series data, making the estimation and forecasting of the time-series model more challenging. In time-series analysis, outliers’ detection and replacement with normal values are the essential step of the data cleaning process that, generally, improves the estimation of parameters and the accuracy of forecasting models. In electricity demand or price data, outliers (also known as extreme values) are generally present due to many factors, including unexpected increased demand, power plant shutdown, and extreme weather conditions. The box plot of hourly electricity demand is plotted in Figure1. One can see that there are some extreme values (outliers) present in the demand series, especially at midnight hours. In general, these extreme values are identified and replaced by many methods [53]. In this research work, they are treated by the moving filter on-demand method suggested by Borovkova and Permana [54].Figure 1
Nord Pool electricity demand: hourly box plots for demand data range from January 1, 2015, to December 31, 2019.The moving filter on-demand (MFD) is a generalization of the standard deviation filter on-demand (SFD) technique. The SFD works on the whole time series, whereas the MFD operates on a rolling window of constant width. At the step, the original time series is divided intow=q/r parts, where r represents the width of the window and q is the total number of observations in the series. Then, MFD is applied to the first window of the original time series. The demand values whose absolute deviation from the estimated mean μ^ is greater than a multiple of estimated standard deviation (SD) are identified as extreme values. The subset of spikes Xto detected with the use of the MFD filter is obtained as follows:(1)Yto=Yt:Yt−μ^≥1.64⋅SDY.The window then movesr observations forward, and the filtering method is repeated in the next window to determine the extreme values. The procedure is repeated until this method treats all w parts. Once the extreme values are identified, they are replaced by the median value of the corresponding window to obtain extreme value free demand series.


### 3.2. General Model

After possible outliers’ identification and replacement with the normal values to obtain the demand series free of outliers, the next step is to model and forecast electricity demand. To this end, the component estimation procedure is used to capture different features of the demand series [55]. Suppose the logarithmic electricity demand series is log Yt,j, where t∈N and j=1,…,24. The demand series is decomposed into a deterministic component denoted by Dt,j and a stochastic component denoted by Xt,j as follows:(2)logYt,j=Dt,j+Xt,j.Here, the deterministic component in (2) is comprised of the complex structure of the seasonal component, such as annual seasonality st,j, weekly seasonality wt,j, bank holiday effect bt,j, and long-trend component ct,j. Mathematically, it can be written as follows:(3)Dt,j=ct,j+st,j+wt,j+bt,j.The deterministic componentDt,j in (3) is estimated through a generalized additive model (GAM). In the case of stochastic component Xt,j, three different models, namely the FAR (p), FARX (p), and AR (p), are used. Once both the components are modeled and forecasted, the one-day-ahead out-of-sample forecast is obtained as follows:(4)Y^t+1,j=expD^t+1,j+X^t+1,j.


### 3.3. Modeling the Deterministic Component

This section will explain the modeling and forecasting of the deterministic component. The generalized additive modeling technique is used to model the long-term trendct,j, annual st,j, weekly wt,j periodicities, and bank holiday effect bt,j. More specifically, the long-run (trend) component ct,j, which is a function of time t, is estimated through the ordinary least-squares (OLS) technique. In case of the yearly seasonality, the series 1,2,…,365,1,2,…,365,…,1,2,…,365 is used and is estimated using the OLS. On the other hand, for weekly periodicities and for the bank holidays, dummy variables are used; i.e., wt,j=∑i=17φict,j with ct,j=1 if t refers to the ith day of the week and zero otherwise and bt,j=∑i=12ωct,j with ct,j=1 if t refers to a bank holiday and zero otherwise. The coefficients φi and ωi are again estimated with the OLS approach. Once the deterministic components are estimated, the one-day-ahead forecast for c^t+1,j=c^t,j and s^t+1,j=s^t,j, as both of these components represent long-term dynamics for our forecast horizon, whereas the forecast of wt and bt is straightforward as both of these are deterministic functions of time or the calendar conditions. Hence, one-day-ahead forecast for the deterministic component can be obtained as follows:(5)D^t+1,j=c^t,j+s^t,j+w^t+1,j+b^t+1,j.To obtain the stochastic component, we subtract the deterministic component fromlogYt,j as follows:(6)X^t,j=logYt,j−D^t,j=logYt,j−c^t,j+s^t,j+w^t,j+b^t,j.The modeling procedure of the stochastic component is discussed in the next section in detail.


### 3.4. Modeling the Stochastic Component

This section describes the estimation of the residual componentXt,j, which is modeled using the functional autoregressive, FAR (p), functional autoregressive with exogenous variables, FARX (p), and the classical univariate autoregressive, AR (p), models. The classical AR (p) is a well-known model; however, before explaining the functional models, some preliminaries are required, which are stated in the next section.


#### 3.4.1. Basics of Functional Data Analysis

This section provides some preliminaries that are essential for building the functional models. The stochastic componentXt,j in (2) is transformed into functional data using some known basis functions. Using suitable basis functions, the functional object corresponding to the daily demand profile for the day t is represented as follows:(7)Xtτ=∑k=1Kakγt,kττ∈J,where ak is the constant parameters and γt,kτ is the Fourier basis functions. An example of functional time-series trajectories for 1826 days is given in Figure 2 where each functional trajectory represents a daily demand profile.Figure 2
Nord Pool electricity demand: functional daily trajectoriesXtτt=1,…,1826, where Xtτ is the log transformed functional observations of hourly electricity demand, based on 15 Fourier basis functions.Assume thatXtτ given in (7) is an arbitrary functional time series (FTS) defined on a common probability space Φ,∑,P, where Φ,∑, and P denote the sample space, σ−algebra on Φ, and the probability measure on ∑, respectively. In general, it is assumed that the functions Xτ are elements of the semi-metric and Hilbert space. Furthermore, it is assumed that functions are the elements of the square-integrable functions X∈ℒ20,1 residing in the Hilbert space ℋ satisfying Xt2=∫X2τdτ<∞ with an inner product x,y=∫xτyτdτ,∀x,y∈ℒ20,1. The notation X∈ℒℋυΦ,∑,P is used to indicate for some υ>0, EXυ<∞. Note that υ=1 results in the population mean curve μτ=EXτ, and when υ=2, then nonnegative definite covariance operator is obtained as follows:(8)Cτ,δ=CovXτ,Xδ,Cτ,δ=EXτ−μτXδ−μδ.The covariance operatorCτ,δ in (8) allows the covariance operator of X, denoted by G, which can be expressed as follows:(9)Gαδ=∫10Cτ,δατdτ.Using Mercer’s lemma, there is an orthonormal sequenceαk of continuous functions in ℒ20,1 and a non-increasing sequence βk of positive numbers, such that(10)Cτ,δ=∑k=1∞βkαkταkδ,where the αkτ denotes the kth functional principal component (FPC) and βk denotes the kth eigenvalue in the decreasing order [56]. By the separability of Hilbert spaces, the Karhunen–Loève (KL) expansion of a random process Xτ can be expressed as follows:(11)Xtτ=μt+∑k=1∞λk,tαkτ,where λk,t denotes the kth principal component score (PCS) defined as λk,t=∫Yταkτdτt. The PCSs constitute an uncorrelated sequence of random variables with zero mean and variance βk.The expansion 3.7 facilitates dimension reduction as the firstd terms often provide a good approximation to the infinite sums, and thus, the information contained in Xtτ can be adequately summarized by the d-dimensional vector β1,…,βd. The approximated processes can be defined as follows:(12)Xtτ=∑k=1dλk,tαkτ+ετ,where d denotes the number of retained principal components, and ετ denotes the error function with a mean of zero and a finite variance, containing the FPCs excluded from the retained d components. For more details about the FPC analysis and its practical demonstration, the interested reader may consult Ramsay and Silverman [57]; Shang [58]; and Aue et al. [59].In practice, the mean curveμτ, FPC’s ατ, and PCS’s λ can only be estimated through realizations of a random process. Suppose X1τ,…,XNτ are fully observed FTS, and then, the standard estimators for the mean function μt and covariance operator Cτ,δ are given by the following sample averages:(13)μ^τ=1N∑t=1NXtτ,(14)C^τ,δ=∑k=1∞β^kα^kτα^kδ,where β^1>β^2>⋯≥0 are the sample eigenvalues of C^τ,δ, and α^1τ,α^2τ,… are the corresponding orthogonal sample eigenfunctions. Hörmann and Kokoszka [36] studied the aforementioned estimators and concluded that they are consistent estimators for weekly dependent process. Then, using the KL expansion, the realizations of the random process X can be written as follows:(15)Xtτ=μ^τ+∑k=1dλ^k,tα^kτ+eτ,t=1,2…,N,where λ^k,t is the kth estimated PCS for the tth observation. Figure 3 provides the cumulative proportion of variation explained by the first three PCs. Note that the first three PCs explained more than 99% of the total variation. Figure 4 represents the first three functional principal components with the amount of variation explained by each FPC. Note that the 1st eigenfunction explains approximately 95.2% of the total variation in the data.Figure 3
Cumulative proportion of variation explained by the first three PCs.Figure 4
First three functional principal components with amount of variation explained by each FPC.The next section provides details about FAR (p) and FARX (p) models used to model and forecast Xtτ.


#### 3.4.2. Functional Autoregressive Model of Order p

The forecasting of the FTS is a complicated task for which several techniques are discussed in the literature. Functional autoregressive (FAR) models are one of the most popular functional models used for FTS forecasting. Consider that the functional trajectoriesXtτ are available for t=1…,N. In the FAR modeling, Xtτ is linearly dependent on its p lagged (past) values and an error term. Mathematically, the FAR (p) model can be written as follows:(16)Xtτ=μτ+∑l=1pΨiXt−lτ+ξtτ,where ξtτ∈N is independently and identically distributed sequence in ℒℋ2 with Eξt=0 and Ψll=1p is bounded linear operator mapping ℒℋ2⟶ℒℋ2 such that the above recursive equation has a unique causal solution. To model and forecast Xtτ in (16) using the FAR(p) model, (13)–(15) are used for modeling and forecasting in the following three steps.(i)
The dimensiond is fixed and the estimated FPC scores are obtained as λ^k,t=∫X^tτα^kτdτ for each observation X^tτ,t=1,…,N,k=1,…,d, and the estimated k-variate FPC score vectors λ^t=λ^1,t,…,λ^d,t′,t=1,…,N.(ii)
For a fixed orderp, the vector autoregressive model (VAR(p)) Yt=∑l=1pΨlYt−l+εt for eigenscore vectors is constructed to produce forecasting λ^N+1=λ^1,N+1,…,λ^d,N+1′. Given the vectors λ^1,…,λ^N, the Durbin–Levinson and innovations algorithm can be readily applied here.(iii)
In the last step, the multivariate time series is re-transformed to functional version using the KL theoremX^N+1τ=μ^τ+λ^1,N+1α^1τ+⋯+λ^d,N+1α^dτ. Based on the estimated FPC scores and sample eigenfunctions, the resulting X^N+1τ is then used as one-step-ahead forecast of XN+1τ.


#### 3.4.3. Forecasting with Exogenous Variables FARX (p)

In the FTS modeling, not only the lagged values of a response variable but also other explanatory exogenous variables are added to the model to increase the forecasting accuracy. These exogenous variables may be scalars, vector-valued, or functional. The objective is then to obtain an empiricalX^N+1 given observations of the trajectories X1,…,XN and the number of exogenous variables Xf1,…,Xfv. The predictor variables need not be defined on the same space. For example, Xf1 could be vector, Xf2 a scalar, Xf3t a function, and Xf4t could contain lagged values of Xf3. The FAR model having functional type of exogenous variables is generally denoted by FARX(p). The functional exogenous variables Xf1t,…,Xfρt are assumed to bestationary, with time-invariant mean functions μf1t,…,μfρt. Suppose Xtτ, t∈N, and the exogenous variables Xfρt, f∈ℕ, are available. Then, the FARX(p) model can be defined as [60]follows:(17)Xtτ−μt=∑l=1pΨlτXt−lτ−μτ+∑m=1ρϕmτXfmτ−μfmt+εtτ,where Ψlτ and ϕmτ are the functional operators of the model and εtτ is the functional error term. It is worth mentioning that the dynamics of the functional response can be represented not only with the lagged functional exogenous variables, but there may also be some causal dependence with scalar exogenous variables. Suppose that ρ scalar exogenous variables denoted by Xf1,…,Xfρ having mean μf1,…,μfρ are available. Then, the FARX(p) model can be written as follows:(18)Ytτ−μt=∑l=1pΨlτXt−lτ−μτ+∑m=1ρϕmXfm−μfm+εtτ.Both FARX(p) models in (17) or (18) corresponding to having functional or scalar exogenous variables, respectively, are the generalization of the model discussed in (16) with different numbers of functional or scalar exogenous variables. The estimation method discussed in Section 3.4.2 is modified and summarized as follows:(1)
(a)
The dimensiond is fixed for t=1,…,N and the data X1,…,Xt are used to compute vectors λ^t=λ^1,t,…,λ^d,t′ that contain the first d estimated FPC scores.(b)
In case of functional exogenous variables, the value ofρ is fixed, and the data X1,…,Xtt=1,…,N are used to compute the vector φ^t=φ^1,t,…,φ^v,t′, which contain the first ρ leading FPC scores. The process for each functional exogenous variable is repeated.(c)
In the next step, all the exogenous variable vectors are combined into a single vector,θ^N=θ^1,N,…,θ^ρ,N′.(2)
λ^1,…,λ^N and θ^N are used to obtain one-step-ahead forecast as follows:(19)λ^t+1=λ^1,t+1,…,λ^d,t+1′.(3)
The KL expansion is used to obtain(20)X^t+1τ=μ^τ+λ^1,t+1α^1τ+⋯+λ^d,t+1α^dτ,as one-day-ahead forecast for Xt+1 using the KL theorem.


#### 3.4.4. Dimension and Order Selection of the Functional Models

The selection of dimensiond and order p is an important step for the FTS forecasting using FAR (p) and FARX (p). In the literature, various methods are used to select the order of FAR. For example, Kokoszka and Reimherr [61] proposed a multistage testing hypothesis method for the selection of the order of FAR. Salmerón and Ruiz-Medina [62] studied the multispectral decomposition to select the order of the FAR (p) model. One common approach is choosing a minimum number of PCs so that the proportion of variance explained by the FPCs exceeds a certain threshold value. In this study, we used a completely automatic method for choosing order p and dimension d for both FAR (p) and FARX (p), which minimizes the mean square error of forecasting studied [63], which is asymptotically equivalent to the functional final prediction error (fFPE) given as follows:(21)fFPEp,d=t+p∗dt−p∗dtrΔ^Y+∑k>dβ^k,where βk is the kth eigenvalue of Cτ,δ and Δ^Y in (21) is the covariance matrix of the random vector Y1,…,Yd and is an unbiased estimator of ΔY. The best d and p are the minimizer of the fFPE function.


#### 3.4.5. Univariate Autoregressive Model

Autoregressive models are the most commonly employed linear models in univariate time-series forecasting problems. In the autoregressive model, the dependent variable is regressed on itsp lagged values. The number of lag values required to forecast the future value is the order of an AR model, which is denoted by the value p. Mathematically, the autoregressive model of order p, AR (p), can be written as follows:(22)Xt=ϕ+Ψ1Xt−1+⋯+ΨpXt−p+εt,where ϕ is the constant (intercept), Ψll=1,…,p is the parameters of AR (p) model, and εt is assumed to follow a white noise process with mean zero and a constant variance. In the literature for demand modeling, an AR (7) is often used. Thus, this study also considers the AR (7) model. The model parameters are estimated using the maximum-likelihood estimation (MLE) method.A flowchart of the proposed modeling procedure is given in Figure5.Figure 5
Flowchart of the proposed modeling framework.


## 3.1. Moving Filter on Demand

The presence of outliers may significantly affect the estimation and the forecasting. The outliers are mostly present in the time-series data, making the estimation and forecasting of the time-series model more challenging. In time-series analysis, outliers’ detection and replacement with normal values are the essential step of the data cleaning process that, generally, improves the estimation of parameters and the accuracy of forecasting models. In electricity demand or price data, outliers (also known as extreme values) are generally present due to many factors, including unexpected increased demand, power plant shutdown, and extreme weather conditions. The box plot of hourly electricity demand is plotted in Figure1. One can see that there are some extreme values (outliers) present in the demand series, especially at midnight hours. In general, these extreme values are identified and replaced by many methods [53]. In this research work, they are treated by the moving filter on-demand method suggested by Borovkova and Permana [54].Figure 1
Nord Pool electricity demand: hourly box plots for demand data range from January 1, 2015, to December 31, 2019.The moving filter on-demand (MFD) is a generalization of the standard deviation filter on-demand (SFD) technique. The SFD works on the whole time series, whereas the MFD operates on a rolling window of constant width. At the step, the original time series is divided intow=q/r parts, where r represents the width of the window and q is the total number of observations in the series. Then, MFD is applied to the first window of the original time series. The demand values whose absolute deviation from the estimated mean μ^ is greater than a multiple of estimated standard deviation (SD) are identified as extreme values. The subset of spikes Xto detected with the use of the MFD filter is obtained as follows:(1)Yto=Yt:Yt−μ^≥1.64⋅SDY.The window then movesr observations forward, and the filtering method is repeated in the next window to determine the extreme values. The procedure is repeated until this method treats all w parts. Once the extreme values are identified, they are replaced by the median value of the corresponding window to obtain extreme value free demand series.


## 3.2. General Model

After possible outliers’ identification and replacement with the normal values to obtain the demand series free of outliers, the next step is to model and forecast electricity demand. To this end, the component estimation procedure is used to capture different features of the demand series [55]. Suppose the logarithmic electricity demand series is log Yt,j, where t∈N and j=1,…,24. The demand series is decomposed into a deterministic component denoted by Dt,j and a stochastic component denoted by Xt,j as follows:(2)logYt,j=Dt,j+Xt,j.Here, the deterministic component in (2) is comprised of the complex structure of the seasonal component, such as annual seasonality st,j, weekly seasonality wt,j, bank holiday effect bt,j, and long-trend component ct,j. Mathematically, it can be written as follows:(3)Dt,j=ct,j+st,j+wt,j+bt,j.The deterministic componentDt,j in (3) is estimated through a generalized additive model (GAM). In the case of stochastic component Xt,j, three different models, namely the FAR (p), FARX (p), and AR (p), are used. Once both the components are modeled and forecasted, the one-day-ahead out-of-sample forecast is obtained as follows:(4)Y^t+1,j=expD^t+1,j+X^t+1,j.


## 3.3. Modeling the Deterministic Component

This section will explain the modeling and forecasting of the deterministic component. The generalized additive modeling technique is used to model the long-term trendct,j, annual st,j, weekly wt,j periodicities, and bank holiday effect bt,j. More specifically, the long-run (trend) component ct,j, which is a function of time t, is estimated through the ordinary least-squares (OLS) technique. In case of the yearly seasonality, the series 1,2,…,365,1,2,…,365,…,1,2,…,365 is used and is estimated using the OLS. On the other hand, for weekly periodicities and for the bank holidays, dummy variables are used; i.e., wt,j=∑i=17φict,j with ct,j=1 if t refers to the ith day of the week and zero otherwise and bt,j=∑i=12ωct,j with ct,j=1 if t refers to a bank holiday and zero otherwise. The coefficients φi and ωi are again estimated with the OLS approach. Once the deterministic components are estimated, the one-day-ahead forecast for c^t+1,j=c^t,j and s^t+1,j=s^t,j, as both of these components represent long-term dynamics for our forecast horizon, whereas the forecast of wt and bt is straightforward as both of these are deterministic functions of time or the calendar conditions. Hence, one-day-ahead forecast for the deterministic component can be obtained as follows:(5)D^t+1,j=c^t,j+s^t,j+w^t+1,j+b^t+1,j.To obtain the stochastic component, we subtract the deterministic component fromlogYt,j as follows:(6)X^t,j=logYt,j−D^t,j=logYt,j−c^t,j+s^t,j+w^t,j+b^t,j.The modeling procedure of the stochastic component is discussed in the next section in detail.


## 3.4. Modeling the Stochastic Component

This section describes the estimation of the residual componentXt,j, which is modeled using the functional autoregressive, FAR (p), functional autoregressive with exogenous variables, FARX (p), and the classical univariate autoregressive, AR (p), models. The classical AR (p) is a well-known model; however, before explaining the functional models, some preliminaries are required, which are stated in the next section.


### 3.4.1. Basics of Functional Data Analysis

This section provides some preliminaries that are essential for building the functional models. The stochastic componentXt,j in (2) is transformed into functional data using some known basis functions. Using suitable basis functions, the functional object corresponding to the daily demand profile for the day t is represented as follows:(7)Xtτ=∑k=1Kakγt,kττ∈J,where ak is the constant parameters and γt,kτ is the Fourier basis functions. An example of functional time-series trajectories for 1826 days is given in Figure 2 where each functional trajectory represents a daily demand profile.Figure 2
Nord Pool electricity demand: functional daily trajectoriesXtτt=1,…,1826, where Xtτ is the log transformed functional observations of hourly electricity demand, based on 15 Fourier basis functions.Assume thatXtτ given in (7) is an arbitrary functional time series (FTS) defined on a common probability space Φ,∑,P, where Φ,∑, and P denote the sample space, σ−algebra on Φ, and the probability measure on ∑, respectively. In general, it is assumed that the functions Xτ are elements of the semi-metric and Hilbert space. Furthermore, it is assumed that functions are the elements of the square-integrable functions X∈ℒ20,1 residing in the Hilbert space ℋ satisfying Xt2=∫X2τdτ<∞ with an inner product x,y=∫xτyτdτ,∀x,y∈ℒ20,1. The notation X∈ℒℋυΦ,∑,P is used to indicate for some υ>0, EXυ<∞. Note that υ=1 results in the population mean curve μτ=EXτ, and when υ=2, then nonnegative definite covariance operator is obtained as follows:(8)Cτ,δ=CovXτ,Xδ,Cτ,δ=EXτ−μτXδ−μδ.The covariance operatorCτ,δ in (8) allows the covariance operator of X, denoted by G, which can be expressed as follows:(9)Gαδ=∫10Cτ,δατdτ.Using Mercer’s lemma, there is an orthonormal sequenceαk of continuous functions in ℒ20,1 and a non-increasing sequence βk of positive numbers, such that(10)Cτ,δ=∑k=1∞βkαkταkδ,where the αkτ denotes the kth functional principal component (FPC) and βk denotes the kth eigenvalue in the decreasing order [56]. By the separability of Hilbert spaces, the Karhunen–Loève (KL) expansion of a random process Xτ can be expressed as follows:(11)Xtτ=μt+∑k=1∞λk,tαkτ,where λk,t denotes the kth principal component score (PCS) defined as λk,t=∫Yταkτdτt. The PCSs constitute an uncorrelated sequence of random variables with zero mean and variance βk.The expansion 3.7 facilitates dimension reduction as the firstd terms often provide a good approximation to the infinite sums, and thus, the information contained in Xtτ can be adequately summarized by the d-dimensional vector β1,…,βd. The approximated processes can be defined as follows:(12)Xtτ=∑k=1dλk,tαkτ+ετ,where d denotes the number of retained principal components, and ετ denotes the error function with a mean of zero and a finite variance, containing the FPCs excluded from the retained d components. For more details about the FPC analysis and its practical demonstration, the interested reader may consult Ramsay and Silverman [57]; Shang [58]; and Aue et al. [59].In practice, the mean curveμτ, FPC’s ατ, and PCS’s λ can only be estimated through realizations of a random process. Suppose X1τ,…,XNτ are fully observed FTS, and then, the standard estimators for the mean function μt and covariance operator Cτ,δ are given by the following sample averages:(13)μ^τ=1N∑t=1NXtτ,(14)C^τ,δ=∑k=1∞β^kα^kτα^kδ,where β^1>β^2>⋯≥0 are the sample eigenvalues of C^τ,δ, and α^1τ,α^2τ,… are the corresponding orthogonal sample eigenfunctions. Hörmann and Kokoszka [36] studied the aforementioned estimators and concluded that they are consistent estimators for weekly dependent process. Then, using the KL expansion, the realizations of the random process X can be written as follows:(15)Xtτ=μ^τ+∑k=1dλ^k,tα^kτ+eτ,t=1,2…,N,where λ^k,t is the kth estimated PCS for the tth observation. Figure 3 provides the cumulative proportion of variation explained by the first three PCs. Note that the first three PCs explained more than 99% of the total variation. Figure 4 represents the first three functional principal components with the amount of variation explained by each FPC. Note that the 1st eigenfunction explains approximately 95.2% of the total variation in the data.Figure 3
Cumulative proportion of variation explained by the first three PCs.Figure 4
First three functional principal components with amount of variation explained by each FPC.The next section provides details about FAR (p) and FARX (p) models used to model and forecast Xtτ.


### 3.4.2. Functional Autoregressive Model of Order p

The forecasting of the FTS is a complicated task for which several techniques are discussed in the literature. Functional autoregressive (FAR) models are one of the most popular functional models used for FTS forecasting. Consider that the functional trajectoriesXtτ are available for t=1…,N. In the FAR modeling, Xtτ is linearly dependent on its p lagged (past) values and an error term. Mathematically, the FAR (p) model can be written as follows:(16)Xtτ=μτ+∑l=1pΨiXt−lτ+ξtτ,where ξtτ∈N is independently and identically distributed sequence in ℒℋ2 with Eξt=0 and Ψll=1p is bounded linear operator mapping ℒℋ2⟶ℒℋ2 such that the above recursive equation has a unique causal solution. To model and forecast Xtτ in (16) using the FAR(p) model, (13)–(15) are used for modeling and forecasting in the following three steps.(i)
The dimensiond is fixed and the estimated FPC scores are obtained as λ^k,t=∫X^tτα^kτdτ for each observation X^tτ,t=1,…,N,k=1,…,d, and the estimated k-variate FPC score vectors λ^t=λ^1,t,…,λ^d,t′,t=1,…,N.(ii)
For a fixed orderp, the vector autoregressive model (VAR(p)) Yt=∑l=1pΨlYt−l+εt for eigenscore vectors is constructed to produce forecasting λ^N+1=λ^1,N+1,…,λ^d,N+1′. Given the vectors λ^1,…,λ^N, the Durbin–Levinson and innovations algorithm can be readily applied here.(iii)
In the last step, the multivariate time series is re-transformed to functional version using the KL theoremX^N+1τ=μ^τ+λ^1,N+1α^1τ+⋯+λ^d,N+1α^dτ. Based on the estimated FPC scores and sample eigenfunctions, the resulting X^N+1τ is then used as one-step-ahead forecast of XN+1τ.


### 3.4.3. Forecasting with Exogenous Variables FARX (p)

In the FTS modeling, not only the lagged values of a response variable but also other explanatory exogenous variables are added to the model to increase the forecasting accuracy. These exogenous variables may be scalars, vector-valued, or functional. The objective is then to obtain an empiricalX^N+1 given observations of the trajectories X1,…,XN and the number of exogenous variables Xf1,…,Xfv. The predictor variables need not be defined on the same space. For example, Xf1 could be vector, Xf2 a scalar, Xf3t a function, and Xf4t could contain lagged values of Xf3. The FAR model having functional type of exogenous variables is generally denoted by FARX(p). The functional exogenous variables Xf1t,…,Xfρt are assumed to bestationary, with time-invariant mean functions μf1t,…,μfρt. Suppose Xtτ, t∈N, and the exogenous variables Xfρt, f∈ℕ, are available. Then, the FARX(p) model can be defined as [60]follows:(17)Xtτ−μt=∑l=1pΨlτXt−lτ−μτ+∑m=1ρϕmτXfmτ−μfmt+εtτ,where Ψlτ and ϕmτ are the functional operators of the model and εtτ is the functional error term. It is worth mentioning that the dynamics of the functional response can be represented not only with the lagged functional exogenous variables, but there may also be some causal dependence with scalar exogenous variables. Suppose that ρ scalar exogenous variables denoted by Xf1,…,Xfρ having mean μf1,…,μfρ are available. Then, the FARX(p) model can be written as follows:(18)Ytτ−μt=∑l=1pΨlτXt−lτ−μτ+∑m=1ρϕmXfm−μfm+εtτ.Both FARX(p) models in (17) or (18) corresponding to having functional or scalar exogenous variables, respectively, are the generalization of the model discussed in (16) with different numbers of functional or scalar exogenous variables. The estimation method discussed in Section 3.4.2 is modified and summarized as follows:(1)
(a)
The dimensiond is fixed for t=1,…,N and the data X1,…,Xt are used to compute vectors λ^t=λ^1,t,…,λ^d,t′ that contain the first d estimated FPC scores.(b)
In case of functional exogenous variables, the value ofρ is fixed, and the data X1,…,Xtt=1,…,N are used to compute the vector φ^t=φ^1,t,…,φ^v,t′, which contain the first ρ leading FPC scores. The process for each functional exogenous variable is repeated.(c)
In the next step, all the exogenous variable vectors are combined into a single vector,θ^N=θ^1,N,…,θ^ρ,N′.(2)
λ^1,…,λ^N and θ^N are used to obtain one-step-ahead forecast as follows:(19)λ^t+1=λ^1,t+1,…,λ^d,t+1′.(3)
The KL expansion is used to obtain(20)X^t+1τ=μ^τ+λ^1,t+1α^1τ+⋯+λ^d,t+1α^dτ,as one-day-ahead forecast for Xt+1 using the KL theorem.


### 3.4.4. Dimension and Order Selection of the Functional Models

The selection of dimensiond and order p is an important step for the FTS forecasting using FAR (p) and FARX (p). In the literature, various methods are used to select the order of FAR. For example, Kokoszka and Reimherr [61] proposed a multistage testing hypothesis method for the selection of the order of FAR. Salmerón and Ruiz-Medina [62] studied the multispectral decomposition to select the order of the FAR (p) model. One common approach is choosing a minimum number of PCs so that the proportion of variance explained by the FPCs exceeds a certain threshold value. In this study, we used a completely automatic method for choosing order p and dimension d for both FAR (p) and FARX (p), which minimizes the mean square error of forecasting studied [63], which is asymptotically equivalent to the functional final prediction error (fFPE) given as follows:(21)fFPEp,d=t+p∗dt−p∗dtrΔ^Y+∑k>dβ^k,where βk is the kth eigenvalue of Cτ,δ and Δ^Y in (21) is the covariance matrix of the random vector Y1,…,Yd and is an unbiased estimator of ΔY. The best d and p are the minimizer of the fFPE function.


### 3.4.5. Univariate Autoregressive Model

Autoregressive models are the most commonly employed linear models in univariate time-series forecasting problems. In the autoregressive model, the dependent variable is regressed on itsp lagged values. The number of lag values required to forecast the future value is the order of an AR model, which is denoted by the value p. Mathematically, the autoregressive model of order p, AR (p), can be written as follows:(22)Xt=ϕ+Ψ1Xt−1+⋯+ΨpXt−p+εt,where ϕ is the constant (intercept), Ψll=1,…,p is the parameters of AR (p) model, and εt is assumed to follow a white noise process with mean zero and a constant variance. In the literature for demand modeling, an AR (7) is often used. Thus, this study also considers the AR (7) model. The model parameters are estimated using the maximum-likelihood estimation (MLE) method.A flowchart of the proposed modeling procedure is given in Figure5.Figure 5
Flowchart of the proposed modeling framework.


## 3.4.1. Basics of Functional Data Analysis

This section provides some preliminaries that are essential for building the functional models. The stochastic componentXt,j in (2) is transformed into functional data using some known basis functions. Using suitable basis functions, the functional object corresponding to the daily demand profile for the day t is represented as follows:(7)Xtτ=∑k=1Kakγt,kττ∈J,where ak is the constant parameters and γt,kτ is the Fourier basis functions. An example of functional time-series trajectories for 1826 days is given in Figure 2 where each functional trajectory represents a daily demand profile.Figure 2
Nord Pool electricity demand: functional daily trajectoriesXtτt=1,…,1826, where Xtτ is the log transformed functional observations of hourly electricity demand, based on 15 Fourier basis functions.Assume thatXtτ given in (7) is an arbitrary functional time series (FTS) defined on a common probability space Φ,∑,P, where Φ,∑, and P denote the sample space, σ−algebra on Φ, and the probability measure on ∑, respectively. In general, it is assumed that the functions Xτ are elements of the semi-metric and Hilbert space. Furthermore, it is assumed that functions are the elements of the square-integrable functions X∈ℒ20,1 residing in the Hilbert space ℋ satisfying Xt2=∫X2τdτ<∞ with an inner product x,y=∫xτyτdτ,∀x,y∈ℒ20,1. The notation X∈ℒℋυΦ,∑,P is used to indicate for some υ>0, EXυ<∞. Note that υ=1 results in the population mean curve μτ=EXτ, and when υ=2, then nonnegative definite covariance operator is obtained as follows:(8)Cτ,δ=CovXτ,Xδ,Cτ,δ=EXτ−μτXδ−μδ.The covariance operatorCτ,δ in (8) allows the covariance operator of X, denoted by G, which can be expressed as follows:(9)Gαδ=∫10Cτ,δατdτ.Using Mercer’s lemma, there is an orthonormal sequenceαk of continuous functions in ℒ20,1 and a non-increasing sequence βk of positive numbers, such that(10)Cτ,δ=∑k=1∞βkαkταkδ,where the αkτ denotes the kth functional principal component (FPC) and βk denotes the kth eigenvalue in the decreasing order [56]. By the separability of Hilbert spaces, the Karhunen–Loève (KL) expansion of a random process Xτ can be expressed as follows:(11)Xtτ=μt+∑k=1∞λk,tαkτ,where λk,t denotes the kth principal component score (PCS) defined as λk,t=∫Yταkτdτt. The PCSs constitute an uncorrelated sequence of random variables with zero mean and variance βk.The expansion 3.7 facilitates dimension reduction as the firstd terms often provide a good approximation to the infinite sums, and thus, the information contained in Xtτ can be adequately summarized by the d-dimensional vector β1,…,βd. The approximated processes can be defined as follows:(12)Xtτ=∑k=1dλk,tαkτ+ετ,where d denotes the number of retained principal components, and ετ denotes the error function with a mean of zero and a finite variance, containing the FPCs excluded from the retained d components. For more details about the FPC analysis and its practical demonstration, the interested reader may consult Ramsay and Silverman [57]; Shang [58]; and Aue et al. [59].In practice, the mean curveμτ, FPC’s ατ, and PCS’s λ can only be estimated through realizations of a random process. Suppose X1τ,…,XNτ are fully observed FTS, and then, the standard estimators for the mean function μt and covariance operator Cτ,δ are given by the following sample averages:(13)μ^τ=1N∑t=1NXtτ,(14)C^τ,δ=∑k=1∞β^kα^kτα^kδ,where β^1>β^2>⋯≥0 are the sample eigenvalues of C^τ,δ, and α^1τ,α^2τ,… are the corresponding orthogonal sample eigenfunctions. Hörmann and Kokoszka [36] studied the aforementioned estimators and concluded that they are consistent estimators for weekly dependent process. Then, using the KL expansion, the realizations of the random process X can be written as follows:(15)Xtτ=μ^τ+∑k=1dλ^k,tα^kτ+eτ,t=1,2…,N,where λ^k,t is the kth estimated PCS for the tth observation. Figure 3 provides the cumulative proportion of variation explained by the first three PCs. Note that the first three PCs explained more than 99% of the total variation. Figure 4 represents the first three functional principal components with the amount of variation explained by each FPC. Note that the 1st eigenfunction explains approximately 95.2% of the total variation in the data.Figure 3
Cumulative proportion of variation explained by the first three PCs.Figure 4
First three functional principal components with amount of variation explained by each FPC.The next section provides details about FAR (p) and FARX (p) models used to model and forecast Xtτ.


## 3.4.2. Functional Autoregressive Model of Order p

The forecasting of the FTS is a complicated task for which several techniques are discussed in the literature. Functional autoregressive (FAR) models are one of the most popular functional models used for FTS forecasting. Consider that the functional trajectoriesXtτ are available for t=1…,N. In the FAR modeling, Xtτ is linearly dependent on its p lagged (past) values and an error term. Mathematically, the FAR (p) model can be written as follows:(16)Xtτ=μτ+∑l=1pΨiXt−lτ+ξtτ,where ξtτ∈N is independently and identically distributed sequence in ℒℋ2 with Eξt=0 and Ψll=1p is bounded linear operator mapping ℒℋ2⟶ℒℋ2 such that the above recursive equation has a unique causal solution. To model and forecast Xtτ in (16) using the FAR(p) model, (13)–(15) are used for modeling and forecasting in the following three steps.(i)
The dimensiond is fixed and the estimated FPC scores are obtained as λ^k,t=∫X^tτα^kτdτ for each observation X^tτ,t=1,…,N,k=1,…,d, and the estimated k-variate FPC score vectors λ^t=λ^1,t,…,λ^d,t′,t=1,…,N.(ii)
For a fixed orderp, the vector autoregressive model (VAR(p)) Yt=∑l=1pΨlYt−l+εt for eigenscore vectors is constructed to produce forecasting λ^N+1=λ^1,N+1,…,λ^d,N+1′. Given the vectors λ^1,…,λ^N, the Durbin–Levinson and innovations algorithm can be readily applied here.(iii)
In the last step, the multivariate time series is re-transformed to functional version using the KL theoremX^N+1τ=μ^τ+λ^1,N+1α^1τ+⋯+λ^d,N+1α^dτ. Based on the estimated FPC scores and sample eigenfunctions, the resulting X^N+1τ is then used as one-step-ahead forecast of XN+1τ.


## 3.4.3. Forecasting with Exogenous Variables FARX (p)

In the FTS modeling, not only the lagged values of a response variable but also other explanatory exogenous variables are added to the model to increase the forecasting accuracy. These exogenous variables may be scalars, vector-valued, or functional. The objective is then to obtain an empiricalX^N+1 given observations of the trajectories X1,…,XN and the number of exogenous variables Xf1,…,Xfv. The predictor variables need not be defined on the same space. For example, Xf1 could be vector, Xf2 a scalar, Xf3t a function, and Xf4t could contain lagged values of Xf3. The FAR model having functional type of exogenous variables is generally denoted by FARX(p). The functional exogenous variables Xf1t,…,Xfρt are assumed to bestationary, with time-invariant mean functions μf1t,…,μfρt. Suppose Xtτ, t∈N, and the exogenous variables Xfρt, f∈ℕ, are available. Then, the FARX(p) model can be defined as [60]follows:(17)Xtτ−μt=∑l=1pΨlτXt−lτ−μτ+∑m=1ρϕmτXfmτ−μfmt+εtτ,where Ψlτ and ϕmτ are the functional operators of the model and εtτ is the functional error term. It is worth mentioning that the dynamics of the functional response can be represented not only with the lagged functional exogenous variables, but there may also be some causal dependence with scalar exogenous variables. Suppose that ρ scalar exogenous variables denoted by Xf1,…,Xfρ having mean μf1,…,μfρ are available. Then, the FARX(p) model can be written as follows:(18)Ytτ−μt=∑l=1pΨlτXt−lτ−μτ+∑m=1ρϕmXfm−μfm+εtτ.Both FARX(p) models in (17) or (18) corresponding to having functional or scalar exogenous variables, respectively, are the generalization of the model discussed in (16) with different numbers of functional or scalar exogenous variables. The estimation method discussed in Section 3.4.2 is modified and summarized as follows:(1)
(a)
The dimensiond is fixed for t=1,…,N and the data X1,…,Xt are used to compute vectors λ^t=λ^1,t,…,λ^d,t′ that contain the first d estimated FPC scores.(b)
In case of functional exogenous variables, the value ofρ is fixed, and the data X1,…,Xtt=1,…,N are used to compute the vector φ^t=φ^1,t,…,φ^v,t′, which contain the first ρ leading FPC scores. The process for each functional exogenous variable is repeated.(c)
In the next step, all the exogenous variable vectors are combined into a single vector,θ^N=θ^1,N,…,θ^ρ,N′.(2)
λ^1,…,λ^N and θ^N are used to obtain one-step-ahead forecast as follows:(19)λ^t+1=λ^1,t+1,…,λ^d,t+1′.(3)
The KL expansion is used to obtain(20)X^t+1τ=μ^τ+λ^1,t+1α^1τ+⋯+λ^d,t+1α^dτ,as one-day-ahead forecast for Xt+1 using the KL theorem.


## 3.4.4. Dimension and Order Selection of the Functional Models

The selection of dimensiond and order p is an important step for the FTS forecasting using FAR (p) and FARX (p). In the literature, various methods are used to select the order of FAR. For example, Kokoszka and Reimherr [61] proposed a multistage testing hypothesis method for the selection of the order of FAR. Salmerón and Ruiz-Medina [62] studied the multispectral decomposition to select the order of the FAR (p) model. One common approach is choosing a minimum number of PCs so that the proportion of variance explained by the FPCs exceeds a certain threshold value. In this study, we used a completely automatic method for choosing order p and dimension d for both FAR (p) and FARX (p), which minimizes the mean square error of forecasting studied [63], which is asymptotically equivalent to the functional final prediction error (fFPE) given as follows:(21)fFPEp,d=t+p∗dt−p∗dtrΔ^Y+∑k>dβ^k,where βk is the kth eigenvalue of Cτ,δ and Δ^Y in (21) is the covariance matrix of the random vector Y1,…,Yd and is an unbiased estimator of ΔY. The best d and p are the minimizer of the fFPE function.


## 3.4.5. Univariate Autoregressive Model

Autoregressive models are the most commonly employed linear models in univariate time-series forecasting problems. In the autoregressive model, the dependent variable is regressed on itsp lagged values. The number of lag values required to forecast the future value is the order of an AR model, which is denoted by the value p. Mathematically, the autoregressive model of order p, AR (p), can be written as follows:(22)Xt=ϕ+Ψ1Xt−1+⋯+ΨpXt−p+εt,where ϕ is the constant (intercept), Ψll=1,…,p is the parameters of AR (p) model, and εt is assumed to follow a white noise process with mean zero and a constant variance. In the literature for demand modeling, an AR (7) is often used. Thus, this study also considers the AR (7) model. The model parameters are estimated using the maximum-likelihood estimation (MLE) method.A flowchart of the proposed modeling procedure is given in Figure5.Figure 5
Flowchart of the proposed modeling framework.


## 4. Analysis of Nord Pool Electricity Demand

This section starts with a description of the Nordic power market and provides details about the data used in this study. Some of the features of the data set used are also given in this section. In addition, different accuracy descriptive statistics used in the study are also introduced in this section.


### 4.1. The Nordic Power Market

The Nordic electricity market includes Finland, Norway, Sweden, and Denmark, which deregulated and integrated their respective electricity markets during the 1990s. The data set consists of hourly electricity demand data obtained from the Nord Pool electricity market (https://www.nordpoolgroup.com/). The electricity markets operate 24 hours a day, seven days a week. The performance of the proposed model is validated with the electricity demand time series ranging from January 1, 2015, to December 31, 2019, corresponding to 43,824 equidistant time points (hourly load periods). For modeling and forecasting purposes, the data set is split into two parts:(i)
Model Estimation Period: from January 1, 2015, to December 31, 2018 (35,064 observations, coveringt  = 1461 days)(ii)
Out-of-Sample Forecast Period: from January 1, 2019, till December 31, 2019 (8760 observations, coveringt  = 365 daysThe models are estimated based on expending window technique, and the one-day-ahead out-of-sample forecast is obtained for the whole year.As described in Section1, the yearly periodicity refers to the cyclical demand behavior from one year to another that is evident in the demand time series. For example, a discrete-time series display of electricity demand is shown in Figure 6. The yearly seasonality is evident from the graph as the electricity demand pattern over a year is similar to the adjacent year. The demand is higher in winter compared with other seasons. The load profile on working days shows different behaviors than the nonworking days resulting in weekly periodicity. The electricity demand daily profiles for a week are plotted in Figure 7, showing demand variation on working and nonworking days.Figure 6
Nord Pool electricity demand: a univariate time-series display of electricity demand for the period January 1, 2015, to December 31, 2019. There are 43,824 discrete time points each representing one load hour.Figure 7
Nord Pool electricity demand: curves represent demand functional profiles for one week.


### 4.2. Forecast Accuracy Measures

Different forecast accuracy measures are computed to compare the models. To simplify the comparison, we convert the forecasted trajectories into hourly electricity demand data. The forecast accuracy is assessed by the mean absolute error (MAE), mean absolute percentage error (MAPE), and root-mean-square error (RMSE).


#### 4.2.1. Mean Absolute Error

The MAE or mean absolute deviation (MAD) is calculated by taking the average of absolute difference between the estimated forecast and the actual value at the same time so that the negative values do not cancel the positive values. Mathematically, it can be written as follows:(23)MAE=meanYt,j−Y^t,j,where Yt,j is the observed (the original observed series) and Y^t,j is the forecasted one-day-ahead electricity demand for ttht=1,2,…,365 day and jthj=1,2,…,24 load period.


#### 4.2.2. Mean Absolute Percentage Error

The MAPE is calculated using the average of the absolute error in each period divided by the observed values multiplied by 100. The MAPE is useful when the size or size of a prediction variable is significant in evaluating the accuracy of forecasting [64]. The MAPE is a relative error indicating how much error is observed in forecast compared with the observed. Mathematically, the MAPE can be expressed as follows:(24)MAPE=meanYt,j−Y^t,jYt,j×100.


#### 4.2.3. Root-Mean-Square Error

The RMSE depends on the scale of the dependent variable. It is used as a measure to compare forecasts for the same series across different models. The RMSE is the square root of the average squared differences between predicted and actual values. Mathematically, it is defined as follows:(25)RMSE=meanYt,j−Y^t,j2.


## 4.1. The Nordic Power Market

The Nordic electricity market includes Finland, Norway, Sweden, and Denmark, which deregulated and integrated their respective electricity markets during the 1990s. The data set consists of hourly electricity demand data obtained from the Nord Pool electricity market (https://www.nordpoolgroup.com/). The electricity markets operate 24 hours a day, seven days a week. The performance of the proposed model is validated with the electricity demand time series ranging from January 1, 2015, to December 31, 2019, corresponding to 43,824 equidistant time points (hourly load periods). For modeling and forecasting purposes, the data set is split into two parts:(i)
Model Estimation Period: from January 1, 2015, to December 31, 2018 (35,064 observations, coveringt  = 1461 days)(ii)
Out-of-Sample Forecast Period: from January 1, 2019, till December 31, 2019 (8760 observations, coveringt  = 365 daysThe models are estimated based on expending window technique, and the one-day-ahead out-of-sample forecast is obtained for the whole year.As described in Section1, the yearly periodicity refers to the cyclical demand behavior from one year to another that is evident in the demand time series. For example, a discrete-time series display of electricity demand is shown in Figure 6. The yearly seasonality is evident from the graph as the electricity demand pattern over a year is similar to the adjacent year. The demand is higher in winter compared with other seasons. The load profile on working days shows different behaviors than the nonworking days resulting in weekly periodicity. The electricity demand daily profiles for a week are plotted in Figure 7, showing demand variation on working and nonworking days.Figure 6
Nord Pool electricity demand: a univariate time-series display of electricity demand for the period January 1, 2015, to December 31, 2019. There are 43,824 discrete time points each representing one load hour.Figure 7
Nord Pool electricity demand: curves represent demand functional profiles for one week.


## 4.2. Forecast Accuracy Measures

Different forecast accuracy measures are computed to compare the models. To simplify the comparison, we convert the forecasted trajectories into hourly electricity demand data. The forecast accuracy is assessed by the mean absolute error (MAE), mean absolute percentage error (MAPE), and root-mean-square error (RMSE).


### 4.2.1. Mean Absolute Error

The MAE or mean absolute deviation (MAD) is calculated by taking the average of absolute difference between the estimated forecast and the actual value at the same time so that the negative values do not cancel the positive values. Mathematically, it can be written as follows:(23)MAE=meanYt,j−Y^t,j,where Yt,j is the observed (the original observed series) and Y^t,j is the forecasted one-day-ahead electricity demand for ttht=1,2,…,365 day and jthj=1,2,…,24 load period.


### 4.2.2. Mean Absolute Percentage Error

The MAPE is calculated using the average of the absolute error in each period divided by the observed values multiplied by 100. The MAPE is useful when the size or size of a prediction variable is significant in evaluating the accuracy of forecasting [64]. The MAPE is a relative error indicating how much error is observed in forecast compared with the observed. Mathematically, the MAPE can be expressed as follows:(24)MAPE=meanYt,j−Y^t,jYt,j×100.


### 4.2.3. Root-Mean-Square Error

The RMSE depends on the scale of the dependent variable. It is used as a measure to compare forecasts for the same series across different models. The RMSE is the square root of the average squared differences between predicted and actual values. Mathematically, it is defined as follows:(25)RMSE=meanYt,j−Y^t,j2.


## 4.2.1. Mean Absolute Error

The MAE or mean absolute deviation (MAD) is calculated by taking the average of absolute difference between the estimated forecast and the actual value at the same time so that the negative values do not cancel the positive values. Mathematically, it can be written as follows:(23)MAE=meanYt,j−Y^t,j,where Yt,j is the observed (the original observed series) and Y^t,j is the forecasted one-day-ahead electricity demand for ttht=1,2,…,365 day and jthj=1,2,…,24 load period.


## 4.2.2. Mean Absolute Percentage Error

The MAPE is calculated using the average of the absolute error in each period divided by the observed values multiplied by 100. The MAPE is useful when the size or size of a prediction variable is significant in evaluating the accuracy of forecasting [64]. The MAPE is a relative error indicating how much error is observed in forecast compared with the observed. Mathematically, the MAPE can be expressed as follows:(24)MAPE=meanYt,j−Y^t,jYt,j×100.


## 4.2.3. Root-Mean-Square Error

The RMSE depends on the scale of the dependent variable. It is used as a measure to compare forecasts for the same series across different models. The RMSE is the square root of the average squared differences between predicted and actual values. Mathematically, it is defined as follows:(25)RMSE=meanYt,j−Y^t,j2.


## 5. Results and Discussion

This section provides details about the out-of-sample forecasting performance for various models described above. In addition, the results from the empirical study are discussed in this section.The main objective of the current research study is to model and forecast electricity demand using functional models assuming that the dynamics of electricity demand can be decomposed into deterministic and stochastic components. For this purpose, in the1st step, the outliers are detected and replaced with normal values by the method discussed in Section 3.1. In the 2nd step, log transformation (logarithms) is applied to the filtered demand series to stabilize the variance of the series. The filtered log demand series is then divided into deterministic and stochastic components. Both of these components are modeled and forecasted as described in Sections 3.3 and 3.4, and one-day-ahead forecasts are obtained for a whole year (365 days) using the rolling window technique.The accuracy of the aforementioned models is measured by the MAE, MAPE, and RMSE values. In Table1, the overall accuracy of the understudy models is listed. From this table, it is evident that functional autoregressive models outperformed the classical univariate autoregressive model. Within functional models, it is noted that the accuracy of the functional model is further improved by the inclusion of the discrete exogenous (daily maximum and minimum demand) variables in the FAR (p) model. The FARX (p) produces lower forecasting errors within the functional models than FAR (p), suggesting that the exogenous variable improves the forecasting. The FARX (p) model obtains the lowest MAPE value of 2.736, whereas FAR (P) and AR (P) produce significantly higher MAPE values of 6.274 and 9.732, respectively. The MAE and RMSE are also lower for FARX (p) and comparatively higher for the other two competitors. Overall, the MAPE values are small, which suggests the usefulness of our proposed forecasting methodology.Table 1
Out-of-sample forecasting: performance comparison among FAR (p), FARX (p), and AR (7) in terms of MAE, MAPE, and RMSE.
ModelMAEMAPERMSEFAR (p)1525.3546.2741918.528FARX (p)1129.8342.7361552.573AR (7)1640.4149.7322163.824The monthly forecasting errors for the models used in this study are listed in Table2. From this table, it is noted that the errors are relatively higher in December, January, and February compared with other months. As the electricity demand increases in these months, it is natural to have higher forecasting errors in these months compared with other months. Again, the functional models perform relatively well compared with the classical AR (p) model. The forecasting performance of the FARX (p) model is evident from this table as it produces comparatively lower forecasting errors than the rest. The lowest and highest MAPE values obtained by the FARX (p) models are 2.189 and 3.51. These MAPE values correspond to August and January, respectively, indicating that lower demand can be more accurately forecasted than higher demand.Table 2
Monthly forecasting errors for the electricity demand using FAR (p), FARX (p), and AR (7) models.
MonthErrorsFAR (p)FARX (p)AR (7)JanuaryMAE1976.2771706.8312143.758MAPE4.0133.5147.303RMSE2571.6412382.3792873.344FebruaryMAE1654.441200.222055.76MAPE3.4972.5235.461RMSE1940.0571474.2272368.534MarchMAE1455.7081269.2831717.423MAPE3.2762.8045.292RMSE1843.2771694.0252818.67AprilMAE1675.4771352.4912081.115MAPE5.3383.4356.442RMSE1802.1991223.6852812.542MayMAE1485.3931099.2511715.463MAPE7.5052.9929.153RMSE1717.731513.1662055.535JuneMAE1492.786996.8071957.983MAPE8.0692.76411.124RMSE1462.9051263.8671795.036JulyMAE1359.679817.0732292.06MAPE8.2082.34310.539RMSE1496.4671074.0371853.465AugustMAE1430.913752.5071876.894MAPE7.6462.18911.899RMSE1322.593959.5752393.628SeptemberMAE1305.892821.4951789.749MAPE6.9122.48110.713RMSE1584.8881066.2272709.673OctoberMAE1332.144950.9211833.398MAPE5.4022.3417.256RMSE1991.1671215.6942765.328NovemberMAE1374.9631017.9361874.082MAPE3.1092.2415.673RMSE1674.7141291.181874.082DecemberMAE1848.7971569.3362081.957MAPE5.0513.4059.978RMSE2314.8472133.822987.564Table3 provides the weekly forecasting errors for the models under study. The weekly errors are generally higher for Monday, Saturday, and Sunday compared with other days of the week. The FARX (p) produces the lowest MAPE value of 2.36, which corresponds to Thursday, for which the electricity demand is considered to be more stable than the rest of the days. On the other hand, the higher MAPE value of 3.017 is observed on Monday, the day after the weekend. On Monday, the demand profile changes compared with the weekend, and hence, forecasting accuracy lowers for this day. As in previous tables, the functional model, especially FARX (p), outperformed the univariate AR (p) model.Table 3
Weekly forecasting errors for the electricity demand using FAR (p), FARX (p), and AR (7) models.
ModelErrorDays of a weekMTWTFSSFAR (p)MAE1765.2161447.9861551.1051356.9511683.9062489.7252291.673FARX (p)1269.8481164.7551190.67991.0721239.6341108.771943.515AR (7)2793.4452228.0322102.742051.9572688.7743951.4733394.153FAR (p)MAPE5.6915.6495.7105.2175.2737.6647.657FARX (p)3.0172.72152.77162.36292.92502.88742.4657AR (7)8.4587.3008.0729.85810.18912.89410.518FAR (p)RMSE2389.492359.1591987.3792245.9292134.5621743.841978.482FARX (p)1742.8631782.1861640.1261349.3371628.8511397.7791235.803AR (7)4116.0574081.3433825.0073535.034231.546420.9285852.956Finally, the hourly forecasting errors for FAR (p), FARX (p), and AR (p) in terms of the MAE, MAPE, and RMSE are listed in Table 4. Note that the errors vary throughout the day as peak and low hour demand forecasting results can be different. For example, the lowest MAPE value of 1.837 is observed for the second load period, while a MAPE value of 3.82 is the highest value obtained from the FARX (p) model for the seventeenth load period. The superior performance of the functional models is evident as their forecasting accuracy is considerably higher than the classical AR (p) model.Table 4
Hourly forecasting errors for the electricity demand using FAR (p), FARX (p), and AR (7) models.
HourErrorsFAR (p)FARX (p)AR (7)1MAE959.445703.9081148.178MAPE4.6221.9626.883RMSE1298.152896.1261570.792MAE882.311651.4841096.252MAPE4.30781.8378.107RMSE1157.836856.8421586.9413MAE972.551694.9331079.363MAPE5.3351.96910.089RMSE1231.743919.4931563.4314MAE1239.432723.0641813.955MAPE6.2482.04410.190RMSE1271.563966.1331912.9505MAE1006.819770.7351518.308MAPE7.2692.1468.461RMSE1484.1611122.1341545.1846MAE907.918877.4361068.968MAPE7.8963.37610.272RMSE1643.7871162.1241980.2387MAE1432.7251198.6241972.094MAPE7.1193.0159.938RMSE1786.4851630.9142363.8818MAE1693.7341428.4192078.932MAPE6.1743.47.917RMSE2299.9921677.4272824.999MAE1870.8951380.5342068.726MAPE5.1773.2146.686RMSE1918.8161317.0622392.00410MAE1727.9771303.5221995.832MAPE4.5773.0076.538RMSE2024.7481791.1352576.14711MAE1641.9291278.961937.907MAPE4.8052.9525.295RMSE1931.6041332.2592379.14812MAE1712.4081279.8192070.237MAPE4.7522.9698.598RMSE1956.9121699.012492.60313MAE1504.9711250.9541835.173MAPE4.6842.925.283RMSE1983.9861642.1442452.33714MAE1689.1841240.9221910.258MAPE4.8772.9036.259RMSE2106.1951657.0132541.69115MAE1369.6191026.7131574.099MAPE5.892.8749.283RMSE1778.6691373.8662319.1216MAE1445.2841020.281737.731MAPE4.7282.8527.031RMSE1768.3461650.8322249.52217MAE1340.768922.1931711.119MAPE5.7313.825.927RMSE1621.0131123.7762061.40918MAE1401.9591204.1521662.161MAPE6.1442.7586.775RMSE1812.9231566.0422338.13419MAE1314.4331151.8991657.877MAPE5.9522.626.732RMSE1790.0241564.7452149.79520MAE1405.251227.0781648.65MAPE3.9552.7535.77RMSE2213.6711685.0732643.68821MAE1526.5281292.2261989.541MAPE4.2732.9527.782RMSE2273.4981773.1762931.27122MAE1540.7231316.0391971.014MAPE4.4323.097.624RMSE1954.121765.1612354.60923MAE1721.7271280.9232252.269MAPE4.8293.1246.978RMSE1851.8481686.1582360.29524MAE1820.8981391.1922567.675MAPE7.793.10210.543RMSE1781.0211532.952145.987


## 6. Conclusions

Electricity demand forecasting in the deregulated electricity market is a very important and challenging task. The forecasting problem is not straightforward as the demand time series contains many specific features, including high frequency, long trend, multiple seasonalities, spikes or jumps, and bank holiday effect. To accurately forecast electricity demand, these features must be addressed by the forecasting model. Therefore, this study revisited the problem of one-day-ahead electricity demand forecasting using functional modeling techniques. To this end, the demand series is filtered first for the extreme values. The filtered series is then split into two components: deterministic and stochastic. Each component is modeled and forecasted separately, and the final forecasts are obtained by adding them back. The generalized additive modeling technique is used to model the deterministic component that contains the long trend, yearly and weekly seasonalities, and bank holiday effect. For the stochastic component, functional autoregressive (FAR), FAR with exogenous variable (FARX), and classical univariate AR models are used. The exogenous variables include the lagged minimum and maximum electricity demand. Data from the Nord Pool electricity market is used, and the one-day-ahead out-of-sample forecasts are obtained for a whole year. The forecast accuracy is evaluated using the mean absolute error (MAE), mean absolute percentage error (MAPE), and root-mean-square error (RMSE). The results suggest that the functional modeling approach produces significantly better results than the classical approach. The FARX model produces lower forecasting errors compared with the FAR and AR.As the inclusion of the exogenous variable significantly improves the forecasting accuracy, the effect of other exogenous variables, such as daily temperature and season effect, can be investigated in the future. Furthermore, as this study only considers linear model comparison, nonlinear models can also be compared with the proposed functional models.



---

*Source: 6709779-2022-06-11.xml*