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.
70,358 characters · 20 sections · 0 citation commands
Probabilistic Multi-Step-Ahead Short-Term Water Demand Forecasting with Lasso
\affil[1]{Faculty of Economics, esp. Economics of Renewable Energies, University Duisburg-Essen, Universitätsstr. 2, 45141 Essen, Germany, with corresponding author email. Email: [email removed]} \affil[2]{Faculty of Economics, esp. Economics of Renewable Energies, University Duisburg-Essen, Universitätsstr. 2, 45141 Essen, Germany.}
The increasing availability of high-frequency data in the water sector brings new opportunities to further refine and optimize the efficiency of water utilities. Here, the use of data-driven, short-term water demand forecasting models to reduce energy costs has become a popular field of application in practice, as noted by \citeN{Alvisi.2007}, \citeN{Brentan.2017b}, and \citeN{Arandia.2016}. Reliable information about the expected demand allows, for instance, for optimizing the control of storage capacities to balance demand peaks and to run and schedule pumps more efficiently.\\ The quantification of the underlying uncertainty of future demand forecasts turned out to be of great importance in the decision-making process as emphasized by \citeN{Donkor.2014} and \citeN{Hutton.2015}. In this regard, \citeN{Alvisi.2017} noted that a distinction has to be made between prediction uncertainty and emulation uncertainty. The former denotes the uncertainty associated with the natural variability of the true water demand process and describes the actual quantity of interest for practitioners. The latter denotes the uncertainty which arises and cascades within the data collection and modelling procedure (e.g., measurement/-data uncertainty, parameter uncertainty, or model structure uncertainty), as outlined by \citeN{Hutton.2014}. In contrast to prediction uncertainty, emulation uncertainty must be quantified but marginalized, so that the probabilistic forecaster issues only the expected natural variability of the true water demand process. Here, \citeN{Alvisi.2017}, \citeN{Gagliardi.2017}, and \citeN{Chen.2018} have published promising approaches.\\ However, to date the need for modelling the correlation structure between single hours in a multi-step-ahead forecast has not been addressed. The same applies to the need for appropriate evaluation criteria, which should likewise be able to penalize errors in the mean, the marginal properties and the correlation structure.\\ The practicality of providing a complete probabilistic multi-step-ahead forecast can be illustrated by the planing and management of storage capacities. Here, decision makers are interested in the aggregated demand, so that a statement can be made about the probability with which a water storage capacity can guarantee the supply over a certain period of time. This information forms the foundation to better balance demand peaks and to better schedule the pumping arrangements to take advantage of the electricity price structure. Furthermore, the evaluation of the correlation structure requires more sophisticated evaluation measures. To illustrate the limitations of evaluation measures used so far in the water demand forecasting literature, four hypothetical point forecasts are introduced in Fig. (ref), as done in similar way by \citeN{Haben.2014} in the field of energy demand modelling. As evaluation criterion the well-known mean absolute error (MAE) is applied. In the context of using the forecasts for water storage optimization to balance the expected demand, forecast (a) provides the best fit. Forecast (b) is moderate, and both forecasts on the bottom (c) and (d) are rather poor. Even so, forecast (b) does not exactly hit the peak; it is only slightly shifted, while both forecasts on the right substantially miss the demand peak. However, by considering the MAE, forecast (b), which was assessed as moderate in terms of storage optimization, achieves together with forecast (c) the worst score. Here, the authors introduce, for the following model evaluation, the energy score (ES) as an appropriate evaluation criterion to adequately account for the correlation structure over time, beside the mean and the marginal properties.\\ Focusing on the applied models and methods in the short-term water demand forecasting literature, it is striking that a vast variety of methods and modelling techniques have already been applied to predict future water demand. Initially linear regression and time series models were used, as outlined by \citeN{Herrera.2010}. However, with advances in the field of machine learning, various types of methods such as artificial neural networks, for example applied by \citeN{Adamowski.2010}, \citeN{Bata.2020}, \citeN{Guo.2018}, \citeN{Ghiassi.2008}, and \citeN{Anele.2017}, support vector regressions used by \citeN{Brentan.2017}, \citeN{Msiza.2008} and \citeN{Shabani.2016}, and random forests, as applied by \citeN{Herrera.2010} and \citeN{Chen.2017}, were successfully introduced. Furthermore, also combinations of the above methods, so-called hybrid-methods have received substantial attention, as applied for example by \citeN{Ambrosio.2019}. Nevertheless, linear regression and time series models are still ranked among the most popular modelling methods, as applied in the recent past by \citeN{Arandia.2016}, \citeN{Chen.2014}, \citeN{Caiado.2010}, and \citeN{Chen.2018}. However, as noted by \citeN{Pacchin.2019} and \citeN{Ghalehkhondabi.2017}, it is still difficult to pick a single method as the overall best, so that the performance of forecasting models based on different forecasting techniques are comparable. \\ Considering the stylized facts of the water demand process, non-stationary and non-linear behavior due to multiple seasonalities, autoregressive and external effects (as for example calendar and weather effects) represent the major challenge for modelling the water demand process and hence, require sophisticated models and methods, as noted by \citeN{Herrera.2010} and \citeN{Adamowski.2012}. Following \citeN{Romano.2014}, in practice parsimonious models, which are able to efficiently adapt to ever-changing operating conditions by applying self-learning ability, are preferred. Hence, in the recent literature, especially non-linear models with a low-dimensional feature space have been applied to model efficiently the complex structure of the water demand process, for example, the artificial neural networks used by \citeN{Pacchin.2019} and \citeN{Alvisi.2017}.\\ In this paper, the authors choose a rather different approach and introduce a high-dimensional feature space in a linear modelling framework. By taking advantage of the least absolute shrinkage and selection operator (lasso), introduced by \citeN{Tibshirani.1996}, the feature space can be automatically tuned, so that an efficient, parsimonious, simple interpretable and fast computable forecasting model is obtained. Hence, the model is well suited for application in real-time operating conditions. Moreover, an appropriate multi-step-ahead forecasting framework, to issue a complete multidimensional probabilistic forecasting distribution, is introduced and appropriate point and probabilistic forecasting measures are used to assess the forecasting performance. In this regard, the strictly proper ES can be highlighted, as it allows for simultaneously penalizing errors in the mean, the marginal properties and the correlation structure of an issued complete probabilistic multi-step-ahead forecast. The paper is organized as follows. Section 2 describes the data and discusses the stylized facts of the water demand process; Section 3 presents the time series model; Section 4 introduces the lasso estimation method; Section 5 outlines the forecasting procedure and presents the applied benchmark models; Section 6 presents appropriate point and probabilistic forecasting evaluation measures and introduces the ES as a strictly proper scoring rule; Section 7 summarizes the results for the calibration and validation period, interprets the proposed model, and discusses the practicality of complete probabilistic multi-step-ahead forecasts; and Section 8 concludes the paper.
The model is fitted to the hourly water demand data (2012-2018) of a water supplier in the western part of Germany. The data cover six years: four years are used for the model calibration (training) and two years for the model validation (testing). The calibration data are visualized in Fig. (ref) (a) and (b). The water distribution system under consideration provides fresh water from one water utility to around one million customers, who use the water for agricultural, private and industrial purposes. The data were not manipulated, so that outliers have not been replaced or transformed. Only the clock change was adjusted to simplify the data analysis.
The water demand process is demand driven and characterized by regional features, which may vary both spatially and temporally from region to region, as outlined by \citeN{Hutton.2015}. As highly influencing factors, the climatic and geographic conditions, as well as the commercial and social conditions of people have been identified, as mentioned by for example \citeN{Anele.2018}. These factors are responsible for the typical observable patterns on a daily, weekly, and an annual basis, as noted by \citeN{Adamowski.2010}. Considering the daily structure illustrated in Fig. (ref) (c), a trend in the mean but also a varying variance in the day course are observable. The weekly structure illustrated in Fig. (ref) (d) is characterized by full working day effects from Tuesday to Thursday and weekend effects on Saturday and Sunday. Friday and Monday are transition days; they are characterized by both, weekend and full working day effects. Considering the annual patterns, two characteristics are striking. First, the occurrence of holidays has to be highlighted: these days influence the behavior of the demand dramatically, as noted by \citeN{Hutton.2015} and depicted in Fig. (ref) (d). The authors classified public holidays in so-called fixed weekday holidays (FWH) and fixed date holidays (FDH), as proposed by \citeN{Ziel.2018}. The former always occurs on the same weekday, but on varying dates, and the latter always occurs on the same date, but on varying weekdays. Second, meteorological effects might also have a significant influence, as the water demand depends on gradual changes of weather conditions in the year, as illustrated in Fig. (ref) (a). Here, a weak seasonal course is observable; especially striking is the high demand in April and the low demand in November. An obvious explanation might be the impact of the growth phase in agriculture on the water demand.
To appropriately account for the high-dimensional complexity of the water demand process, the authors propose a model which is characterized by a huge feature space. However, by applying an automatic shrinkage and selection estimation method, the feature space can be tuned, so that a sparse, simple interpretable and fast computable forecasting model is obtained.\\ As inputs, the authors consider autoregressive effects, multiple seasonal patterns and calendar effects. Weather inputs are not considered, as the residual diagnostic indicates no significant effect. This might be explained by the fact that autoregressive effects to some extent already capture the influence of past external inputs such as the gradual increase of temperature over time. However, as indicated by \citeN{Bakker.2014}, an improvement of the forecasting performance by introducing weather inputs is feasible. Therefore, more sophisticated modelling approaches and the inclusion of additional weather forecasts might be required. This, in turn, would imply a significant increase in model complexity and in computational time, which might not be balanced out by the expected benefits, as the proposed model already indicates a good forecasting performance. Hence, the authors refrain in doing so and leave this issue for further research. The proposed model can be defined as
where $Y_t$ is the water demand at hour $t$ and $\mathbb{Y}_t = (Y_{t-1},Y_{t-2},...) $ denotes the past realizations. The model contains a deterministic component $\mu(t)$, an autoregressive component $\Phi(\mathbb{Y}_t)$, a time-varying autoregressive component $\Psi(t,\mathbb{Y}_t)$ and a zero mean noise process, so that $\mathbb{E}[\epsilon(t)]=0$. The $\epsilon(t)$ component accounts for the stochastic nature of the process, especially for the time-dependent variance structure. To provide a more thorough understanding of the proposed forecasting model, each component of equation ((ref)) is described in the following subsections.
As noted in section 2, the water demand process is characterized by strong periodic patterns on a daily, weekly and an annual scale. To control for these non-stationary features, the proposed model includes a deterministic component $\mu(t)$, which varies with time $t$. It can be defined as
where $\mu_0$ denotes a constant. The daily and weekly structure is modelled by the hour of the day $HoD$ and the hour of the week $HoW$ dummy functions, respectively. As the lasso approach is especially sensitive to changes over time, cumulative dummies are also introduced. For illustration purposes, the $HoD$, $HoD^{\text{cum.}}$, and $HoW$ dummies are illustrated in Fig. (ref) for a period of two weeks.\\ The annual structure is modelled by a public holiday component $\Pi(t)$ and a meteorological component $\Upsilon(t)$, which are defined in the two equations below. The modelling of public holidays poses a major challenge in water demand forecasting, as each holiday is characterized by an individual structure and simultaneously by a rare occurrence. Hence, the authors defined the public holiday component $\Pi(t)$, so that the individual daily structure is reasonably covered, but the risk of over fitting is rather low. The public holiday component is defined as
where $HD(t)$ denotes a dummy function, which models every single public holiday as a function of the day rather than a function of the hour. This allows for the individual characteristic of each day. To adequately capture the hourly characteristics, two additional dummy functions are introduced, namely the $FWH^{\text{cum.}}(t)$ and the $FDH^{\text{cum.}}(t)$ dummy function. The former is the cumulative hour of the day dummy function of the fixed weekday holidays and the latter of the fixed date holidays. By grouping the holidays into two classes, the authors are able to model the hourly patterns on similar holiday types. To illustrate the public holiday component $\Pi(t)$, each sub-component is visualized for an example period of two weeks in Fig. (ref).\\ The meteorological component is modelled by a linear combination of $K$ basis functions, as proposed by \citeN{Ziel.2016b} in the framework of electricity demand forecasting. This component can be defined as
where $\upsilon_{i}$ denotes the parameter for each cumulative basis function $B_i^{\text{cum.}}$. This approach enables the modelling of recurring events in a smooth manner. It increases the flexibility but simultaneously guarantees interpretability of the model. Focusing on the meteorological cycle $K = 4$, cumulative basis functions are used to model each season of the year in a local set-up, as a cold winter, for example, does not necessarily indicate a cold summer, and vice versa. For a detailed explanation of how the basis functions are computed, compare Appendix (ref) and \citeN{Ziel.2016c}, respectively. To illustrate the considered standard and cumulative basis functions, an example period of one year is visualized in Fig. (ref). Here it can be noted that other approaches such as Fourier series or wavelet compositions might lead to similar results, as applied by \citeN{Alvisi.2007} and \citeN{Adamowski.2012}.
The autoregressive component is the short-term memory of the process and for the proposed model the most valuable component, as highlighted in Fig. (ref), it can be defined as
where $\phi_k$ is the parameter of the past water demand $Y_{t-k}$ from the index set $\mathbb{I}$. Although the estimation method described below is able to automatically select relevant features from a huge feature space, a pre-selection of index sets is required. The pre-selected lags are presented in Table (ref).\\ To account for the linear dependency structure of the recent past, the first 361 lags are included. Moreover, seasonal autoregressive components are likewise included to account for the influence of previous weeks and months. The pre-selection of lags is based on applying statistical instruments such as the autocorrelation function (ACF) and partial autocorrelation function (PACF).
The time-varying autoregressive component $\Psi(t,\mathbb{Y}_t)$ is in fact an interaction term, containing the most relevant features of the deterministic and autoregressive part and can be defined as
where $\psi_k(t)$ denotes the interaction term of lags of index set $\mathbb{S}$, $\mathbb{S} \subseteq \mathbb{I}$ and the deterministic dummy functions of $HoD(t)$ and $HD(t)$. The selection of the interaction terms was based on the findings of the automatic parameter selection and shrinkage algorithm of the lasso estimator for the deterministic and autoregressive components.
As noted in section (ref) and depicted in Fig. (ref) (c), the water demand process suffers from heteroscedasticity. To account for the conditional time-varying variance, the error term $\epsilon(t)$ is modelled as a function of time, assuming that the variance in time point $t$ depends on the variance of past hours. Here, a time-varying ARCH(p) process is used, as suggested by \citeN{Ziel.2015} in the context of energy market modelling. Then the error term $\epsilon(t)$ can be written as\\
\centerline{ $\epsilon(t) = \sigma(t)Z_t$ where $Z_t$ is i.i.d. with $\mathbb{E}(Z_t)=0$ and $\mathbb{V}$ar $(Z_t)=1.$}
The conditional variance $\sigma^2(t)$ is modelled as
where $\alpha_k$ is $\geq0$. This is required as the variance can not become negative. To model the autoregressive structure, the authors use the lags from index set $\mathbb{K}$. As $\mathbb{K} \subseteq \mathbb{I}$, the lag pre-selection is based on the pre-selection in equation ((ref)). The deterministic part $\alpha_0(t)$ is defined as
where $\alpha_0(t)$ equals the deterministic component of equation ((ref)), adjusted for the cumulative components. This is understandable, as changes over time require the existence of negative parameters, which in turn are ruled out by the non-negativity constraint. Hence, $\Gamma(t)$ is the non-cumulative public holiday component, which is defined as
and $\Xi(t)$ is the non-cumulative meteorological component, which is defined as
As mentioned above, the authors use the lasso estimator and its properties for handling huge feature spaces to obtain a parsimonious and efficient model. The lasso algorithm is able to distinguish between relevant and irrelevant features with the so-called selection ability. Moreover, the lasso algorithm is also able to weight the features in accordance to their explanatory power for the independent variable. This is denoted as the shrinkage ability. Hence, irrelevant features will be fully excluded and less important features lose influence, as outlined by \citeN{Hastie.2015}. The lasso is the standardized version of the well-known ordinary least squares (OLS) estimator, extended by a penalty term. The lasso optimization problem can be defined as
where $\boldsymbol{\tilde{Y}}$ denotes the standardized water demand vector, $\boldsymbol{\tilde{X}}$ the standardized feature matrix and $\widehat{\boldsymbol{\tilde{\beta}}}$ the estimated standardized lasso parameter vector. For illustration purposes, in Table (ref) the feature matrix for the conditional mean estimation and in Table (ref) the feature matrix for the conditional variance estimation are depicted. The tuning parameter $\lambda$ regulates the impact of the penalty term. In case $\lambda$ converges to zero, the penalty term becomes meaningless and the standard OLS solution is obtained. However, as $\lambda$ approaches one, the penalty term forces the estimation parameters one by one to become exactly zero. This solution would imply that no parameters are included in the model. To choose the appropriate tuning parameter $\lambda$, a selection algorithm is required. Here, the authors use the Bayesian information criterion (BIC), which is especially suitable for huge feature spaces as it is more conservative than, for example, the well-known Akaike information criterion (AIC). This is a valuable property, as especially the risk of over fitting increases with a high model parametrization. For further details about the BIC, see \citeN{Neath.2012}.
For computation, the authors apply the glmnet package in R. It is based on the coordinate descent algorithm, which enables a fast and efficient parameter estimation. For further explanation, see \citeN{Friedman.2010} and \citeN{Hastie.2015}.
As noted by \citeN{Alvisi.2017} the quantity of interest in probabilistic forecasting is the natural variability of the true water demand process and not the uncertainty arising and cascading in the modelling procedure. Based on this, complete probabilistic multi-step-ahead forecasts are constructed, which consider beside the mean and the marginal properties, especially the correlation structure between each issued hour $h \in H$, so that a joint distribution rather than a marginal distribution is issued. This implies that the expected water demand for a fixed forecasting horizon $H>1$, can be seen as a multivariate random variable $\boldsymbol{Y} \in R^H$, where $\boldsymbol{Y}$ follows an unobservable distribution $G$, so that $\boldsymbol{Y} \sim G$.\\ To appropriately model the multivariate random variable $\boldsymbol{Y} \in R^H$, ensemble forecasts reveal preferable properties. In the water demand forecasting literature, ensemble forecasts are commonly used, for example by \citeN{Tiwari.2013} and \citeN{Hutton.2014}. However, these forecasts have predominately been applied to account for uncertainties arising in the emulation process.\\ In this research study, ensemble forecasts are used to model the prediction uncertainty. As forecasting horizon, $H=24$ is chosen and an ensemble is created by recursively solving the corresponding forecasting model in a Monte-Carlo simulation with a total of $M=1,000$ sample paths. To illustrate the meaning of the correlation structure, different model simulations are depicted in the left column of Fig. (ref). In addition to the standard procedure (a), three manipulated correlation structures are introduced. The comonotone model simulation is characterized by a perfect positive path-dependency (b), so that the sample paths are non-intersecting, that is, the path that is highest in the first hour remains highest for the following hours and so on. The countermonotone model simulation (c), by comparison, considers perfect negative pairwise dependency, so that at each hour the sample paths run contrary to the hour before, that is, the path that is highest in the first hour is lowest in the second and so on. Finally, the independent model simulation (d) is characterized by an independent path-dependency, so that each hour of a sample path is arbitrarily connected to the adjacent hours.\\
The forecasting study is based on $1,000$ equally distributed forecasting time points drawn from the validation data to mimic real-time conditions. To account for ever changing operating conditions a rolling window approach is used to ensure a recurring parameter re-estimation at every forecasting time point, as suggested e.g. by \citeN{Pacchin.2017}. Here, it might be worth noting that the forecasting time point is not starting every corresponding day at midnight but instead on varying hours over the day.
As benchmarks, different forecasting models for hourly water demand data from the literature are applied. First, two SARIMA(0,1,4)(0,1,1) models with a seasonal period of 24 and 168 are considered, based on \citeN{Arandia.2016}. As the computational time with a calibration data length of up to four years distinctly exceeds the requirements for real-time applications, the authors decided to shorten the length of the calibration data for all considered competitors from the literature. Here, \citeN{Arandia.2016} have suggested various data lengths; for the corresponding data a length of 28 days reveals the best results. The time series models are estimated and predicted by applying the forecast package in R, for further explanations see \citeN{Hyndman.2008} and \citeN{Hyndman.2019}.\\ Second, three typical machine learning methods are used, namely an artificial neural network $ANN_{Her}$, a support vector machine $SVM_{Her}$ and a random forest $RF_{Her}$. All three models are based on \citeN{Herrera.2010}. As features, the authors use the lags $1,2,24,168$ and a Fourier series of order $1$ for the daily and annual cycles, respectively. The model tuning parameters and the applied length of the calibration data is chosen in accordance with \citeN{Herrera.2010}, whereby the best model is chosen ex post. The support vector regression is estimated by applying the e1071 package, the neural network is estimated by applying the nnet package, and the random forest is estimated by applying the svm package in R, for further explanations see \citeN{Meyer.2019}, \citeN{Venables.2002}, and \citeN{Liaw.2002}, respectively.\\ Third, a neural network $ANN_{Pac}$ based on \citeN{Pacchin.2019} and \citeN{Alvisi.2017} is applied. A log sigmoid transfer function is used in the hidden layer and a pure linear one in the output layer. The length of the calibration period and the features are chosen based on \citeN{Pacchin.2019}. The number of hidden neurons is set equal to 72 in accordance to \citeN{Alvisi.2017} and the learning rate varied between 0.0001, 0.001, 0.01, and 0.1. Here, the best results are obtained expost for a learning rate of 0.01. The neural network is also estimated by applying the nnet package in R.\\ Fourth, three naive models are applied, namely a naive mean model $Naive_{Mean}$ as used by \citeN{Pacchin.2019} and \citeN{Gelazanskas.2015}, a naive mean model conditioned on the type of the day $Naive_{FM}$, and a naive mixed random walk model $Naive_{MRW}$. The $Naive_{Mean}$ model computes the values of the forecast by taking the arithmetic average of each hour of the day. The $Naive_{FM}$ computes likewise the $Naive_{Mean}$ model the forecast based on the arithmetic average of each hour of the day, however, with distinction in the type of the day. Here, the authors have distinguished in Mondays, Tuesdays to Thursdays, Fridays, Saturdays, Sundays, and Holidays. The $Naive_{MRW}$ is the only non-parametric model and is defined as
where $\epsilon_t$ denotes the error term. Finally, also two $AR(p)$ processes are applied. As demonstrated by \citeN{Ziel.2015} in the field of electricity price forecasting, such models are easy to implement and their predictive accuracy is reasonable good. To provide an easily implementable and fast computable forecasting benchmark model in the field of short-term water demand forecasting, the authors have attached the corresponding implementation in R in Appendix (ref).\\ The $AR(p)$ benchmark models differ in the mean adjustment, so one part includes a mean adjusted hour of the day $D$ process and the other a mean adjusted hour of the week $W$ process. The considered $AR(p)$ processes are estimated by the Yule-Walker estimator and computed by applying the stats package in R, for further explanation see \citeN{R.2019}. The $AR(p)$ models are given by
where $i \in \{D,W\}$. The autoregressive structure is captured by including lags of order $p$, which in turn are chosen by applying the BIC. Here, the maximum order $p$ is set equal to 1,500. All introduced benchmark models are based on a bootstrap distribution assumption.
\Section{Evaluation criteria} As the evaluation procedure is of particular relevance to identify the best forecasting model, the choice of an appropriate criteria should always be done in accordance with the forecasting purpose at hand. In the framework of point forecasting evaluation, a broad variety of measures is available. Here, the water demand forecasting community is mainly focused on measures such as the MAE, the mean squared error, the Nash-Sutcliffe model efficiency coefficient (NS) and generalizations of the previously mentioned measures as applied by, for example \citeN{Anele.2018}, \citeN{Donkor.2014}, \citeN{Herrera.2010}, and \citeN{Brentan.2017}. Here, it might be worth noting that although the issued forecasts have a probabilistic distribution, they can easily be reduced to simple point forecasts. Here, the MAE and the root mean squared error (RMSE) can be highlighted, as they reveal preferred properties, as discussed by \citeN{Franses.2016}. The MAE is a strictly proper criterion for the median and the RMSE is a strictly proper criterion for the mean. In this context, "strictly proper" refers to the ability that only the perfect forecast minimizes the named criterion, as discussed by \citeN{Gneiting.2014}.\\ As noted in the introduction, probabilistic forecasting is gaining acceptance. Hence, probabilistic evaluation measures are also required. As pointed out by \citeN{Gneiting.2007}, the aim of probabilistic forecasting is to maximize the sharpness of probabilistic forecasts subject to calibration. Where, calibration denotes the statistical consistency between the issued distribution and the events that materialize and sharpness refers to the concentration of the predictive distribution. To adequately assess the properties of a probabilistic forecast, strictly proper scoring rules are well suited. They assign a numerical score based on the predictive distribution and the events that materialize to the corresponding forecast. Here, the authors introduce the pinball score (PB) as an appropriate scoring rule to evaluate the marginal distribution of an issued forecast. Moreover, as illustrated in Fig. (ref), the dependency structure is likewise relevant, especially if the forecast is to be used in the field of storage capacity optimization. Here, the authors introduce the ES as an appropriate scoring rule.
As point forecasting measures the MAE, RMSE, and the NS are applied. They can be defined as
where $Y_{t_j+h}$ denotes the true water demand at forecasting time point $t_j$ and hour $h$, $\bar{Y}_{t_j}$ the mean of the true demand vector $Y_{t_j} = Y_{t_j+h},...,Y_{t_j+H}$, and $\hat{Y}_{t_j+h}$ the estimated water demand at forecasting time point $t_j$ and hour $h$, $H$ is the number of hours issued in a forecast and $N$ the number of forecasting time points $t_j$ in the forecasting study. These named measures are suited to point forecasting, but are not appropriate in terms of probabilistic forecasting, so that more sophisticated criteria are required.
A promising scoring rule, which appropriately considers the marginal properties of probabilistic forecasts is the PB. It is defined as
where $\hat{q}_{t_j+h,\tau}$ is the issued quantile and $\tau$ the corresponding quantile level. $Y_{t_j}$ is the true water demand at forecasting time point $t_j$, $L$ denotes the number of quantile levels in the dense grid $\boldsymbol{\tau}$, $N$ is the number of forecasting time points in the forecasting study, $H$ is the number of hours issued in a forecast and $\mathbbm{1}$ is an indicator function. Note that the PB converges to the continuous ranked probability score (CRPS) for an infinitely dense and equidistant grid $\boldsymbol{\tau}$. The CRPS in turn is a widely-accepted strictly proper scoring rule for the full distribution function and a generalization of the well-known MAE, as introduced in equation ((ref)). Hence, for the special case $\tau=0.5$, the PB is a scaled version of the MAE. For further information on the PB, see \citeN{Nowotarski.2018} and \citeN{Steinwart.2011}.\\ However, beside the marginal properties in a multi-step-ahead forecasting framework, also the correlation structure between single hours in a corresponding forecasting horizon is highly relevant, as already noted. In this regard, the PB has limitations, as the evaluations of the issued quantiles of each hour are assessed separately. Currently, only a few scoring rules are known, which are also able to account for the path-dependency between hours. In this regard promising candidates are the logarithmic score, the David-Sebastiani score, the variogram score, and the ES, as highlighted by \citeN{Scheuerer.2015} and \citeN{Gneiting.2007}. Focusing on the forecasting purpose of this research, the ES has preferable properties, as it is a strictly proper scoring rule for probabilistic multi-step-ahead forecasts and it is well suited to deal with the issued sample paths of the ensemble forecast. The ES can be defined as
whereby $F$ denotes the issued multivariate distribution of the forecaster, $\boldsymbol{\mathrm{P}}$ and $\boldsymbol{\mathrm{P}}^{\prime}$ are independent random draws of $F$, $\boldsymbol{\mathrm{y}}$ is the considered observation vector of the considered process, and $\|\cdot\|_2$ denotes the Euclidean norm. The ES is like the PB, a generalization of the CRPS and equals the CRPS in the case of $H=1$. For detailed explanations, see \citeN{Gneiting.2007}, \citeN{Gneiting.2008}, and \citeN{Scheuerer.2015}. In this setting, the ES is estimated by
where $\mathbf{Y}_{t_j,H} = (Y_{t_j +1},...,Y_{t_j+H})$ denotes the water demand vector at the forecasting time point $t_j$. $\mathbf{\widehat{Y}}_{t_j,H}^{[i]}$ and $\mathbf{\widehat{Y}}_{t_j,H}^{[l],\prime}$ denote the corresponding $i$-th and $l$-th simulation vector of the simulation matrices $\mathbf{\widehat{Y}}_{t_j,H}$ and $\mathbf{\widehat{Y}}_{t_j,H}^{\prime}$, respectively. $M$ denotes the number of simulations, $H$ the number of hours issued in a forecast and $N$ the number of forecasting time points in the forecasting study.\\ In the following the ES is applied as the determining measure. However, also the PB, the MAE, RMSE, and NS are listed. Although the ES is able to discriminate errors in the mean, the marginal properties and the correlation structure simultaneously, pinpointing specific causes of a poor performance is rather difficult. Here, the PB, the MAE, and the RMSE might provide additional information in terms of diagnostic checking.
To identify differences in the forecasting performance of competitors in a statistically reliable way, the Diebold Mariano (DM) test is applied. As noted by \citeN{Nowotarski.2018}, the DM test is more popular in the framework of point forecasting, but is also applicable to probabilistic forecasting. The DM test is simply an asymptotic $z$-test with the hypothesis that the mean of the loss differential series is zero. On the assumption that the loss differential is covariance stationary, the DM test is asymptotically standard normal.\\ The authors calculated two one-sided DM tests with a significance level of 5 $\%$, to test for significant differences from zero in both directions. For further details of the DM test, see \citeN{Diebold.1995}.
In this section the forecasting performance of the proposed forecasting model $ARXARCHX_{lasso}$ from equation ((ref)) is presented in comparison to the introduced benchmark models for the calibration and validation data. The results are presented in Table (ref) and (ref). As benchmark model the $AR(p)^W$ model is chosen to compute the improvements in $\%$. Moreover, the diagnostic checking of the $ARXARCHX_{lasso}$ is shown for the calibration period to pinpoint further improvements and to highlight the most influential parameters. Finally, in the discussion section the practicality of complete probabilistic multi-step-ahead forecasts is outlined and illustrated by a storage management problem.
By examining the forecasting results in terms of the MAE, RMSE, and NS for the calibration period in Table (ref), the forecasting results of the $SARIMA(0,1,4)(0,1,1)$ models are striking. Here, the low scores in terms of the MAE and RMSE and the high scores in terms of the NS are suspicious and suggest model over fitting. This is most likely caused by the fact that a too short calibration period was chosen. Hence, the extension of the calibration period for the named models can be recommended. Focusing on the remaining forecasting models, the proposed $ARXARCHX_{lasso}$ is ranked the best, followed by the proposed simple $AR(p)$ models. The $Naive$ models are ranked the worst. The rather poor performance of the $ANN_{Pac}$ in comparison to the other sophisticated models is caused by the chosen 24-dimensional output layer of the network, which results in a 24-dimensional forecast in the calibration period instead of a one-step-ahead forecast as issued by the remaining sophisticated models. This allows the $ANN_{Pac}$ to better approximate the forecasting error of the validation period.
The obtained scores in terms of the ES in Table (ref), show that the proposed forecasting model $ARXARCHX_{lasso}$ dominates all considered benchmarks significantly. For illustration, also the comonotone ($^{*}$), the countermonotone ($^{**}$), and the independent ($^{***}$) model simulations of the $ARXARCHX_{lasso}$ are presented. Here, it can be demonstrated, that neither the MAE, nor the RMSE, nor the NS or the PB are able to detect the manipulation in the dependency structure. However, although the ES is the determining measure in terms of evaluating a complete multidimensional distribution, pinpointing specific causes of a poor performance is rather difficult. Here, the PB, MAE, and RMSE allow a more thorough analyses. Focusing on the forecasting accuracy with expanding forecasting horizon from one up to 24 hours in Fig. (ref) (a)-(c), it is observable, that with expansion of the horizon the forecasting error for all considered models increases notably. Here, the $ARXARCHX_{lasso}$ is constantly ranked the most accurate and reveals a moderate growth rate over time, compared to the benchmark models for all three criteria. Hence, the authors conclude that also for smaller forecasting horizons, the $ARXARCHX_{lasso}$ model is well suited. Focusing on Fig. (ref) (d), the $PB(\tau)$ is plotted across the applied probability grid $\boldsymbol{\tau}$. Here, the $PB(\boldsymbol{\tau})$ of equation ((ref)) has been computed without averaging across the quantile levels, so that explicitly the marginal properties in terms of the probabilities can be evaluated. It can be seen that the $ARXARCHX_{lasso}$ achieves the best score for each probability. Hence, the superior performance of the $ARXARCHX_{lasso}$ can be verified at once in terms of the ES, but also for the mean and the marginal properties separately by considering the MAE, RMSE, and PB.\\ Focusing on the considered benchmark models, the forecasting accuracy of the $SARIMA(0,1,4)(0,1,1)$ models has distinctively dropped compared to the calibration period. Here, beside a too short chosen calibration period, the authors suspect that the application of both an ordinary and a seasonal differencing operator might cause the rather poor performance. As most accurate competitors the proposed simple $AR(p)$ models can be highlighted and hence, the authors suggest them as easy implementable and fast computable benchmarks in the field of water demand forecasting.
As the proposed forecasting model $ARXARCHX_{lasso}$ is constructed in a linear modelling framework, the authors are able to make inferences and hence can pinpoint the importance of each considered parameter for the water demand process. Moreover, by doing residual diagnostic shortcomings of the model can be outlined, so that further improvements can be suggested.\\ Considering the ACF of $\widehat{\epsilon}$ in Fig. (ref) (a), it is observable that the autocorrelation structure is sufficiently captured by the conditional mean estimation. Thus, without taking extra external information into account, no major forecasting improvements can be expected. Considering the ACF of $\widehat{\epsilon}^2$ in Fig. (ref) (b), it becomes obvious that volatility clusters exist. Therefore, the authors introduced equation ((ref)), to account for the time-varying variance. As illustrated in the ACF of the residuals of the conditional variance estimation $\widehat{\omega}$ in Fig. (ref) (c), the volatility clusters could be mostly captured. However, the strong first lag of the ACF of the squared residuals of the conditional variance estimation $\widehat{\omega}^2$ in Fig. (ref) (d), is suspicious and indicates that still further improvements are achievable. In this context, the authors would suggest introducing an absolute instead of a squared function for modelling the volatility cluster. This might be a promising task, as the residuals of the water demand process turned out to be heavy-tailed. Moreover, an improvement in accuracy might also be achievable by applying an iteratively reweighted adaptive lasso algorithm estimation procedure, as done by \citeN{Ziel.2016}, to better deal with the conditional heteroscedasticity in the high-dimensional setting.\\ Concentrating on the huge feature space, the linear forecasting framework allows for pinpointing the most influential features. As an automatic shrinkage and selection operator is applied, only a small fraction of the initially introduced 1,468 features is selected and only a sub fraction is considered as highly influential, as illustrated in Fig. (ref) and Fig. (ref). In this context, it has to be noted that the parameters are scaled, so that the influence of the features is illustrated in relation to each other and not their absolute effects.\\ For the conditional mean estimation in Fig. (ref), only 147 out of 856 features have been selected. Here, it is striking that all components contain influential features, what in turn might justify the consideration of the rather high-dimensional feature space to model the water demand process. By highlighting only the most important features, the authors can name hour 23 of the day and the change over time of hour 4 of the day; hour 127 and 129 of the week; hour 7 of the fixed weekday and fixed date holiday, respectively, as well as the autoregressive lags 1,2,23,24,25,167,168, and 169 and the interaction of lag 2 and hour 7. Overall, the autoregressive component can be identified as the most influential. This might also explain the superior performance of the proposed $AR(p)$ models. For instance, the $AR(p)^W$ model has selected in average an order of $p=233$, so that the previous 233 hours have been used. Considering the applied forecasting models in the literature, it can be said that most models have considered autoregressive effects, for example done by \citeN{Pacchin.2019}, \citeN{Anele.2017}, \citeN{Arandia.2016}, and \citeN{Chen.2018}. Moreover, deterministic patterns on the daily, weekly, and annual scale have also regularly been considered, for example done by \citeN{Alvisi.2017}, \citeN{Gagliardi.2017b}, and \citeN{Brentan.2017}. However, as mostly models with a low-dimensional feature space have been used the effects have not been modelled in a high resolution. Here, the authors conclude that the difference in forecasting accuracy can be explained by the difference in the amount of initially considered information and the model itself to deal with the provided information. Besides, it might be worth noting that the number of parameters does not necessarily corresponds to the number of features. Hence, although the $ARXARCHX_{lasso}$ seems to be rather complex, the number of active parameters used in the final model is only moderate compared to the number of parameters of the introduced machine learning based models, as shown in Table (ref).\\ For the conditional variance estimation in Fig. (ref), only 61 out of 612 features are selected. Here, only the hour of the day, the hour of the week, the fixed date holiday and the autoregressive component are considered as influential. By highlighting the most important features of each component, the authors can name hour 8 of the day; hour 7 of the week; hour 7 of a fixed weekday holiday and the autoregressive lags 1,24,48, and 192. By examining the existing literature, only a few modelling approaches have addressed the time-varying variance structure, for example done by \citeN{Caiado.2010} and \citeN{Hutton.2015}.
Besides the mean and the marginal properties is the dependency structure of probabilistic multi-step-ahead forecasts of considerable relevance for practical applications, for example, with respect to time-dependent optimization problems in the field of water planing and management. For illustration, the proposed $ARXARCHX_{lasso}$ model and the manipulated versions of the $ARXARCHX_{lasso}$ model of Table (ref) are considered. They are identical in the mean and the marginal properties but differ in the simulated dependency structure, as highlighted in Table (ref) and displayed in Fig. (ref). Focusing on the accuracy of the simulated dependency structures, the $24\times24$ correlation matrix of the corresponding model simulation (lower triangular) compared with the correlation matrix of the true water demand (upper triangular) is displayed in the middle column of Fig. (ref). It is obvious, that only the standard model simulation (a) covers the true dependency structure of the water demand process reasonably. To verify the practicality, the water storage optimization as an actual operative management is considered. It forms the foundation to better balance demand peaks, to increase the security of supply, and to better schedule the pumping arrangements to take advantage of the electricity price structure. Here, decision makers are especially interested in the expected aggregated or cumulative demand, so that a statement can be made about the probability with which a water storage capacity can guarantee the supply over a certain period of time. The impact of the dependency structure on the cumulative water demand is illustrated on the right column of Fig. (ref). Here, the density of the cumulative sample paths of each model simulation is compared with the true density of the cumulative water demand process.\\ Moreover, assuming a hypothetical water storage tank with a capacity of 290,000 m³ and a considered time period of 24 hours, the actual expected probability for each model simulation is computed and compared with the true probability, as shown in Table (ref). Here, the true cumulative water demand exceeds with a probability of 0.0935 the storage capacity within the corresponding 24 hours. The $ARXARCHX_{lasso}$ with a standard model simulation predicts a probability of 0.0934, the comonotone model simulation a probability of 0.1524, the $ARXARCHX_{lasso}$ with a countermonotone model simulation a probability of 0.0782, and the $ARXARCHX_{lasso}$ with an independent model simulation a probability of 0.0000. Hence, although all model simulations are simulated with identical mean and marginal properties, their expected cumulative water demand distributions differ distinctively and only the standard model simulation (a) of the proposed $ARXARCHX_{lasso}$ provides a reasonable forecast.\\ Finally, to get an impression of the forecasting performance of the $ARXARCHX_{lasso}$ model, in Fig. (ref), three forecasts of the validation period are presented. Beside two regular days (on the right and left side) also one fixed weekday holiday is illustrated. With respect to the presented holiday, it is striking that especially the early morning hours are well captured, as the demand distinctively drops below the usual demand on a weekday and even below the demand of the preceding Saturday and Sunday. Concerning the overall performance of the proposed model, the authors determined that even on average the forecasting performance is satisfactory; however special days such as specific holidays and other rarely occurring events are extremely difficult to capture. This is also valid for some periods during the summer. Here, the extension of the calibration period with a focus on rarely occurring events and holidays or the introduction of weather forecasts might lead to improvements. Nevertheless, on normal days the proposed model performs very well and is therefore well suited for a large part of the year.
This paper proposed a forecasting model with a high-dimensional feature space in a linear framework to capture the complex structure of the water demand process. By applying the least absolute shrinkage and selection operator (lasso), the model could be automatically tuned, so that a parsimonious, simple interpretable and fast computable forecasting model could be obtained. The model clearly outperformed existing benchmarks from the water demand forecasting literature. Here, the authors concluded that the complex water demand process can be more accurate modelled by automatically tuned linear models with a high-dimensional feature space than by applying rather computational complex non-linear models with a low-dimensional feature space.\\ Moreover, an appropriate multi-step-ahead forecasting framework to issue a complete multivariate probabilistic forecasting distribution, which is able to account also for path-dependencies, has been introduced and their practicality in the field of water storage management has been verified. Furthermore, the need for more sophisticated evaluation measures has been outlined and the ES as a strictly proper scoring rule has been proposed. The ES allows for penalizing errors in the mean, the marginal properties and the correlation structure of the corresponding forecast. Nevertheless, there is still much research to be done. Focusing on the proposed forecasting model, it might be worthwhile to include forecasts of weather inputs or further external deterministic processes such as the announced water withdrawal in advance of bulk customers. Concerning the already considered features, the extension of the time-varying autoregressive component and applying more advanced modelling techniques for holidays and the time-varying variance might be promising avenues to explore.
Some or all data, models, or code generated or used during the study are proprietary or confidential in nature and may only be provided with restrictions (e.g. anonymized data).