Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
88,839 characters · 7 sections · 65 citation commands
Efficient Modeling and Forecasting of Electricity Spot Prices
With the ongoing liberalization of electricity markets over the past decades, the volume of electricity traded via exchanges has greatly increased. This in turn led to an increasing transparency of the price for electricity. Due to the substantial dependence of companies and private households on this price, modeling electricity prices has become one of the cornerstones of research into the energy markets. But such modeling turns out to be at the edge of many disciplines in research. For instance, the analysis of the underlying trade mechanisms can be allocated to economics. But the energy production itself is a process which can be related to engineering and the rules governing the exchange of energy, especially renewable energy production, are determined by law and politics. Hence, modeling electricity prices can be a complex issue. This is also reflected in the time series, where many unusual but already stylized facts can be observed. The model proposed in our paper tackles this complexity in several ways. Its distinctive features compared to the existing literature can be boiled down to seven key facts. Our approach: 1. Models the electricity price without any data manipulation, 2. Incorporates every established stylized fact of electricity prices, 3. Provides insights for the structure of the leverage effect, 4. Proves the effect of wind and solar power on price, 5. Accounts for specific holiday effects and daylight saving time effects in the wind and solar generation, 6. Does not need any future information to provide accurate forecasts. Finally, 7, it uses efficient and rapid state of the art estimation techniques.
We will fit our proposed model to the hourly electricity price of the European Power Exchange (EPEX) for the period of 28.09.2010 up to 01.05.2014. The data sets are obtained from the EPEX at www.epexspot.com for the hourly day-ahead spot price data of Germany/Austria, from the European Network of Transmission System Operators for Electricity at www.entsoe.eu for the hourly load data of Germany, and from the Transparency Platform of the European Energy Exchange (EEX) at www.transparency.eex.com for the hourly wind and solar power feed-in for Germany.\footnote{We remark that the load and renewable data for Austria is not included in the data as our main focus is the German energy market.} Our paper is organized as follows. In Section (ref) we give an overview about the setting in which electricity exchange takes place and name the distinct challenges occurring in modeling the time series. The subsequent section sets up our model, which aims to face those challenges. Section (ref) presents the estimation procedures for our model. In Sections (ref) and (ref) we fit the model to the hourly electricity price time series of the European Power Exchange and apply a comprehensive forecast study. The last section concludes by discussing our findings.
Energy markets are a rapidly changing field of the economy. The liberalization of these markets and the subsequent development of the energy mix account for that fact. But as different countries worldwide face different preconditions, for instance politically or climatically, their energy markets tend to have a very heterogeneous structure. Hence, the findings for one country may not, or may only partly, be used for another country or region. In case of the German energy market, a substantial amount of the daily demand is traded via an exchange. Spot market trading takes place by continuous trading and auctions. Prices for this market can be obtained either from intraday trading or from day-ahead auction prices. The latter is represented by the EPEX spot auction price for Austria and Germany. EPEX is a member of the EEX group. Since 2008, the EEX has not prohibited negative prices keles2012comparison. By considering all products of the exchange, there are currently 250 market participants at the EEX. Considering only the Austrian/German day-ahead auction of the EPEX results in 197 participating traders. The non-negligible amount of market participants is also an aftermath of liberal energy laws in Germany, especially when renewable energy is considered. The Erneuerbare-Energien-Gesetz (EEG) and its corresponding enactments are governmental regulations which embody this idea of liberalization. As a consequence, the supply of energy from renewables rose significantly within the past years. Incentives, like feed-in-tariffs for renewables granted by the government, catalyzed this development. But the growing amount of power plants for renewables had a direct influence on the price for electricity, edenhofer2013economics, as such energy can be produced at a cost of almost zero wurzburg2013renewable. Changes in the EEG in 2012 had direct impact on the marketing of renewable energy.\footnote{Even more recent changes of the EEG in August 2014 are not considered within this article, as our dataset ends in May 2014.} The EEG allows for marketing the produced renewable energy not only to TSOs or other market players but also directly at the exchange. According to \S 33g EEG, producers of such energy receive a market premium in addition to the market price for selling their electricity. Hence, also producers of renewable energy have an incentive to directly trade at the electricity exchange or to sell it via OTC business. Nevertheless, other regulations have also affected the market directly. For instance, transmission system operators (TSO) are obliged to sell their electricity only on the day-ahead or intraday spot market, if they decide to use an exchange.\footnote{According to \S 2 AusglMechV} Hence, data since the publication of this regulation in 2009 are of special interest. In the economic theory of competitive markets, the price of electricity should equal its marginal cost. As the feed-in of renewable energy with zero marginal cost replaces every other energy source with higher marginal cost, the price of electricity should decline keles2013combined. This is known as the merit-order effect. Furthermore, the demand curve for electricity can be assumed to be inelastic sensfuss2008merit as a certain amount of power is needed regardless of the price. This implies that modeling the impact of renewables on the price of electricity is highly necessary, as those energy sources will always lead to a modified marginal cost structure as the share of renewable energy sources is increasing. Empirical evidence for the reduction of electricity prices caused by the emergence of renewable energy has been shown by many authors in the recent literature. For instance, woo2011impact use a regression analysis for the Texas electricity price market to examine the effect of wind power generation. Another multivariate regression approach was applied to German and Austrian electricity prices by wurzburg2013renewable, where, among other time series, wind and solar power were examined and also led to a reduction in price. huisman2013renewable obtained equivalent results for the Nord Pool market by modeling energy supply and demand. This relation is also illustrated in Figure (ref), which shows the price, load, and the patterns of solar power and wind power for two weeks of October 2012. Whenever the combined effect of wind and solar power is high, the line representing the price seems to decline considerably.
But the introduction of renewables not only leads to a price reduction effect, it can also increase the goodness-of-fit for modeling the time series of prices. Concerning this, a comprehensive study was done by cruz2011effect. They combined sophisticated models like Holt--Winters, ARIMA, and neural networks, with, e.g., wind power, and provided evidence for this inclusion's being beneficial to the price modeling. Also yan2013mid, liebl2013modeling and kristiansen2012forecasting obtained competitive goodness-of-fit statistics by including wind power in their approaches. Using a simulation study for the dependencies of the EEX electricity price and the wind power generation keles2013combined were able to show the advantages of including wind power. The price of electricity is also heavily dependent on the day of the week and also changes over a year, which last is primarily governed by the four seasons of the year weron2006modeling. The reason for this is twofold. First, especially solar energy production depends by a law of nature on the period of sunshine: reduced, for instance, in winter. Second, the daily demand for electricity is dependent on the working days, e.g. whether industrial machines are running and require energy or not. An example of two observed months for illustrative purposes is given in Figure (ref). The light-brown shaded area represents the electric load pattern of the first two weeks of October in 2012. Comparing the price time series with the load time series indicates that both variables are positively related. As on weekends the load is usually lowered, weekends and weekdays\footnote{We use the term weekday to refer to the days from Monday to Friday} within the chosen period tend to exhibit different price patterns. This is also true for holidays, e.g. German unity day on the 3rd of October, where the price and the load are both lower than on other weekdays. A more detailed investigation of special days and phases of the day is done in Section (ref). The consideration of these effects is also an important topic in the literature. The weekly dependence is usually modeled in time series analysis by incorporating the equivalent lagged value, e.g. lag 168 for hourly data. (e.g. in kristiansen2012forecasting) Nevertheless, the consideration of special days or special phases of the day is done only rarely. cancelo2008forecasting use a combined ARIMA model to forecast load and provide evidence for the inclusion of special days within the estimation. guthrie2007electricity calculate a periodic autoregression for the half-hourly electricity price of the New Zealand Electricity Market. They divide the day into five different time periods and show that they differ from each other. Another appealing approach is to let the coefficients of the proposed model vary over time. For instance, karakatsani2008forecasting present a comparative analysis of models with and without time varying parameters. They show for British half-hourly electricity prices that models with time varying parameters dominate other autoregressive approaches with constant parameters. A combined model for the inclusion of holidays and time-varying coefficients was introduced by koopman2007periodic. They use a Reg-ARFIMA-GARCH model for the EEX, Powernext and AAX price data. Electricity is also distinct from most other commodities as it is not efficiently storable with the existing technology. This, therefore, adds another source of risk. knittel2005empirical Nevertheless, to a certain extent it can be regarded as indirectly storable because some energy sources such as fuel or gas can be stored and therefore used equivalently to fulfill obligations from, e.g., derivatives. huisman2012electricity The combination of almost non-storability, inelastic demand, and the strong dependence on highly fluctuating energy supplies causes the price of electricity to be endowed with unique characteristics. These are known as the stylized facts of electricity prices. Most authors refer to three or even more characteristics. The most frequently mentioned are seasonality, mean-reversion, and high heteroscedastic volatility with extreme price spikes. (E.g. weron2006modeling and eydeland2003energy.) For illustrative purposes, the EPEX electricity price of the first two months of the year 2012 are depicted in Figure (ref). The pattern shows high price spikes in both directions, which usually last for several hours. After the impact of such a shock, the price reverts to its usual daily level.
Empirical consideration of negative prices in the German/Austrian electricity market is quite rare as they were not present until 2009. Moreover, modeling negative prices is also a difficult task, as some fitting and transformation approaches like logarithmic transformation are not possible. Hence, some authors simply cut off or shift the time series. Nevertheless, keles2012comparison showed via simulation study for the EEX that the inclusion of negative prices yielded better results for their investigated models. Also fanone2013case argue for the inclusion of negative prices. Price spikes can be modeled in different ways. According to christensen2012forecasting there are mainly three different common approaches. These are specific autoregressive, Markov-regime-switching, and jump diffusion models. Other authors take advantage of a combination of these models, e.g. bordignon2013combining and escribano2011modelling. Besides these obvious characteristics of electricity prices, the recent literature has argued for an appreciation of the leverage effect, which is well-known in financial economics. A leverage effect accounts for the variance of a time series's reacting asymmetrically to negative and positive price shocks. black1976leverage The standard leverage effect has higher variance for negative than for positive price shocks. By analyzing hourly California electricity prices knittel2005empirical were able to detect an inverse leverage effect for the time series. Their paper showed that the volatility responses were more intense to positive than to negative shocks. This was later acknowledged for instance by bowden2008short and liu2013applying. cifter2013forecasting, however, provided evidence for the standard leverage effect by examining the daily returns of the Nord Pool electricity price market. As a matter of fact, modeling electricity prices is a complex issue. In addition to classic autoregression approaches, many researchers have applied methods from related scientific fields, e.g. heuristic or a combination of autoregression and other techniques. As a detailed overview of all those approaches would be outside the scope of this paper, the interested reader is referred to the works of weron2006modeling or aggarwal2009electricity who provide a comprehensive overview of the recent literature. In the following section we will address the challenges which occur when electricity prices are examined. We introduce our own model, which is tailor-made for the hourly EPEX electricity price data. It will be shown in detail how the challenges influence our approach and how they are tackled. We will refer to this model as periodic VAR-TARCH. It includes hourly data of the electricity load, and of wind and solar power feed-in. The model is capable of forecasting the price itself and also its random dependent variables.
Let $\boldsymbol Y_t$ be the considered multivariate process with dimension $d=3$ . So $\boldsymbol Y_t = \left(Y_t^{{\mathcal{P}}}, Y_t^{{\mathcal{L}}}, Y_t^{{\mathcal{R}}} \right)$ is a vector of the price of electricity, the load, and the renewable power feed-in at time $t$. Note that $Y_t^{{\mathcal{R}}}$ is the sum of the wind and solar power feed-in. Let ${\mathcal{I}} = \{{\mathcal{P}}, {\mathcal{L}}, {\mathcal{R}} \}$ be the corresponding index set for the three univariate processes.
The basic model is almost equivalent to a simple vector autoregressive model. It is given by the autoregressive model
for $i\in {\mathcal{I}}$, where $\boldsymbol \mu(t) = (\mu^{\mathcal{P}}(t), \mu^{\mathcal{L}}(t), \mu^{\mathcal{R}}(t))$ is a trend component, $\boldsymbol \varepsilon_t = \left({\varepsilon}^{\mathcal{P}}_t, {\varepsilon}^{\mathcal{L}}_t, {\varepsilon}^{\mathcal{R}}_t \right)$ are noise terms with zero mean and time dependent covariance matrix $\Sigma_t$, and they are assumed to be independent. Note that the autoregressive parameters $\phi^{i,j}_k(t)$ may also depend on time, such as the trend $\mu(t)$ and the matrix $\Sigma_t$. The considered lags are given by the index sets $I_{i,j}$. These index sets are crucial and determine the autoregressive dependency structure within $\boldsymbol Y_t$. Here it is important to know that the estimation method described below contains an automatic model selection procedure. Hence, there is always a trade-off in the pre-selection of index sets. Giving the model large index sets $I_{i,j}$ can lead to better models as there is more possible information provided. But on the other side this increases the computation time and the probability of modeling noise. Thus, especially the index sets $I_{i,j}$ have to be chosen carefully.
The index sets $I_{i,j}$ we used are given in Table (ref).
On first view the chosen lags in $I_{i,j}$ seem more or less random, albeit every lag is chosen by employing statistical facts of the data.
For most index sets $I_{i,j}$ the first lags $1,2, \ldots 361 $ simply model the linear dependence of the past 361 hours, so 2 weeks plus a day and an hour. For $I_{{\mathcal{P}}, {\mathcal{R}}}$, $I_{{\mathcal{L}}, {\mathcal{R}}}$ we chose a dependence of the past 49 hours. But we also added some larger lags like $504$ and $505$ to $I_{{\mathcal{P}}, {\mathcal{P}}}$, $I_{{\mathcal{L}},{\mathcal{L}}}$, and $I_{{\mathcal{P}}, {\mathcal{L}}}$ which account for a possible dependence on the price of the past weeks.
These lags allow a structure in the process that is similar to a multiplicative seasonal AR model. For example, the multiplicative structure of the lag polynomials in a seasonal AR(1)$\times$(2)$_{168}$ is given by $(1-a_1B)(1-a_{168}B^{168}- a_{336}B^{336}) = 1- a_1B - a_{168}B^{168} + a_{1}a_{168}B^{169} + a_{336}B^{336} + a_{1}a_{336}B^{337} $ with $B$ as backshift- resp. lag operator.
Thus it contains lag 1, 168, 169, 336, 337. Using this evauation of the lag-polynomial for higher order seasonal AR processes, we receive that lags like 504, 505, 672, 673, 840, 841 or 1008,1009 should be contained in the model. Such seasonal lags were also considered by different authors in the past, such as liu2013applying.
As mentioned above, we could also increase index sets, e.g. replacing $I_{{\mathcal{P}},{\mathcal{P}}}$ by $\{1,\ldots,1009\}$. However, we noticed during our study that many parameters between $(k-1)168+1$ and $k168$ are zero for large $k$, i.e. $k\geq3$. Another simple but useful approach for selecting appropriate index sets is based on evaluating the sample partial autocorrelation function of $\boldsymbol Y_t$; the lags that correspond to large absolute values should be included. Indeed this approach leads to very similar index sets as we use, especially we observe large PACF values at lags $504$ and $505$ as discussed above. The pre-selection of lags therefore provides some flexibility to the applicant, as lags, which are just assumed to be significant can be implemented very easily and will still undergo a statistical reliable lag selection process.
Note that we choose $I_{{\mathcal{P}}, {\mathcal{R}}}$ and $I_{{\mathcal{L}}, {\mathcal{R}}}$ as being small, as we assume that there is no impact of renewable energy feed-in to the price to last for more than two days. Moreover, the renewable energy feed-in does not have any weekly seasonal structure, only a daily and a yearly one. Thus, $I_{{\mathcal{R}},{\mathcal{R}}}$ does not contain higher order weekly based lags like $504$ and $505$ as e.g. $I_{{\mathcal{P}},{\mathcal{P}}}$.
The trend component $\boldsymbol \mu(t)$ is modeled by a linear combination of $M(\mu^i)$ basis functions plus a linear trend, so we have
for $i\in {\mathcal{I}}$ with parameters $\mu^i_j$ and basis functions $\widetilde{B}^{\mu^i}_j(t)$. The structure of the used basis functions is inspired by periodic basis functions. This enables modeling recurring events like specific days of the week. Nevertheless, as some events like holidays are not periodically, we are allowing $\widetilde{B}^{\mu^i}_j(t)$ to be not strictly periodic. As basic concept for the construction of $\widetilde{B}^{\mu^i}_j(t)$ we consider periodic B-splines as used in harvey1993forecasting. However, this approach using local B-splines was mainly chosen to increase the flexibility and the interpretability of the model. We therefore want to point out, that other, e.g. Fourier approaches, might lead to estimation results of comparable quality.
The chosen B-splines are denoted as $B_{T,d_{\mathcal{K}}}(t) =B_{T,d_{\mathcal{K}}}(t; {\mathcal{K}}(d_{\mathcal{K}},T,D), D)$, where $D$ represents the degree of the spline, ${\mathcal{K}} = {\mathcal{K}}(d_{\mathcal{K}},T,D) = \{k_{0},\ldots, k_{D+1}\}$ the set of equidistant knots $ k_{0} < \cdots < k_{D+1} $, that have equal distance $d_{\mathcal{K}}$ and are centered around $T$. Hence, we have $k_{0}= T - d_{\mathcal{K}} \frac{D+1}{2}$, $k_{D+1} = T + d_{\mathcal{K}} \frac{D+1}{2}$ and if we select an odd degree $D$ we get $k_{\frac{D+1}{2}} = T$. The distance $d_{\mathcal{K}}= k_1-k_{0}$ is important for applications as it accounts for the degree of approximation of the time dependent effects and affects the amounts of parameter in the model, the smaller $d_{\mathcal{K}}$ the better the corresponding approximation, but the larger the considered parameter space. For $D$ we choose $D=3$, which provides the popular cubic splines which are twice continuously differentiable.
For mirroring the movement of the electricity price data we want that $\mu^i(t)$ smoothly changes over the days and during the week as well. So, for example, the coefficient for Sunday at 8 am is different to the coefficients at Sunday 9 am or Monday 8 am, but it is equal for every Sunday at 8 am.
However, we also observe that the mean behavior of at least some days is very similar to that of other days, see Figure (ref).
This finding is embedded in our model by proposing the assumption of some days having equal hourly mean components, which in turn reduces the number of parameters in the model. For example, the mean at Tuesday 8 am is assumed to be the same as on a Wednesday 8 am. Finally we consider five groups of parameter restrictions: Saturday, Sunday, Monday from 0 am to 12 am, Friday from 12 am to 12 pm and the core week from Monday 12 am to Friday 12 am. The Monday morning and Friday afternoon and evening has to be considered to get an appropriate phase-in and phase-out of the weekend. Friday, represented as the cyan colored line, departs from the typical structure of the other days (excluding the weekend) approximately at 12 am. It is almost the same for Monday, with the difference that Monday needs 7 to 8 hours to reach an equivalent structure comparable to the other weekdays. In order to maintain equidistant phase-in and phase-out periods, we set both periods to 12 hours.
Furthermore, there are significant impacts of public holidays, as mentioned by several authors. We consider every German public holiday as a Sunday. In addition to that there are some regional public holidays in Germany which are celebrated only by some areas. Every regional public holiday that concerns at least 25% of the population and is not on a Sunday, is treated as a Saturday as well. Moreover, we consider Christmas Eve and New Years Eve as regional public holidays. Consequently we take the 12 hours after a (regional) public holiday as Monday-am and the 12 hours before every (regional) public holiday as Friday-pm. Of course, this applies only if the day itself is not a Sunday, Saturday or (regional) public holiday. An overview of the explained day-grouping is given in Table (ref). In addition, the Table provides the number of unique hours within each group. Unique refers to an hour having one unique parameter. For instance, in case of the group normal, there are only 24 unique hours, even though its contained days, from Monday lunch time to Friday lunch time, has 96 hours in total. This approach accounts for the mean of hours of different days within this group being very similar.
To employ our model we denote the corresponding five groups of Table (ref) as ${\mathbb G} = \{1, \ldots, 5\}$ and create a time set ${\mathbb T}_{g}$. Every group $g$ is an element of ${\mathbb G}$, namely $g\in {\mathbb G}$. The combined time set ${\mathbb T}_{g}$ for all groups contains the integer value of the first hour of every group and day. The subset for the first group ${\mathbb T}_{1}$, for instance, contains only those integer values of hours, which lay inside the group full off and represent the hour 0 am of this day. For the group phase in and normal however, the subsets ${\mathbb T}_{3}$ and ${\mathbb T}_{5}$ contain the values of the hour 12 am for each day. The differences of the elements of ${\mathbb T}_{g}$ are therefore usually 12,24 or 168, but can also be a different number, when a public holiday occurs. An example for such a set ${\mathbb T}_{g}$ could therefore be ${\mathbb T}_{1}=\lbrace\ldots,0,168,336,504,528,672,\ldots\rbrace$, where at hour 528 (Monday) a holiday occurred. Within each group $g$ we set a specific number of basis functions. This number is related to the unique hours every group has. We determine the vector $N_{\mathbb G}$ as $N_{\mathbb G} = (d_{\mathcal{K}}^{\text{weekly}})^{-1} (24,24,12,12,24)'$. Please notice that we divide the amount of unique hours by $d_{\mathcal{K}}^{\text{weekly}}=4$, so that for instance the group full off contains $24/4=6$ different basis functions. In order to label these basis functions distinctively, we introduce the vector $C_{\mathbb G} = (d_{\mathcal{K}}^{\text{weekly}})^{-1} (0,24,48,60,72,96-d_{\mathcal{K}}^{\text{weekly}})' $, which contains the cumulative sum of the hours in $N_{\mathbb G}$ as entries. The subtraction of $d_{\mathcal{K}}^{\text{weekly}}$ in the last element of $C_{\mathbb G}$ is to avoid the singularity issue that occurs as the sum of all basis functions is constant. Then we define the function $G(\iota,C_{\mathbb G}) = g$ for $ C_{\mathbb G}(g) < \iota \leq C_{\mathbb G}(g+1)$. The index variable $\iota$ is integer valued, starts from 1 and contains in our case every natural number up to $(96-4)/4=23$. Hence, the function $G(\iota,C_{\mathbb G})$ matches every $\iota$, which represents the label of a basis function, to its specific group.
Moreover, we have to introduce the set ${\mathcal{C}} = \{C_{\mathbb G}(g)+1 | g\in {\mathbb G} \}$ which will contain the index of the first basis function of each group. In our case this turns out to be ${\mathcal{C}} = \{1,7,13,16,19\}$. The resulting basis functions $\widetilde{B}^{\mu^i}_j(t)$ can then be constructed as follows:
for $j\in {\mathcal{C}}$. The other basis functions can be obtained by shifting these basis functions. For $j \in \{ C_{\mathbb G}(g) +2, \ldots, C_{\mathbb G}(g+1)\}$ they are iteratively defined as $\widetilde{B}^{\mu^i}_j(t) = \widetilde{B}^{\mu^i}_{j-1}(t-d_{\mathcal{K}}^{\text{weekly}}) $. These basis functions apply for the price and the load data, so $i\in \{{\mathcal{P}}, {\mathcal{L}}\}$. Figure (ref) displays a practical application of equation ((ref)) to the time from 22.09.2012 to 05.10.2012. It shows the first basis function of each of the five groups. The model automatically detects that the 03.10.2012 is a public holiday and therefore breaks the usual periodicity for some of the groups, whereas the other groups remain stable. The group normal for instance leaves out the basis functions for Tuesday and Wednesday, even though they would usually be considered within this group. However, the group full off does now include the Wednesday as public holiday and therefore gets a basis function for this day.
Furthermore, we can observe an annual pattern within the data, most distinctively for the load and solar data. Hence, we added B-splines with a yearly periodicity of $S_{\text{annual}}=24\times365.24 = 8765.76$. As we assume that the meteorological processes have a smoother change over the year we added only $\eta^{\text{annual}}_{{\mathcal{R}}} = 6$ basis functions, whereas we take $\eta^{\text{annual}}_{i} = 12$ basis functions for the price and load, so $i\in \{{\mathcal{P}}, {\mathcal{L}} \}$. This leads to $d_{{\mathcal{K}},{\mathcal{R}}}^{\text{annual}} = 1460.96$ and $d_{{\mathcal{K}},i}^{\text{annual}} = 730.48$. For the wind and solar components we assumed no weekly structure, only a daily and a yearly one. Therefore we construct basis function for the daily structure with a periodicity of $S_{\text{daily}} = 24$ and annual basis functions. We choose $d^{\text{daily}}_{{\mathcal{K}}}=4$, thus there are $\eta^{\text{daily}} = 6$ basis functions added. We also take into account interactions between the yearly and daily components, i.e. how the daily cycle behavior changes over the year. This seems important as, for instance, the length of the sunshine period of the day changes over the year and should have a strong influence. This is modeled by a multiplication of both components.
In detail the additional basis functions for the price and load are given by $$ \widetilde{B}^{\mu^i}_{j}(t) = \sum_{k \in {\mathbb Z} } B_{k S_{\text{annual}} , d_{{\mathcal{K}},i}^{\text{annual}}}(t) $$ for $j = (d_{\mathcal{K}}^{\text{weekly}})^{-1} 96$ and $i\in \{{\mathcal{P}}, {\mathcal{L}}\}$. The other basis functions can be obtained by shifting these basis functions. For $j \in \{ (d_{\mathcal{K}}^{\text{weekly}})^{-1} 96 +1, \ldots, (d_{\mathcal{K}}^{\text{weekly}})^{-1} 96 +\eta_i^{\text{annual}} -1 \}$ they are iteratively defined as $\widetilde{B}^{\mu^i}_j(t) = \widetilde{B}^{\mu^i}_{j-1}(t-d_{{\mathcal{K}}, i}^{\text{annual}}) $. Hence the number of basis functions $M(\mu^i) = (d_{\mathcal{K}}^{\text{weekly}})^{-1} 96 + \eta^{\text{annual}}_{i} -1$ for $i\in \{{\mathcal{P}}, {\mathcal{L}}\}$.
For the wind and solar component we again define the basis functions for the daily and annual cycle as:
$$ \widetilde{B}^{*,\mu^{\mathcal{R}}}_{\text{daily},2}(t) = \sum_{k \in {\mathbb Z} } B_{k S_{\text{daily}}, d_{\mathcal{K}}^{\text{daily}}}(t) \text { and } \widetilde{B}^{*,\mu^{\mathcal{R}}}_{\text{annual},2}(t) = \sum_{k \in {\mathbb Z} } B_{k S_{\text{annual}} , d_{\mathcal{K}}^{\text{annual}}}(t). $$ Please remark that we set $\widetilde{B}^{*,\mu^{\mathcal{R}}}_{\text{daily},1}(t) = \widetilde{B}^{*,\mu^{\mathcal{R}}}_{\text{annual},1}(t) = 1$ which will become necessary in a later step, where we model the interactions. The subsequent B-splines are again iteratively defined with $\widetilde{B}^{*,\mu^{\mathcal{R}}}_{\text{daily},h_1}(t) = \widetilde{B}^{*,\mu^{\mathcal{R}}}_{\text{daily},h_1}(t - d_{\mathcal{K}}^{\text{daily}})$ for $h_1 \in \{ 3, \ldots, \eta^{\text{daily}} \}$ and $\widetilde{B}^{*,\mu^{\mathcal{R}}}_{\text{annual},h_1}(t) = \widetilde{B}^{*,\mu^{\mathcal{R}}}_{\text{annual},h_1}(t - d_{\mathcal{K}}^{\text{annual}})$ for $h_2 \in \{ 3, \ldots, \eta^{\text{annual}}_{\mathcal{R}} \}$. For modeling the interactions of the daily and annual cycle we can now introduce a modified version of basis functions for the renewables: $$ \widetilde{B}^{*,\mu^{\mathcal{R}}}_{j}(t) = \widetilde{B}^{*,\mu^{\mathcal{R}}}_{\text{daily},h_1}(t) \widetilde{B}^{*,\mu^{\mathcal{R}}}_{\text{annual},h_2}(t),$$
where $h_1 = j \operatorname{mod} \eta^{\text{daily}} +1 \text{ and } h_2 = j \operatorname{div} \eta^{\text{annual}}_{\mathcal{R}} + 1 $ for $j\in \{1, \ldots, \eta^{\text{daily}} \eta^{\text{annual}}_{{\mathcal{R}}} -1\}$. The term mod stands for the modulo operator and div is the integer without the remainder of a division. In this situation our previous definition of the first daily and annual basis functions is used, as the total set of B-spline now contains not only the simple daily and annual basis functions but also the interactions of them. Hence, we have included an amount of $M(\mu^{\mathcal{R}}) = \eta^{\text{daily}} \eta^{\text{annual}}_{{\mathcal{R}}} -1$ basis functions.
Furthermore, in Equation (ref), the parameters $\phi^{i,j}_k(t)$ might also depend on time. So the autoregressive behavior can also change over time, especially in a seasonal manner. A related approach was also used by koopman2007periodic and bosco2007deregulated for modeling daily electricity prices, and by guthrie2007electricity for intra-day prices. We assume that this dependence has a structure as in (ref), so for $\phi^{i,j}_k(t)$ it is given by $$ \phi^{i,j}_k(t) = \phi^{i,j}_{k,0} + \sum_{j=1}^{M_i(\phi_k^{i,j})} \phi_{k,j} \widetilde{B}_j^{\phi^{i,j}_k}(t).$$ Again we have to choose the basis functions $\widetilde{B}^{\phi^{i,j}_k}_j(t)$, but from an application point of view it is more practical to take the same basis functions as for $\mu^i$, which is implemented in our approach.
Unfortunately it would expand the parameter space enormously if we assume that for every possible lag $k\in I_{i,j}$ the coefficient $\phi^{i,j}_k(t)$ changes over time. Therefore, we let only those coefficients $L_{i,j}\subseteq I_{i,j}$ to change over time, which we assume to be most important. These are summarized in Table (ref), the other ones are empty.
Furthermore, there is a special impact on $Y_t$ that has not been yet modeled. Due to daylight saving time, we have a time shift in March and October. As we removed the added hour in October and extrapolated the missing hour in March, there is no problem in modeling the electricity spot price as the population adjusts their every day behavior very quickly. The same holds true for the electricity load if we assume that the influence of regenerative power feed-in is negligible. But for the renewable energy feed-in this looks different, here we can improve our model. In winter the sun peak is at 12 am, in summer it is at 1 pm. This has an impact especially on the solar power. Figure (ref) illustrates this effect on the observed solar power feed-in in the morning.
It can be observed that right after the time shift in March there is a downwards jump in the solar power feed-in present, and after the time shift in October a jump upwards occurs. We will model this behavior by a shift in the considered basis functions. Therefore, let ${\mathbb T}_{\text{DST}}$ denote the set of all time points with summer time. Then the considered time adjusted basis functions are given by $$\widetilde{B}^{\mu^{\mathcal{R}}}_{j}(t) = \boldsymbol 1_{\{t\notin {\mathbb T}_{\text{DST}} \}} \widetilde{B}^{*,\mu^{\mathcal{R}}}_{j}(t) + \boldsymbol 1_{\{t\in {\mathbb T}_{\text{DST}} \}} \widetilde{B}^{*,\mu^{\mathcal{R}}}_{j}(t+1).$$ This shift is applied to all basis function of the wind and solar power feed-in. For the mean component $\mu^{{\mathcal{R}}}(t)$ the effect is clear, here this is equivalent to a shift in the hour labels. In contrast, the effect on the autoregressive parameters is not that obvious. The feed-in $Y_{t}^{\mathcal{R}}$ depends e.g. on $\phi^{{\mathcal{R}},{\mathcal{R}}}_{24}(t)$ before the time shift and is different to $\phi^{{\mathcal{R}},{\mathcal{R}}}_{24}(t)$ after the time shift. But $\phi^{{\mathcal{R}},{\mathcal{R}}}_{24}(t)$ would be the same as $\phi^{{\mathcal{R}},{\mathcal{R}}}_{24}(t-25)$ in March or $\phi^{{\mathcal{R}},{\mathcal{R}}}_{24}(t-23)$ in October, due to the time shift in daylight savings time. Nevertheless, in our approach this is not the case, as we explicitly model the daily and yearly interaction effects. The coefficients therefore truly embed the effects of the time shift.
But also the variance structure of the model is relevant. Here we assume that $\Sigma_t = \text{diag}(\boldsymbol \sigma_t) =\text{diag}(\sigma^{\mathcal{P}}_t,\sigma^{\mathcal{L}}_t,\sigma^{\mathcal{R}}_t) $. The error term $\boldsymbol \varepsilon_t$ can be written as $${\varepsilon}^i_t = \sigma^i_t Z^i_t \text{ where } (Z^i_t)_{t\in {\mathbb Z}} \text{ is i.i.d. with } {\mathbb E}(Z^i_t)=0 \text{ and } \operatorname{{\mathbb V} ar}(Z^i_t) =1 .$$ As modeled by keles2012comparison we want to explain conditional heteroscedasticity by an ARCH and GARCH type model. This approach is suitable to model the price spikes, especially if the underlying distribution of $Z_t$ is heavy tailed.
But as the price process $\boldsymbol Y_t$ is heavy tailed due to spikes, the residual process $\boldsymbol \varepsilon_t$ is likely to be heavy tailed with a small tail index, too. We noticed that it is not suitable to consider the square of ${\varepsilon}^i_t$, as done in a typical GARCH process, because they become even more heavy tailed. This halves the corresponding tail index and the convergence of the estimators will get worse. Instead of taking squares, we use the absolute value $ |{\varepsilon}^i_t|$, which keeps the tail index on the same level. Some GARCH-type processes already use this structure, e.g. the TARCH model by rabemananjara1993threshold. Additionally, a TARCH and TGARCH model can handle the so called leverage effect, so that negative and positive past residuals have different influences on the volatility. We assume a TARCH process for $\sigma_t$ so that
holds with $\alpha^{+,i}_k,\alpha^{-,i}_k\geq 0$ where ${\varepsilon}^{+,i}_t = \max\{{\varepsilon}^{i}_t, 0\}$ and ${\varepsilon}^{-,i}_{t} = \max\{-{\varepsilon}^{i}_t, 0\}$. For the TARCH process we will assume that only the trend component $\alpha^i_0(t)$ varies seasonally over time. Thus the full parameter vector $\boldsymbol \alpha_i$ is given by $( \alpha^i_{0,0} , \ldots, \alpha^i_{0,M(\alpha^i_{0})} , (\alpha^{+,i}_k,\alpha^{-,i}_k)_{k\in J_{i}} )$. For $\alpha^i_0(t)$ we consider the same basis functions as for the corresponding $\mu^i(t)$, but without the linear trend, as it could conflict the positivity constraint to the parameters. Here we make use of another advantage of the periodic B-spline approach in contrast to a Fourier approximation. The positivity of the considered basis functions gives us automatically the restriction that the corresponding parameters have to be positive as well as $\sigma_t$ must be positive.
The index sets $J_i$ that we use are given in Table (ref), they contain the typical discussed lags.
By multiplying $\gamma^i = {\mathbb E}|Z^i_t|$ and adding $v^i_t = \sigma^i_t(|Z^i_t| - \gamma^i)$ on both sides of equation (ref) we get
where $\tilde{\alpha}^{i}_0(t) = \gamma^i \alpha^{i}_0(t)$, $\tilde{\alpha}^{+,i}_k = \gamma^i \alpha^{+,i}_k$, and $\tilde{\alpha}^{-,i}_k = \gamma^i \alpha^{-,i}_k$. Here $v^i_t$ is a weak white noise process with ${\mathbb E}(v^i_t) = 0$. The auxiliary model in (ref) is useful as it allows us to estimate $\boldsymbol \alpha_i$ given estimates of $\boldsymbol \varepsilon^i$. The fitted values $\tilde{\sigma}^i_t$ of equation (ref) are proportional to $\sigma^i_t$ up to the constant $\gamma^i$. $\gamma^i$ is the first absolute moment of ${\varepsilon}^i_t$ which might help to characterize the distribution of ${\varepsilon}^i_t$. If ${\varepsilon}^i_t$ follows a normal distribution it is exactly $\sqrt{2\pi^{-1}} \approx 0.798$, whereas, e.g., the standardized t-distributions have a larger first absolute moment.
Given reasonable estimates $\widehat{{\varepsilon}}^i_t$ and $\widehat{\tilde{\sigma}}^i_t$ for ${\varepsilon}^i_t$ and $\tilde{\sigma}^i_t$ we can use the plug-in estimator for $\gamma^i$ to estimate it. With $\kappa^i_t = (\widehat{\tilde{\sigma}}^i_t)^{-1} \widehat{{\varepsilon}}^i_t$ this is given by
as $\operatorname{{\mathbb V} ar}(\kappa_t^i)=({\gamma^i})^{-2}$. $\sigma^i_t$ can now be estimated with $\widehat{\sigma}^i_t = (\widehat{\gamma}^i)^{-1} \widehat{\tilde{\sigma}}^i_t $.
Finally, we can rewrite (ref) in matrix notation with
with $\boldsymbol Y^i = (Y_1^i, \ldots, Y^i_n)$, $\boldsymbol X^i = (\boldsymbol X^{i}_1, \ldots, \boldsymbol X^i_{n})$, $\boldsymbol \varepsilon^i = ({\varepsilon}^i_1, \ldots, {\varepsilon}^i_n)$ and parameter vector $\boldsymbol \theta_i$ for $i\in {\mathcal{I}}$ given time points $t\in \{1,\ldots, n\}$. The vectors $\boldsymbol X_t^i$ are given through equation ((ref)) by $$\boldsymbol X_t^i = \left( 1, t, \boldsymbol B_t(\mu^i), \boldsymbol B^{Y^i}_t(I_{i,{\mathcal{P}}}) , \boldsymbol B^{Y^i}_t(I_{i,{\mathcal{L}}}) ,\boldsymbol B^{Y^i}_t(I_{i,{\mathcal{R}}}) , \boldsymbol B^{Y^i, \text{per}}_{t}(L_{i,i}) \right)$$ where we have $\boldsymbol B_t(\mu^i) = (\widetilde{B}^{\mu^i}_1(t), \ldots, \widetilde{B}^{\mu^i}_{M(\mu^i)}(t) )$, $\boldsymbol B^{Y^i}_t(I_{i,j}) = (Y^i_{t-k})_{k\in I_{i,j}}$, for $j \in {\mathcal{I}}$, and $\boldsymbol B^{Y^i,\text{per}}_t(L_{i,i}) = (Y^i_{t-k} \boldsymbol B_t(\mu^i) )_{k\in L_{i,i}} $. In the following we will refer the parameters that correspond to $(1, t, \boldsymbol B_t(\mu^i))$ as deterministic component, $\boldsymbol B^{Y^i}_t(I_{i,{\mathcal{P}}})$ as price component, $\boldsymbol B^{Y^i}_t(I_{i,{\mathcal{L}}})$ as load component, $\boldsymbol B^{Y^i}_t(I_{i,{\mathcal{R}}})$ as wind and solar component and $\boldsymbol B^{Y^i, \text{per}}_{t}(L_{i,i})$ as periodic component given an $i\in {\mathcal{I}}$.
All in all our model has approximately 3500 possible parameters, which are given in Table (ref). They match the mentioned deterministic, price, load, wind and solar, and periodic component defined above. Please notice that these parameters only represent the maximum amount of possible parameters, as our estimation procedure, which will be explained in the next section, will eliminate parameters, which are not significantly different from zero.
For the estimation of the parameters we will consider a variation of the iteratively weighted least squares approach as used for example in mbamalu1993load or mak1997estimation. So basically we apply a weighted least squares estimation of model (ref). Then we use the corresponding residuals to estimate the volatility which is used to compute new weights for re-estimating model (ref).
The considered model has a lot of regressors, and it is likely that some of them have no significant impact on the model, but rather increase spurious effects. The suggested lasso (least absolute shrinkage and selection operator) approach can handle this in an efficient way. It was introduced by tibshirani1996regression in the context of shrinkage and selection in regression models and was recently applied by hsu2008subset and ren2010subset in the context of vector autoregressive models.
The weighted lasso estimator of the parameter vector $\boldsymbol \theta^i$ of equation (ref) or (ref) is given by
with a weight vector $\boldsymbol w^i = (w^i_1, \ldots, w^i_n)$, $p_i$ as length of $\boldsymbol \theta^i$ and tuning parameter $\lambda_{i,n}$. In the initial step $\boldsymbol w^i$ is chosen to be $\boldsymbol 1$. For the estimation we use the least angle regression (LARS) estimation algorithm of efron2004least.
As selection criterion we consider the Akaike information criterion (AIC) to estimate . This is asymptotically equivalent to the cross-validation in regression analysis as suggested in efron2004least
In the next step we want to compute a new weight matrix $\boldsymbol W = (\boldsymbol w^1, \ldots, \boldsymbol w^d)$ based on the volatility estimates for $\boldsymbol \sigma^i_t$. Thus we have to estimate the volatility parameters from equation (ref) using equation (ref) and (ref).
Here we have to ensure that all parameters take non-negative values, as the regressors are all positive. To solve this non-negative least squares problem, we use the NNLS algorithm as described by lawson1974solving. Note that this estimation technique can be seen as a parameter selection approach as well, as some parameters are likely estimated to be 0. This sparsity effect in a NNLS problem was recently analyzed by meinshausen2013sign. slawski2013non showed that the non-negative least squares approach is potentially superior to the positive lasso, a lasso approach with a positive parameter constraint.
The general estimation scheme is given by
\framebox{ \parbox{.89\textwidth}{
} }
Note that before performing the LARS algorithm on the weighted regressor matrix, we standardize their columns so that the columns of $\boldsymbol X^i$ from equation (ref) have zero sample mean and sample variance $1$, so that LASSO can act as a proper model selection algorithm. Additionally we could also replace $\widehat{\boldsymbol \sigma}_t$ by $\widehat{\tilde{\boldsymbol \sigma}}_t$ as they are the same upto the constant $\widehat{\gamma}_i$. This does not change the result of the estimation algorithm, but we can improve the computation time slighty as we do not have to estimate $\gamma^{i}$ by equation (ref).
For the stopping criteria we suggest to look at the convergence of the $\boldsymbol \sigma^i$ within the algorithm. A reasonable criteria is to iterate until $\delta_K^i = \| \widehat{\boldsymbol \sigma}^i_{K-1} - \widehat{\boldsymbol \sigma}^i_{K}\| <\epsilon$ for all $i\in {\mathcal{I}}$. For $K=1$ the $\boldsymbol \sigma^i_{1}$ is the homoscedastic estimate from the homoscedastic lasso estimation. For our application we choose $\delta_K^i$ as $n^{-1} \|\widehat{\boldsymbol \sigma}^i_{K-1} - \widehat{\boldsymbol \sigma}^i_{K}\|_1$ and $\epsilon=0.001$. In the empirical study we noticed that this algorithm converges quickly (see Figure (ref)), so a small number of iterations $K_{\text{max}}$ seems to be sufficient for practitioners. In Figure (ref) we can observe that in the fifth iteration the change in $\widehat{\boldsymbol \sigma}^i$ is less than $\epsilon$. Hence, we decide to stop after $K_{\text{max}} = 4$ iterations in the rest of our study.
The estimation technique used has the huge advantage of being fast, even if hundreds of potential regressors are included in the model. The computational complexity of the dominating LARS algorithm is the same as for a common OLS estimation, which is ${\mathcal{O}}(K_{\text{max}} dp^2n)$ for $p<n$ where $K_{\text{max}}$ is the number of maximal iteration steps, $d$ the dimension of the process, $p$ the maximal number of regressors in (ref), and $n$ the number of observations. Here we see that the number of parameters $p$ in the models has a quadratic impact, thus it is important to keep the parameter space manageable regarding the computation time. In our computations the iteratively reweighted lasso with $K_{\text{max}}=4$ takes only a few minutes on a 3 GHz computer, with $n$ about 31,000, $p$ about 3500, and $d=3$.
The consistency and asymptotic behavior of $\widehat{\boldsymbol \theta}_i$ and $\widehat{\boldsymbol \alpha}_i$ for the process considered is not obvious. However, the lasso-type estimators in a heteroscedastic regression setting were recently analyzed in severien2012shrinkage, wagener2012bridge, and wagener2013adaptive.
There are also some asymptotic results for the lasso estimator in an autoregressive setting. Important results are given for univariate ARX processes by wang2007regression and for the VAR model by hsu2008subset.
Even though there is a lack of theoretical justification for the asymptotic properties of the used estimators, we will assume asymptotic normality for $\widehat{\theta}_i$ and $\widehat{\boldsymbol \alpha}_i$. For $\widehat{\theta}_i$ we assume the asymptotic behavior as in wagener2012bridge. Following their result we have $\sqrt{n} (\widehat{\boldsymbol \theta}_i(1) -\boldsymbol \theta_i^*(1)) \to N(\boldsymbol \zeta_i(\boldsymbol \theta_i^*), \boldsymbol \Gamma_i(\boldsymbol \theta_i^*) )$ in distribution where $\boldsymbol \theta_i^*$ is the true parameter vector, $\widehat{\boldsymbol \theta_i}(1)$ and $\boldsymbol \theta_i^*(1)$ the corresponding non-zero parts, $\boldsymbol \zeta_i(\boldsymbol \theta_i^*)$ as asymptotic mean vector, and $\boldsymbol \Gamma_i(\boldsymbol \theta_i^*)$ the asymptotic covariance matrix, if $n^{-\frac{1}{2}} \lambda_{i,n} \to \lambda_{i,0} \geq 0$. For the mean and the covariance matrix we have $\boldsymbol \zeta_i(\boldsymbol \theta_i^*) = - \frac{\lambda_{i,0}}{2} \boldsymbol G_i^{-1} \text{sign}( \boldsymbol \theta_i^*(1) ) $ and $\boldsymbol \Gamma_i(\boldsymbol \theta_i^*) = \boldsymbol G_i^{-1} \boldsymbol G_i^{\boldsymbol W} \boldsymbol G_i^{-1}$ with $\boldsymbol G_i$ as limit of the Gramian $\boldsymbol G_{n,i} = n^{-1}\boldsymbol X_i^\top \boldsymbol X_i$, $\boldsymbol G_i^{\boldsymbol W}$ as limit of the weighted Gramian $\boldsymbol G^{\boldsymbol W}_{n,i} = n^{-1}\boldsymbol X_i^\top \boldsymbol W_i \boldsymbol X_i$, weight matrix $\boldsymbol W_i = \text{diag}({(\sigma_1^i)}^{-2},\ldots, {(\sigma_n^i)}^{-2})$, and $\text{sign}$ as signum function.
We see that the asymptotic behavior depends linear on the unknown limit $\lambda_{i,0}$. For analyzing the asymptotic in our situation we estimated the model for many different $n$ ranging from $2000$ to our maximum sample size and computed the corresponding $\lambda_{n,i}$. In Figure (ref) we plotted $n$ against $n^{-\frac{1}{2}}\lambda_{i,n}$.
There we can observe that $n^{-\frac{1}{2}}\lambda_{i,n}$ goes clearly to zero as $n$ increases. Thus the assumption $\lambda_0=0$ seems to be appropriate, and the zero vector $\boldsymbol 0$ is a reasonable choice as an estimator for the mean vector $\boldsymbol \zeta_i(\boldsymbol \theta_i^*)$. In addition the covariance matrix can be estimated by the common plug-in estimator $ \widehat{\boldsymbol \Gamma}_i = n( \boldsymbol X_i^\top \boldsymbol X_i )^{-1} \boldsymbol X_i^\top \widehat{\boldsymbol W}_i \boldsymbol X_i ( \boldsymbol X_i^\top \boldsymbol X_i )^{-1}$.
We performed the mentioned estimation technique on the full data set. For the model diagnostic the sample autocorrelation function (ACF) of the estimated standardized residuals and their absolute values are given in Figure (ref).
The left 3x3 matrix depicts how the standardized residuals of one time series correlate with the lagged standardized residuals of the same or any other of the three time series. The right 3x3 matrix shows the same relation for the absolute standardized residuals. For instance, the left upper ACF shows the correlation of the standardized residuals of the price with its own lagged standardized residuals. The illustration right next to it, Price & Load, shows the correlation of the standardized residuals of the price with the lagged standardized residuals of the load. Hence, we can see that the independence assumption for $\boldsymbol Z_t$ seems to be satisfied for most of the relations. However, the strongest serial correlation structure can be observed for the standardized residuals $Z_t^{{\mathcal{L}}}$ of the load. There could be several reasons for the remaining correlation structure. For example, we may simply did not choose enough parameters to cover the complex dependence structure, so an enlargement of $I_{{\mathcal{L}},{\mathcal{L}}}$ or $L_{{\mathcal{L}},{\mathcal{L}}}$ or a reduction of $d_{\mathcal{K}}$ might help. It could be also possible that there might be non-linear or interactions effects that were not considered so far.
Further, the sample autocorrelation of $|\boldsymbol Z_t|$ is mainly zero, so that $\boldsymbol Z_t$ seems to be quite homoscedastic. Thus the TARCH approach seems to be appropriate to explain the conditional heteroscedasticity. Nevertheless, there is still space for improvements.
Given the estimated model (ref) we can evaluate the t-values of every coefficient given by their estimate over its standard error which can be estimated under the asymptotic normality assumption. Roughly speaking we can say, that the larger the absolute t-value of a coefficient, the more significant is its impact on the process. The evaluated t-values are given in Figure (ref), distinguished into the groups introduced at the end of section (ref).
It is worth mentioning that the load and wind and solar feed-in seem to have a strong impact on the actual electricity price, $Y_t^{\mathcal{P}}$. The most recent past load values (lag 1 and 2) seem to be very important as well as the values one week ago (lag 168, 169). In contrast to this, wind and solar power has a vast impact only in the short-run at lag 1. For the load dependency structure, $Y_t^{\mathcal{L}}$, we receive that the first autoregressive lag and the subsequent daily lags tend to be of great importance. Interestingly, lag 337, which accounts for effects two weeks ago, seem to also have a special importance. But there is also a strong relationship concerning the renewable energy feed-in, which is especially high for lag 1 and 2. For the solar and wind power feed-in, $Y_t^{\mathcal{R}}$, it is interesting that the first two autoregressive lags clearly dominate all the other larger lags. However, the daily lagged variables still play an important role.
In general we can observe that the dependency structure is quite complex. Therefore we conclude that the chosen set of lags $I_{i,j}$ and $L_{i,i}$, though having thousands of parameters, seem to be reasonable. This is also backed by the results of the parameter selection algorithm. For the price process $82.2\%$ of the chosen parameters were included. The load process was set to $92.8\%$ and wind and solar to $76.6\%$ of the maximum amount of parameters. Especially the high percentage for the load indicates that an enlargement of the parameter set is likely to improve the results.
Assuming asymptotic normality allows us to compute several confidence intervals and bands. For example, for $\boldsymbol \mu(t)$ and time varying autoregressive parameters $\phi^{i,j}_k(t)$. We illustrate the behavior of $\mu^{\mathcal{P}}(t)$ and $\phi^{{\mathcal{P}},{\mathcal{P}}}_1(t)$ in Figure (ref). From this illustration it can be seen that all coefficients have distinct weekly patterns and a seasonal structure.
Moreover, analyzing the volatility structure determined by $\boldsymbol \sigma_t$ yields promising results. Using the TARCH model for $\boldsymbol \varepsilon_t$ allows us to decompose $\boldsymbol \sigma_t$ into three components. The first one is determined by the periodic coefficient $\boldsymbol \alpha_0(t)$. It describes the deterministic part of the volatility and is automatically a lower bound for $\boldsymbol \sigma_t$. The second and third components describe the residual impact of the positive and negative residuals, respectively. Note that in an ARCH model, both components are equal. The decomposition of the estimated volatilities $\widehat{\sigma}_t^i$ is given in Figure (ref).
For illustrative purposes, we selected approximately two weeks of October/November 2011. It can be seen that the volatility of each of the three time series is influenced by a deterministic seasonal component. Moreover, for the price process, we cannot observe a strong leverage effect, the impact of past negative residuals is almost equal to the impact of the past positive ones. This is represented by the area marked as negative in the left picture being approximately equal to the area marked as positive. This relation is obviously different for load time series. Here we observe a strong leverage effect. The impact of negative residuals is higher than the impact of positive residuals. For the wind and solar feed-in an inverse leverage effect is present, as the positive area clearly dominates the negative area.
Furthermore, we can perform tests to check the significance of the leverage effect. This will be especially interesting for the price series, as judging from our subset of October/November 2011 it was not clear whether a significant effect occurred. Therefore we evaluate several tests $H_0: A_{i,k} = \sum_{j=1}^k \alpha_j^{+,i} - \alpha_j^{-,i} = 0$ for $i \in {\mathcal{I}}$. The $k$ indicates the lags upto which the leverage effect is taken into account.
Figure (ref) provides the t-values of the corresponding test for all possible $k$. In order to construct the different values for $k$ we simply cut the remaining model parameters of the TARCH-part after $k$ and performed the test described above. We therefore explicitly state that we did not re-estimate the model for every $k$. It is obvious that the wind and solar energy has clearly positive t-values with its peak at about a day. Those positive t-values indicate a strong inverse leverage effect. For the load time series the leverage effect seems to point in the opposite direction. Including all model parameters we can state, that there is a significant standard leverage effect within the data. This is also supported by the t-value for the load in Table (ref). However, the picture becomes fuzzy when the electricity price is considered. Judging by the Figure it can be obtained, that the electricity prices firstly experiences a positive leverage effect, which diminishes after a certain amount of time. When all parameters up to $k$ are considered, the leverage effect even becomes insignificant, as can be seen in Table (ref). However, applying a variance model which can account for a leverage effect is still important also for electricity prices, as only those models can account for such a complex structure with different signed impacts of different lags on the volatility.
In the same way we can perform an approximate test as to whether an increasing load or wind-solar feed-in leads to a significantly increasing or decreasing price in the long run. These tests have the null hypothesis $H_0 : \Phi_{{\mathcal{P}},{\mathcal{L}}} = \sum_{k\in I_{{\mathcal{P}},{\mathcal{L}}}} \phi^{{\mathcal{P}}, {\mathcal{L}}}_{k} = 0$ and $H_0 : \Phi_{{\mathcal{P}},{\mathcal{R}}} = \sum_{k\in I_{{\mathcal{P}},{\mathcal{R}}}} \phi^{{\mathcal{P}}, {\mathcal{R}}}_{k} = 0$. Applying this test to the data (see Table (ref)) it can be concluded that there is only a small long run impact of the load to the electricity price, which is statistically not significant. In contrast, an increasing in the amount of wind or solar power feed-in causes a decrease in the price. Taking the standardization of $Y_t^i$ by its mean and standard deviation before the estimation of (ref) into account we can quantify these results. An increase of 1 GWh in the load changes the price by $-0.108 (\pm 0.499)\footnote{\text{The the term in parentheses gives the symmetric $90%$ confidence interval}} \frac{\text{EUR}}{\text{MWh}}$ and an increase of 1 GWh wind or solar power feed-in changes the price by $-2.031 ( \pm 0.375) \frac{\text{EUR}}{\text{MWh}}$. This result concerning the wind and solar feed-in is consistent with the recent literature on quantifying the impact of renewable energy in Germany, see e.g. wurzburg2013renewable.
Moreover, we estimated a linear trend in $Y_t^{\mathcal{P}}$ which is negative and significantly different from zero, so over time, the price seems to decrease. But this effect is not that strong: over one year, the price changes by $-0.0200 (\pm 0.0069)\frac{\text{EUR}}{\text{MWh}}$. Nevertheless, we have to be careful with extrapolating a linear trend far in the future, as it is economically impossible for such a negative trend to last forever.
Given the estimated model (ref) of the sample $(\boldsymbol Y_1, \ldots, \boldsymbol Y_n)$ we can easily carry out a forecast. Let $\boldsymbol Y_t = g(\boldsymbol Y_{t-1}, \boldsymbol Y_{t-2}, \ldots)$ be the representation of (ref). Then we can compute $\widehat{\boldsymbol Y}_{n+h}$ iteratively by $$\widehat{\boldsymbol Y}_{n+h} = g(\widehat{\boldsymbol Y}_{n-1+h}, \widehat{\boldsymbol Y}_{n-2+h}, \ldots) $$ and defining $\widehat{\boldsymbol Y}_{t}:= \boldsymbol Y_t$ for $t\leq n$.
We performed a forecasting study, where we choose subsequences of the $n=31465$ observations. So we choose a data part $H_l$ that is approximately two years long (exactly 18481 observations = 110 weeks + 1 hour), starting at observation $l$. The $l$ is chosen so that the observed sample ends with an observation for the electricity price from 23h-24h. Finally, we performed a complete estimation on the data sample $H_l$ and estimate the next $h=672$ hours, that is four weeks. This method provides a fair estimation technique, as we use the same amount of past observations for every forecast. Given our sample the procedure provides $N=506$ observations, while $\widehat{\boldsymbol Y}_{h,k}$ and $\boldsymbol Y_{h,k}$ for $k \in \{1,\ldots, N\}$ denote the corresponding predicted values and observation, respectively. Note that we are forecasting the three-dimensional process $\boldsymbol Y_t$. However, in this paper we are interested in the price process $Y_t^{\mathcal{P}}$, therefore the further discussion will focus on this process.
The first benchmark we consider is the homoscedastic solution of model (ref). We get its forecasts automatically by taking the estimated parameters after the first iteration step within the iterative estimation algorithm. Thus, we can directly see the impact of the considered TARCH part.
As other own benchmarks we consider the weekly persistent process $Y_t^{\mathcal{P}} = Y_{t-168}^{\mathcal{P}}$ as persistent model, such as two AIC selected (V)AR processes with time varying mean. Their model is given by $$\boldsymbol X_t = \boldsymbol \mu_j \boldsymbol 1_{\{t \in (168 {\mathbb N} + j) \}} + \sum_{k=1}^p ( \boldsymbol{\Phi}_k \boldsymbol X_{t-k}- \boldsymbol \mu_j \boldsymbol 1_{\{t \in (168 {\mathbb N} + j) \}}) + \boldsymbol \varepsilon_t$$ where $j \in \{1, \ldots, 168\}$. We consider the univariate choice $\boldsymbol X_t = Y_t^{\mathcal{P}}$ and the two-dimensional incorporating the load as well, so $\boldsymbol X_t = (Y_t^{\mathcal{P}}, Y_t^{{\mathcal{L}}})$. We also tried to include the wind and solar feed-in but the out-of-sample performance got worse. We estimated the models in a two step approach, first removing the weekly mean and second estimating a (V)AR process via Gaussian AIC selection, where the $168$ parameters for the mean are ignored. For the estimation process we solve the Yule--Walker equations that provide a guaranteed stationary solution. As maximal possible order for the univariate process, we choose $p_{\max}=1210$ and for the two-dimensional one $p_{\max}=555$. The estimation of the processes is very fast and done in a few seconds. However, in our empirical results, the order $p$ of the AR process is usually automatically chosen to be about 800, whereas the in the two-dimensional VAR case, it is usually about 400, so both cover the weekly mean and conditional mean behavior.
Moreover, we tried to use as many models from the recent literature for a benchmark as possible. But unfortunately finding the right competitor is difficult for several reasons. First of all many authors only consider positive observations or delete some chunks of the data (e.g. outliers, holidays, etc.). Second, some of them consider information from random regressors such as the load as known, which in a real world situation would not be the case. And finally there are some models, especially from machine learning, that simply require an enormous computational time. Furthermore, note that feed-forward neural networks with one hidden layer and lagged process $\boldsymbol Y_t$ as input acts very similar to an AR($p$) process if the chosen lags in the input layer are covered in an AR($p$) process. But, as mentioned, neural networks are extremely time consuming, due to the learning phase, which is often based on random selections.
However, we consider three benchmarks from the literature that have been applied to electricity prices. First, an ARMA(5,1) model with trend as well as annual, weekly and daily cycles as suggested as one of the best models in keles2012comparison. Note that their daily cycles vary over the seasons of the year. Moreover, we consider the functional data analysis approach from liebl2013modeling. But he also modifies the data: he removes outliers and excludes holidays and weekends. In fact, he models the prices by estimating the merit order curve under consideration of the load subtracted by the wind power feed-in. But for our data we noticed that his model gains better results when only the load is used. Furthermore, in his studies, two functional principal components were sufficient to model the data well. But we got better results by using three principal components. In the original paper the loading coefficients were predicted using a SARIMA model with periodicity five, whereas we use the same model with a periodicity seven as we are not ignoring the weekends. Moreover, we remark that the usage of the basic model of liebl2013modeling as a benchmark to our model is limited, as his approach was tailor-made for the exclusion of special days like weekends and heavy outliers. It may therefore be the case that our reported model performance for his model is only due to the violation of some of his basic assumptions. As a third benchmark from the existing literature we employ the wavelet-ARIMA approach from conejo2005day. Such a model is often used as a benchmark in electricity spot price forecasting. In our application we use the Daubechies 4 wavelet. For modeling the coefficients of the wavelet decomposition we choose ARIMA(12,1,1) processes to capture their autoregressive structure.
As a performance measure for our forecasting study, we do not use MAPE or WMAPE as is often done in the literature, because the $\text{MAPE}^{\mathcal{P}}_h = \frac{1}{N} \sum_{k=1}^N \big| \frac{ Y^{\mathcal{P}}_{h,k} - \widehat{Y}^{\mathcal{P}}_{h,k} } {Y^{\mathcal{P}}_{h,k}} \big| $ is obviously pointless for data that can take zero values. Instead, we consider the mean absolute forecast error for the prediction time $h$ ($\text{\textbf{MAE}}_h$) as well as the mean of the mean absolute forecast error up to a forecast horizon of $h$ ($\text{\textbf{MMAE}}_h$). They are given by $$ \text{\textbf{MAE}}_h = \frac{1}{N} \sum_{k=1}^N |\boldsymbol Y_{h,k} - \widehat{\boldsymbol Y}_{h,k}| \ \ \ \text{ and } \ \ \ \text{\textbf{MMAE}}_h = \frac{1}{hN} \sum_{j=1}^h \sum_{k=1}^N |\boldsymbol Y_{j,k} - \widehat{\boldsymbol Y}_{j,k}| .$$ As main criterion to compare different models we suggest the $\text{\textbf{MMAE}}_{24}$, as it estimates the expected average absolute error for the whole next trading day, which is especially for the practical application of relevance. The computed $\text{MAE}_h^{\mathcal{P}}$ and $\text{MMAE}_h^{\mathcal{P}}$ for the considered models with selected prediction horizons are given in Table (ref). The evolution of $\text{MAE}_h^{\mathcal{P}}$ and $\text{MMAE}_h^{\mathcal{P}}$ with increasing $h$ is presented in Figure (ref).
It is remarkable that our proposed model can outperform every other one, especially these that are used in the literature. Remember that we compare pure out-of-sample methods, that use no information from the prediction time during the estimation. Moreover, we noticed that many models are not able to outperform simple AIC selected AR($p$) or VAR($p$) models, so e.g. keles2012comparison and liebl2013modeling. Such benchmarks have not not been considered in most of the literature so far, even though they are simple and fast to estimate. The relatively good performance is likely due to the highly carefully chosen model order, so that hundreds of variables are included in the model, which can cover well the behavior. This amount of chosen parameters is in the same range as for our reweighted lasso procedure.
Furthermore, we use Monte Carlo methods to compute the prediction bands conditioned on the given data. Thus, we can perform a residual based bootstrap on the standardized residuals $\widehat{\boldsymbol Z}_t$ to simulate $\boldsymbol Y_{t}$. So we resample from $\widehat{\boldsymbol Z}_t$ to simulate $\Sigma_t$ and afterwards $\boldsymbol Y_t$ to compute the prediction intervals.
For illustrative purposes, we performed a prediction for the sample that is given in Figure (ref), which includes a public holiday. So the considered estimation period is about two years and the prediction horizon is 192 hours. The corresponding point estimates, such as the estimated conditional mean and the $90\%$ and $99\%$ prediction intervals are given in Figure (ref). The prediction bands are computed by evaluating the symmetric conditional value at risk (VaR) resp. quantile levels.
Interestingly, the forecasting performs quite well. For example, Wednesday the 3rd October is an official holiday, so the forecasted values are significantly smaller than those, for example, for the 4th October, even though in this case the forecast still overestimates the real value.
The paper presented a model for the hourly electricity price of the European Power Exchange for Germany and Austria. A periodic VAR-TARCH approach was proposed in order to capture the specific price movements. This model is able to deal with the most difficult challenges in modeling electricity prices. Namely, it considers: the mean-reversion of prices, the seasonality of the data, the possibility of vast price spikes, negative prices, periodicity, time dependent variance, the leverage effect, the impact of renewable energy, the impact of the electricity load, calender effects like holidays, and time shifts due to daylight saving time. For the necessary simultaneous modeling of the price, load, wind and solar power, we used a modern estimation approach, which is efficient and turned out to enable rapid estimations. By fitting our model to the before mentioned data sets we were able to show insides towards the leverage effect of our considered time series. Moreover, we provided evidence for the effect of renewable energy in reducing the price of electricity. Our study showed that an increase in combined wind and solar power of 1 GWh leads to a decrease in price of 2.03 $\frac{\text{EUR}}{\text{MWh}}$. An extensive forecasting study showed that in MAE and MMAE our model outperforms every model which was used as a benchmark, including recently published models in the literature. For the purpose of forecasting, no future knowledge was necessary: every forecast was made exclusively with out-of-sample data. Due to its high efficiency and fast computing time, our model may be of interest for researchers who want to use it as a benchmark for their own model as well as for energy companies who need to forecast electricity prices. Nevertheless, there is still much research left for the future. For instance, as the markets of some European countries are strongly integrated, measuring the import and export feeds may lead to better estimation results and provide guidance for political decisions. It is also debatable whether our model can be reasonably applied to other electricity markets. We notice that the characteristics of the energy portfolio and working day structure of Germany had significant impacts on the price process. But those characteristics are not necessarily possessed by other markets. Therefore, introducing a more comprehensive model may be an appropriate future task. Another direction for future work could also be a more detailed examination of the German energy portfolio, by considering the time series of prices for other energy resources, e.g. coal or gas. Detecting and using the cointegration between those time series may also lead to an improvement in the estimation results. From the statistical point of view the considered model can be extended as well. So it is possible to include non-linear effects and interactions or change points within the variables. Further one might allow more parameter to vary periodically, as, for instance, it might be that the TARCH parameters change over time.