HOME | ÍNDICE POR TÍTULO | NORMAS PUBLICACIÓN

Espacios. Vol. 37 (Nº 22) Año 2016. Pág. 5

A study of the performance of individual techniques and their combinations to forecast urban water demand

Un estudio del rendimiento de las técnicas individuales y sus combinaciones para la previsión de demanda de agua urbana

José DONIZETTI de Lima 1; Gilson ADAMCZUK Oliveira 2;Marcelo Gonçalves TRENTIN 3; Dayse Regina BATISTUS 4; Cassiane Bayerl POZZA 5

Recibido: 25/03/16 • Aprobado: 13/04/2016

Contents

1. Introduction

2. Theoretical framework

3. Methodology

4. Results and Discussion

5. Conclusions

Acknowledgements

References

Appendix I - Complementary results

1. Introduction

Nowadays, approximately 36% of the world population or around 2.4 billion people live in water-scarce regions and 22% of the global gross domestic product (GDP) is produced in areas with problems related to water shortage, stress, or crisis. According to the analysis of the International Food Policy Research Institute (IFPR), if the current water management practice and consumption tendencies are maintained, 45% of the predicted GDP for 2050 will be at risk. This would be equivalent to 1.5 times the size of the current global economy. In addition, approximately 4.8 billion people (or 52% of the entire world population) will be exposed to serious water shortages by 2050 (VWNA, 2012). 

The search for a sustainable water consumption policy is a way to avoid those pessimistic projections, or at least to lessen their impact in the world.  The challenge of this remains in the alliance between the growth of the population and economy and the rational use of water resources. Operations research will play a major role in this task, and the authors regard the use of time series adjustment for monthly water demand in this context.

In comparison with the world's state of affairs, Brazil is in a comfortable situation regarding its water resources. The availability of water indicates an adequate situation when compared with the values from other countries listed by the United Nations (UN). However, in spite of this apparently comfortable situation, the spatial distribution of water resources in Brazil is unequal. More than 80% of the water resources are concentrated in the Amazon basin region, which has the lowest human population density in the country, along with a reduced demand for urban, industrial, and irrigation use. (ANA, 2012).

The water withdrawal volume corresponds to the amount of water obtained by water companies from rivers, large natural reservoirs, water storage dams, or other sources of water to be explored. Water withdrawal is performed by means of pumping, and the streamline is directed to the water treatment stations where chemical products are added to process the water.

Water companies need an adequate forecasting system to allow one to develop short-term, medium-term, and long-term planning. With this in mind, the manager can plan the company's premises, and allow the water-supplied community to adopt actions and develop its own plans to implement a sustainable use of this resource.

The goal of this paper is to propose a systematic approach for evaluating the performance of individual techniques and their combinations to forecast the withdrawal of water volume by water companies. Even considering that there are losses in the distribution of water, this volume represents the water demand in urban centers. In this work, the time horizon used in the study was 12 months.   Section 2 provides a review of the forecasting techniques applied to the water demand (2.1) and the growing tendency that is represented by the combination of forecasting techniques (2.2), although this is only used in urban water demand forecasting occasionally.  The combination of techniques aims at decreasing the residual variance and improving the forecasting accuracy (2.3). In section 3, a detailed systematic approach is proposed, being the original main contribution of this work. Another contribution of this paper is the wide range of combinations that are tested and evaluated comparatively, something unusual in time series literature (Martins & Werner, 2012). In section 4, the systematic approach is applied with a good level of details using a data set of 10 years from 10 cities of the southwestern region of the state of Paraná, Brazil.  Section 4 also provides a discussion of the results. Conclusions and considerations are presented in section 5.

2. Theoretical framework

A time series is a set of observations obtained in a sequential and systematic way over an interval of time. Among the time series, techniques presented in the literature to analyze a set of values obtained during a time interval, two classes of models play an important role: the exponential smoothing models (ES) and the autoregressive integrated moving average (ARIMA) models (Makridakis et al., 1998; Armstrong, 2001; Morettin & Toloi, 2006). On the other hand, the artificial neural networks used in the computational intelligence approach for time series have shown promising results (Calôba et al., 2002; Dias, 2006; Martins & Werner, 2012; Wu et al., 2013).

Regarding the techniques of artificial intelligence, Wu et al. (2013) recommend artificial neural networks (ANN) with a radial basis function due to their potential to approximate nonlinear behavior. According to these authors, the performance of ANN is influenced by the calibration of its parameters. Among these are: (i) the number of occurrence in the input vector (n), which characterizes the unit delays along the time series; (ii) the eg value, which is used to give an objective to the neural network; the network's training uses Root Mean Squared Error (RMSE) as a calibration parameter; and (iii) the width of the receptive fields (sc), which interferes with the outputs of the units contained in the hidden layer of the network.

According to Martins and Werner (2012), the use of a Box-Jenkins methodology (ARIMA) is reasonable since this class is formed by robust models which are consolidated and have a wide range of applications. Moreover, the class of ES models is important because of its simplicity, ease of computational implementation, and for the associated history of promising results, as in the study entitled "M-Competition" (Makridakis & Hibon, 2000; Nikolopoulos et al., 2011a).

The forecasting methods can be understood as the way one or more of such techniques are used. Therefore, a combination of forecasting methods is a part of a forecasting methodology. In the following sections, techniques and methods applied to the specific case of water demand (2.1) and a detailed evaluation of the forecasting combinations (2.2 and 2.3) are  given, since they are the core element of the forecasting systematic approach proposed in this work.

2.1 Forecasting urban water demand

Water demand studies have been considering hours, days, weeks, and months as short-term predictions, whereas years are considered long-term predictions. This cannot be understood using the terms used for forecast horizons in time series studies, where "short" and "long" are related to the distance from the origin of forecasts.

Short- and long-term water demand results have been studied over the last 50 years and such studies can be divided into five categories: regression analysis, time-series analysis, computational intelligence approach, hybrid approaches, and Monte-Carlo simulation approach (Qi & Chang, 2011). Among the five categories, only regression analysis gives the same importance to the short- and long-term ranges. The remaining categories have their applications concentrated in the short-term demand. In recent years, such categories have been used for the prediction of urban water demand, and therefore they are listed in Table 1. A growing tendency in the usage of hybrid strategies is also observed, in which at least two techniques or approaches are applied in the elaboration of forecasting methodologies. It is also worth noticing the growing importance of neural networks in the water demand forecasting field, either in an isolated manner or combined with other techniques.

Table 1: Works about urban water demand in the last five years

Source: Elaborated by the authors.

Regarding the accuracy criteria used in the validation or model fitting, Table 1 shows that there is a predominant set of criteria that are widely used in time series analysis and computational intelligence approaches. Mean Absolute Percentage Error (MAPE), Mean Squared Error (MSE), Root Mean Squared Error (RMSE) and Mean Absolute Error (MAE) appear to be consistently used in recent researches. Other criteria, such as the Standard deviation of the Absolute Percentage Error (SAPE), the Normalized Root Mean Squared Error (NRMSE) and the Fractional Standard Error (FSE) can be considered as extensions or derivations of the previously mentioned criteria. The correlation coefficient (R) and the determination coefficient (R2) are typically used in techniques involving regression, even considering that the R2 can be used in time series models. In the only work that uses Monte Carlos simulation, Almutaz et al. (2012) uses the standard deviation (s) to verify the accuracy of the predictions.

Still regarding the accuracy, the utilization of unusual measurements in time series analysis, such as the TSx, that can be seen in Firat et al. (2010), and E and B criteria used by Campisi-Pinto et al. (2012), can be represented by:

In this equation, TSx is the threshold statistic for a level of x%, where n is the total number of forecasts having a relative forecast error smaller than x%, and N is the total number of forecasts. According to Firat et al. (2010), TSx quantifies the consistency of the forecast errors for a particular model.

The parameter E is the Nash-Sutcliffe model efficiency, Y_i are the real values in period i, Y_p are the forecasted values in this period, and Y_m represents the average value of the series of data. E is widely used in hydrology because it measures the ability of the model to forecast values which are different from the average. The index B provides a good measure of whether the model is overestimating (B > 1) or underestimating (B < 1) compared to observed values. B = 1 indicates non-biased model performance (Campisi-Pinto et al., 2012).

Finally, it can be seen in Table 1 that the studies on water demand for urban use are more concentrated in short-term forecasts (75%), from hours to a few months. Short-term forecasts are important from the standpoint of the control of the  water company premises or equipment. One could even identify losses at a schedule point which shows a demand that is much different from what was expected. On the other side, from the strategic standpoint, the long-term projections (in years) correspond to a relatively low number (25%), and they may relate to the planning of the water companies' expansion or to the study of the possibility of water  shortage occurrences in urban centers (domestic/industrial).

There are a considerable number of applications using the combination of forecasting time series techniques (ES, ARIMA and ANN) in other areas (Calôba et al., 2002; Stock & Watson, 2004; Andrawis et al., 2011; Martins & Werner, 2012; Wu et al., 2013). For urban water demand, the adoption of this strategy is still incipient and constitutes an approach to be explored, especially if used for long-term forecasts (years). Qi & Chang (2011) have commented that such univariate models have low reliability to represent all the involved factors that affect the consumption of water. As a consequence, its use is concentrated in the short-term forecasts.

2.2 The combination of forecasts

Bates and Granger (1969) introduced the combination of forecasts when they were stimulated by the search for a technique that presented the lowest possible error. This method incorporates forecasts obtained from different techniques, with the goal of improving the accuracy of the forecasts while extracting useful and independent information from the available data (Panagiotopoulos, 2011).

The utilization of  forecast combinations is justified by the  facts that: (i) there is no perfect forecasting technique, since there is no way to acquire the reality exactly as it is; (ii) it's not known if a forecasting technique  will always present a better performance than other techniques  in all periods and forecast horizon intervals using the whole series of data (Flores & White, 1988); (iii) the accuracy of an individual forecasting technique is sensitive to several factors (Makridakis et al., 1998); (iv)  the choice of an individual forecasting technique from a set of available techniques is riskier than choosing a combination of them (Hibon & Evgeniou, 2005); and (v) the combination of forecasts can present forecasts which are more robust than those obtained by individual techniques (Souza, 2008).

Several studies have contributed to reinforce that  forecast combinations can bring significant benefits in terms of  forecast accuracy (Clemen, 1989; Makridakis et al., 1998; Andrawis et al., 2011). However, although there may be a concentration of efforts to verify which forms of combinations would be the most cohesive for the improvement of forecasts, there are no unanimous results (Martins & Werner, 2012).

According to Flores and White (1988), in order to use a combination of forecasts, the selection of the base forecast techniques is necessary to determine which forecasts to include in the combination.  One must then choose the combination methods, which is how the individual techniques will be combined.  For Werner and Ribeiro (2006), the choice of the individual techniques depends on five factors, namely: accuracy, forecast horizon, costs, complexity, and availability of data. This confirms the difficulty of executing such a task.

However, combined forecasts have been useful when it is difficult to select the most accurate forecasting method (Panagiotopoulos, 2011). In addition, combination methods that are considered more complex not always demonstrate a superior performance in comparison with simpler methods when the forecast accuracy is concerned (Stock & Watson, 2004; Konig et al., 2005).

For the choice of the best forecast, the analysis of a variety of techniques (including combinations) and the verification of which one provides the best accuracy value is necessary. For this purpose, one can use relative or absolute measures like MAPE, MSE or MAE. The incorporation of more than one comparison criteria is recommended, because different criteria may provide different indications of performance (Werner, 2012). 

Joint analysis and absence of concurrence among the forecasting techniques are supported in the literature, since the proposed use of a final model that adds the forecasting power of more than one technique is being successfully utilized in different areas of research (Clemen, 1989; Makridakis et al., 1998; Andrawis et al., 2011). In the following section, proposed configurations for the forecasting combinations are presented.

2.3 Methods for the combination of forecasts

The combination of forecasts can be represented by:

The main methods that estimate weights in forecast combinations are presented in the following paragraphs.

The simplest way to combine individual forecasts is to attribute the same weight to all of them. In this case, one has the arithmetic mean expressed by:

According to Werner (2012), the use of a forecast combination with arithmetic mean is a frequent practice because it is one of the most known procedures and is also easy to implement in a computer. This type of average, even considering that it does not use optimal weights, can generate better results than more sophisticated methods. Other similar forms of combination are the geometric mean (Faria & Mubwandarikwa, 2008), exposed in equation 6, and the harmonic mean (equation 7).

The approach in which the weights are values that are inversely proportional to a certain accuracy measure was the most sophisticated method presented by Bates and Granger (1969). For such a purpose, forecasts are combined in agreement with their individual performances. More specifically, the weight of a forecast is established in agreement with the accuracy of a forecast that was compared to the sum of all the forecasting methods. For Newbold and Bos (1994), methods that have a better performance must have a larger weight. In this context, weights are assumed to be inversely proportional to the accuracy indexes of all the techniques. The formulas that calculate the weights based on values that are inversely proportional to MAPE (equation 8), MSE (equation 9), and MAE (equation 10), are:

In which the parameters are  as previously defined.

In the approach using the minimum variance method, the weights are calculated in a way that the errors variance of the combined forecasts may be minimized.  In this case, the two forecasts are assumed to be unbiased, with variances in the forecasting errors equal to  and , respectively (Werner & Ribeiro, 2006; Andrawis et al., 2011; Martins & Werner, 2012). Assuming  as the covariance of the forecasting errors from the two techniques and assuming that the weights add up to 1, it is demonstrated that the optimal weights are:

According to Andrawis et al. (2011) the error in the covariance estimation can have an impact on the combined forecasting accuracy. In this context, if the covariance between the two forecasting techniques is not considered, the optimal weights are determined by:

Stock and Watson (1999) introduced a method in which the weights are calculated using a function that is inversely proportional to the mean square error (MSE). This approach is very similar to the uncorrelated minimum variance method presented in equation (9).

When there are two forecasts to be combined, an alternative approach to calculate the weights for the forecast is to use a multiple linear regression model. According to Granger and Ramanathan (1984), the conventional methods that combine forecasts can be viewed as a structured form of regression. These authors state that the methods are equivalent to the ordinary least squares method (OLS), where the combined forecast is the response variable and the individual forecasts are the explanatory variables.

In accordance with Andrawis et al. (2011) there are typically three equations that change according to the flexibility that is allowed in the weights, being given by:

The model expressed by equation (13) contains an intercept that can be useful to correct any possibility of biasing. Secondly, in the models presented in equations (14) and (15) it is assumed that the individual models are not biased.

In the method that uses weights based on absolute errors, weights are designated in accordance with the number of times that the corresponding technique is better than the others, and it is executed by means of a relative frequency (Panagiotopoulos, 2011). One way to use this approach for two techniques is to use a logical function, like:

The determination of weights in combinations can be also obtained by using linear programming (LP). Relevant studies dealing with this theme are detailed in Panagiotopoulos (2011). The mathematical elaboration of weight estimation to minimize MAPE, MSE and MAE, restricted to weights in a "0" to "1" interval, can be represented by:

In equations (17), (18), and (19), the number of observations in the procedure is given by .

There are many methods available in the literature to determine the weights used in the combination of individual forecasting techniques. However, the 14 methods presented in equations (5) to (19) have been elected as the most adequate in the scope of this study. Based on the above considerations, the next section presents the methodology proposed for this research.

3. Methodology

Because there is no unanimous agreement among authors about which individual techniques or methods of combination provide more accurate forecasts, in this study three individual techniques (ES, SARIMA, and ANN) and fourteen combination methods are tested, both with objective approaches. By doing this, the reduction of risks and confidence intervals are established as goals, as well as the use of independent information from different methods.

Major problems, given their magnitude and complexity, must be broken down into phases and stages. This is the case during the implementation of a forecasting system for the volume of withdrawn water in cities aiming at the public supply. In this context, this study used the propositions by Makridakis et al. (1998) as a basis, recommending the following steps: (i) problem definition; (ii) information gathering; (iii) preliminary analysis of data; (iv) model choice and fitting; and (v) use and evaluation of the forecasting model. Recommendations elaborated by Armstrong (2001) have also been taken into account in this work, namely: (i) structuring of the problem; (ii) information gathering; (iii) forecasting methods implementation; and (iv) evaluation of the forecasting method.

As stated by Martins and Werner (2012), it is necessary to analyze a wide range of techniques and verify which one provides the best measures of accuracy for obtaining a better forecast. For this purpose, in this study three individual techniques and fourteen combination methods have been comparatively evaluated using three different measures of accuracy.

The proposed systematic approach can be presented as a sequence of steps for its implementation. Figure 1 presents a flowchart of the systematic approach aimed at solving the problem of identifying individual forecasting techniques and their combinations, which are capable of adequately forecasting the water withdrawal volume for urban use.

Source: Elaborated by the authors.
Figure 1 – Structuring Phases of the Proposed Systematic Approach

The proposed systematic approach is divided into four phases. Each phase is divided into stages, with a noticeable characteristic of phase 3 that is divided into 14 stages. In the preliminary phase, named Phase 1, data from a data base were gathered, filtered, organized in time and validated such that one is capable of obtaining forecasting models. The individual techniques, the combination methods, and the adopted measures of accuracy for the identification of model(s) that are adequate for each series of study are defined in this phase. Finally, still in this phase, the series characteristics are identified as tendency and seasonality, and the sample is divided into fitting and validation.

Phase 2 consists in the identification of the models for the exponential smoothing (ES), SARIMA and ANN. Akaike information criterion (AIC) was chosen among the different criteria that can be used for the process of selecting the best model in ES and ARIMA classes (Snyder & Order, 2009). On the other hand, RMSE was used to select  the best ANN structure. Phase 2 is divided into three stages:

1. Stage 1:  Verify if there is an amount of historical data which is sufficient for the statistical modeling using Box-Jenkins methodology (a minimum of 50 observations). If there is,  an adequate SARIMA(p,d,q)x(P,D,Q)s model should be identified by AIC  to estimate its parameters and verify if the residues are uncorrelated. In addition, calculate the individual forecast, use the produced forecast to verify the accuracy, and keep the results to apply them in the model combinations. In steps: (i) identify an adequate SARIMA (p,d,q)x(P,D,Q)s model by AIC; (ii) estimate the model's parameters; (iii) check the adequacy of the model's coefficients: one verifies if the estimated coefficients are significant, that is, if they are related to a P-value < 5%; (iv) apply the randomness residues tests: ACF, PACF and the Box-Pierce test; and (v) obtain an individual forecast for the model and use the results to generate the combinations. In item (v), one seeks to identify if the residues form a sequence of random numbers. This sequence is frequently denominated as "white noise", since it has equal components at several frequencies. In the event that there is not enough data, only the models of ES (Stage 2) are used.
2. Stage 2: The identification of the best model of ES and the generation of forecasts are used in forecast combinations, analogously to the case of the SARIMA model. In this phase, only individual forecasts are employed.
3. Stage 3: Identifying the appropriate setting for ANN. For network learning, supervised training is used, whose synaptic weights are determined by minimizing RMSE. In order to obtain a proper setting for ANN, these parameters are calibrated by means of adaptation and scaling the application developed in MATLABÒ, by Dias (2006). This application allows you to set a small value for the parameters and gradually increase them at each network training. The value that has  the lowest RMSE in the training will be fixed as the definitive value of the parameter. Initially, the eg parameter is calibrated. The value of eg is fixed and then the sc parameter is calibrated. With eg and sc fixed, the value of the n parameter is calibrated. In this way there is a configuration suitable for ANN to perform the time series forecasting. The ranges explored in this work are: (i) the number of occurrences in the input vector: ranging from 2 to 12, increasing at the rate of 1 unit per training; and (ii) the value of eg and sc: variation of 0.01 to 4, increasing at the rate of 0.02 unit per training.

In phase 3, forecast combination methods are executed with values obtained from phase 2, assigning values to the individual techniques. The first stage of this phase consists in the choice of the best Statistical Technique - TS (ES or SARIMA) by the lowest MAPE in validation for each time series. The combination methods have been detailed in section 2.2, and the weights are: (i) of equal values (arithmetic mean, Comb AM); (ii) geometric mean (Comb GM); (iii) harmonic mean (Comb HM); (iv) inversely proportional to MAPE (Comb iMAPE); (v) inversely proportional to the MSE (Comb iMSE); (vi) inversely proportional to the MAE (Comb iMAE); (vii) obtained through the minimum variance with correlation (Comb MV ); (viii) obtained through the uncorrelated minimum variance (Comb MV ); (ix) obtained through Linear Programming (LP) and by minimizing the MAPE (Comb LP min MAPE); (x) obtained through LP while minimizing the MSE (Comb LP min MSE); (xi) through LP and by minimizing the MAE (Comb LP min MAE); (xii) through the Multiple Linear Regression Model (Comb MLRM) with intercept (Comb MLRM1); (xiii) through MLRM without intercept (Comb MLRM2); and (xiv) through relative frequency with the least absolute error (Comb FR).

Finally, in phase 4, by using the last part of the historical data reserved for the validation of each modeling (technique or method), the most accurate model is elected by the performance measures using MAPE, MSE, or MAE. It is thus possible to evaluate if the forecasts performed by individual models or their combinations can generate better forecasts.

The use of an evaluation criterion is necessary to assess and compare the predictive performance of the forecasts. In this study, three types of quantitative evaluation criteria are used: MAPE, MSE and MAE. It is highlighted that each one of those accuracy statistics is based on errors of one-step-ahead forecasting, that are calculated by the difference between the observed value at time  and the predicted value at time  for the time .

The structuring of the data base was performed in MS-Excelâ software. The software Statgraphics Centurion XVIâ was used for the modeling of data using ES and SARIMA techniques, as well as the statistical analysis of parameters significance of the model, the fitting quality (residues analysis) and the accuracy measures. On the other hand, the application developed in software MATLABÒ identified a suitable structure for ANN.

MS-Excelâ scripts to process spreadsheets have been written for the several proposed combinations. In addition, resources are used, like (i) the function "solver" to estimate weights with LP for the three methods (min MAPE, min MSE, and min MAE); and (ii) the function "data analysis – regression" to estimate the weights using MLRM for both methods (with and without intercept). The next section presents the main results from the systematic application of the proposal to evaluate the performance of forecasting models to forecast the water volume withdrawal in Brazilian cities.

4. Results and Discussion

To validate the  proposed systematic approach, the methodology was applied  to data from 10 cities of the micro-region in the  surrounding area of Pato Branco  City, in the southwestern region of the state of  Paraná, Brazil. The main phases used to analyze each time series following the proposed methodology were: (i) data gathering in the Water and Sanitation Company of Paraná (SANEPAR), which guaranteed the reliability of the information; (ii) representation and analysis of each series being studied by inspecting the behavior of the time plot, confirming that the most prominent features of the series are the growing and seasonal tendency over the years; (iii) individual modeling by ES, SARIMA and ANN; (iv) execution of the forecast combinations; and (v) comparison of the accuracy measures. Thus, three measures of performance for each modeling were calculated, resulting in a total of 510 (3x10x17).

In sections 4.1 to 4.3, the phases of the described systematic approach of section 3 are presented. As each city has its own distinguishing characteristics,  forecasts were directed to each one in an individual manner, trying to represent the tendency behavior and seasonality by means of ES, SARIMA, ANN models, and the 14 resulting combinations. In section 4.4, the obtained results are discussed.

4.1 Phase 1:  Time series and preliminary treatment

During the last decade it  has been observed from the available data for the 10 cities that they had different human development, social, and economic indexes, as shown in Table 2. The results of City's Human Development Index (CHDI) and FIRJAN Index of City Development (FICD) show the differences among the studied cities. For the analyzed period, the annual geometric population growth rate was positive for 5 cities and negative for the others. Population in those cities varied from 3,293 to 72,370 inhabitants. The annual growth rate of water connections was positive for the cities studied in the research.

Table 2: Characterization of the cities in the case-study

Source: Elaborated by the authors with data from PNUD1(2013), FIRJAN2(2012), IBGE3(2014) and IPARDES4(2012).

Data over a time span used in this study were gathered in the regional unit of the Water and Sanitation Company of Paraná (SANEPAR). These data were chosen because they represent real values being organized in time without the need to fill the data set with intermediate values. Figure 2 presents the historical behavior of the withdrawn volume of water in each of the 10 cities that served as a basis for the application of the proposed systematic approach. As expected, the plot of the historical time series revealed a growing tendency and a seasonal behavior in the withdrawn volume of water, even  though the population growth rate was negative in some cities.

For each one of the series, both the historical average level and the linear fit were determined with a regression equation (y = ax + b) and the coefficient of determination (R2). Moreover, the figure places an emphasis on the division that occurred in the gathered samples, which had one data set to fit each one of the models in the scope of this work and another data set for its validation. The sample is formed by 132 observations (11 years) for each city and was divided into a fitting set (120 observations) and a validation set (12 observations). Following the classification used by Makridakis and Hibon (2000), the planning horizon is a medium-term one if it is a 4- to 12- step-ahead forecast, and in this case, its horizon is twelve months.

Source: Elaborated by the authors.
Figure 2: Characterization of  time series of the withdrawn volume in the cities

4.2 Phases 2 and 3: Application of individual techniques and combination methods

The statistical software determined the most accurate models of ES and SARIMA classes. The individual models used in the combinations were considered adequate, as they have significant coefficients and uncorrelated residues. The best statistical technique (ST) was chosen for each city by MAPE validation. The computational application developed in MATLABÒ has identified an adequate structure of ANN for each assessed city. The main results for the individual techniques are shown in Table 3, which highlights the model identified for each series. The programmed spreadsheet permitted the proposed combinations. Moreover, Table 4 presents a summary of results obtained with the combinations. These results are analyzed in Section 4.4.

An initial observation is that all the individual techniques presented appropriate values for MAPE (below 8%). The methods of ES (ESHW or SES, the only ones elected from among those evaluated) were more accurate by 50% in the evaluated cities and also obtained the best average performance measured by validation MAPE. Moreover, SARIMA models also showed promising results, being superior in 40% of the cities, with average MAPE close to those obtained by ES models. On the other hand, the technical artificial intelligence (TAI) evaluated (ANN) was more accurate only in one city and had the lowest average MAPE in comparison with TS.

Table 3:  Identified Model and comparison of the performance of individual techniques by MAPE

  Source: Elaborated by the authors.

Table 4:  Performance of combinations measured by MAPE

Source: Elaborated by the authors.

4.3 Phase 4: Determination of the most accurate model

To evaluate the individual performance of each model (technique or method), the best model for each city was identified according to adopted accuracy measures. These results are listed in Table 5. For the proposed phase application, the MAPE criterion of the validation was adopted to choose the most accurate individual model.

Table 5:  Determination of the most accurate model by the MAPE criterion

Source: Elaborated by the authors.

According to MAPE criterion, the individual techniques were improved in 40% of the evaluated cases. On the other hand, the used combination methods were more accurate in 60% of the evaluated series. In addition, while the individual techniques  produced MAPE lower than 8%, the choice of the most accurate model (individual technique or combination method) generated MAPE lower than 4% for the 10 cities.

Section 4.3 ends the application cycle of the proposed methodology and shows how the managing board in a water and sanitation company could have forecasts for a whole micro-region in a horizon of 12 months with MAPE below 4%. The managing plan for the following year could be done based on a good forecasting system. However, the company should have an automatized systematic approach or there should be an employee with specific knowledge to select the presented individual models. One question may be raised in this situation: "Would there be a technique/method that could be directly applied and produce good results?" In the following section this discussion is developed on an analysis using the descriptive statistics of the results.

4.4 Discussion: Statistical Analysis of the Results

The first analysis consisted in the identification of frequency with which each model was selected (individual technique or combination method) for the 10 cities and each measure of accuracy: MAPE, MSE and MAE, as demonstrated in Table 6.

Table 6:  Summary of the performance in the validation of the evaluated techniques and methods of combination

 Source: Elaborated by the authors. 

The ES technique and the combination methods "iMAPE, iMSE, iMAE, MV (r ¹ 0) and MV (r = 0)" reached a better average performance in the validation, measured by MAPE, and an even higher  average performance was obtained in the fitting. On the other hand, validation had a worse performance than the fitting for SARIMA and ANN models, Comb AM, Comb GM, Comb HM, Comb LP min MAPE, Comb LP min MSE, Comb LP min MAE, Comb MLRM1 and Comb MLRM2.

Next, an analysis of absolute and relative frequencies was performed. It had the goal of identifying the frequency at which a forecast model presents the best average accuracy evaluated for each series of data. Because of that, the number of times that each model obtained the best performance in the validation was registered, with the results listed in Table 7.

Table 7: Frequencies and percentages of the measures of accuracy in the validation

 Source: Elaborated by the authors.

On the comparison between the individual techniques and the combination methods for the three criteria (MAPE, MSE and MAE), we highlight: (i) individual techniques (10/30), with SARIMA (3/30) and ES (7/30); (ii) combinations (20/30); (iii) best TS (ES or SARIMA): 10/30 (33.33%); (iv) Comb MRLM1, that is, with intercept: 9/30 (30%); (v) Comb VM (p = 0): 4/30 (13.33%); (iv) Comb AM: 2/30 (6.67%); (vii) Comb iMAPE: 2/30 (6.67%); (viii) Comb FR: 2/30 (6.67%); (ix) Comb iMSE: 1/30 (3.33%); and (x) other evaluated techniques (ANN and other combinations) that were  less accurate in all the cases studied.

These results show that the combinations of methods generally produce more accurate results than the individual predictors. The unanimity occurred in the cities of Pato Branco (SARIMA), Palmas and Mangueirinha (Comb MLRM1), Coronel Vivida (Comb MV p = 0) and Bom Sucesso do Sul (SE). In this case, the Comb MLRM1 method (30%) was the one which presented a better performance concerning the other three adopted accuracy measures. It is also worth mentioning the ES technique (23.33%).

We observe that MAPE and MSE accuracy measurements always pointed to the same model. For MAPE and MAE criteria, the individual techniques were superior in 40% of the evaluated cases with a proportion of 30% to ES and 10% to SARIMA, while the combinations were better in 60% of the experiments. On the other hand, for MSE criterion, ES and SARIMA techniques were each successful in 10%  of the cases, while the combinations were more accurate in 80% of the cases.

As previously cited, three measures were used in this study as a way to measure the performance of the evaluated modeling. However, the frequency at which a model is chosen does not demonstrate unanimity. For comparison of different series, MAPE statistics are the most used measures of accuracy (Kahn, 1998; Goodwin & Lowton, 1999).

In this sense, Table 8 presents the results for this specific measure of accuracy. However, analysis will be performed for all the referred measures. This table presents descriptive statistics of central tendency and dispersion of validation MAPE to evaluate the global performance observed in the modeling and combinations. The best results for each statistics are in bold format.

Observing results shown in Table 8 allows one to infer that the "worst" model (technique or method) had a maximum MAPE which was less than 8% (7.853%) in all the 10 studied cities. For all the evaluated measures of accuracy, the best TS had a superior average performance for the evaluated statistics: minimum, maximum, mean and standard deviation.

Table 8:  Performance measured by the validation MAPE

  Source: Elaborated by the authors.

According to MAPE, the best TS (ES or SARIMA) had the best average performance, the lowest standard deviation, and the lowest coefficient of variation. In the evaluated context, it is the best forecast technique, even when compared to fourteen combinations. On the other hand, TAI (ANN) had the worst performance  in average and standard deviation categories.  Globally, the combinations did not generate, on average, more accurate results than the best TS (ES or SARIMA). However, the average performance of individual methods was inferior to those of combinations as shown in Table 8.

The best general performance in those 10 cities, according to MAPE variability criterion in validation, was the most accurate model of ES or SARIMA, for it obtained the lowest standard deviation (0.667%) and had the best average (3.297%). On the other hand, the ANN model obtained MAPE maximum average (4.509%) and the highest average standard deviation (1.366%), resulting in the "worst" evaluated model in this aspect.

The best techniques measured by average MSE of the validation had the most accurate model of ES or SARIMA with the best performance (Appendix I). Similar results occurred for MAE accuracy measurements of validation statistics concerning mean and standard deviation (Appendix I).

Comparing the results obtained in this study with data from M-Competition (Makridakis & Hibon, 2000; Nikolopoulos et al., 2011b), the following points have been ratified:

1. The most sophisticated methods do not necessarily have better performance than simpler ones. Moreover, in several situations the performance of simpler methods was also better. In Table 8, it is observed that the individual ES technique overcame the SARIMA technique and ANN;
2. It was confirmed that the choice of an accuracy verification criterion in the forecasting has an influence on the method classification regarding the forecasting performance. In the analysis provided in this work, when the set of time series is available, this work could use the variability of validation MAPE as a criterion.

Small  differences, which were evaluated by the position and dispersion statistics among the evaluated models, indicate that any one of such models is likely to be used. However, given the simplicity, performance, and significance, the ES model is the most usable to estimate the withdrawn volume. This is particularly important if there is no specialist or dedicated professional in the water company with knowledge to apply the proposed methodology.

5. Conclusions

The operationalization of strategies in the planning and control of production is strongly attached to the existence of an effective forecasting system. Sustainability, knowledge of historical reality, and future developments must be known for planning to be implemented. In this sense, all the evaluated models generate adequate forecasts in agreement with the scope of this work. Furthermore, the combinations increase the accuracy of forecasts.

As a first result, predominant frequency is highlighted, at which the results from combinations of forecasts presented a superior performance when compared with individual forecasts, since the worst MAPE was reduced from 8% to 4%.

From the results, it was demonstrated that all the evaluated models are adequate to produce a volume withdrawal forecast in a horizon of 12 months. The best general performance in the 10 cities, according to the variability of validation MAPE, was the best statistical techniques (MAPE = 3.297% and standard deviation = 0.667%). In the practical context of water and sanitation companies, they can use an individual ES technique with good results (MAPE = 3.569% and standard deviation = 0.758%).

The classification and choice of the best method depends on the forecast horizon; that is, a few methods have a better performance for a determined horizon and are worse for others (Makridakis & Hibon, 2000). Thus, the obtained results should not be used for planning horizons which are not contemplated in this study. However, a systematic approach can be used without significant changes in it.

The research conclusion allows one to list some future works like: (i) the application of the proposed systematic approach in a set with more cities, in different regions and with larger populations; and (ii) the performance evaluation for a wider forecast horizon from 2 to 5 years.

Acknowledgements

Authors thank the Water and Sanitation Company of Paraná (SANEPAR) for providing data and support in the development of this research.

References

[1] Adamowski J. 2008. Peak Daily Water Demand Forecast Modeling Using Artificial Neural Networks. Journal of Water Resources Planning and Management, 134:119-128.

[2] Almutaz I, Ajbar A,  Khalid Y & Ali E. 2012.  A probabilistic forecast of water demand for a tourist and desalination dependent city: Case of Mecca, Saudi Arabia. Desalination, 294:53–59.

[3] ANA. Agência Nacional de Águas. Conjuntura dos recursos hídricos no Brasil. Special Edition. Brazil. Available at http://arquivos.ana.gov.br/imprensa/arquivos/Conjuntura2012.pdf. Accessed at 24 nov. 2012.

[4] Andrawis RR, Atiya AF & El-Shishiny H. 2011. Combination of long term and short term forecasts, with application to tourism demand forecasting. International Journal of Forecasting, 27(3):870-886.

[5] Armstrong JS. 2001. Principles of Forecasting: A handbook for researchers and practitioners. New York: Kluwer Academic Publishers.

[6] Azadeh A, Neshat N & Hamidipour H. 2012. Hybrid Fuzzy Regression – Artificial Neural Network for Improvement of Short-Term Water Consumption Estimation and Forecasting in Uncertain and Complex Environments: Case of a Large Metropolitan City. Journal of Water Resources Planning and Management, 138:71-75.

[7] Bates JM & Granger CWJ. 1969. The Combining of Forecasts. Operational Research Quarterly, 20(4):451-468.

[8] Caiado J. 2010. Performance of combined double seasonal univariate time series models for forecasting water demand. Journal of Hydrologic Engineering, 15(3): 215-222.

[9] Calôba GM, Calôba LP & Saliby E. 2002. Cooperação entre redes neurais artificiais e técnicas clássicas para previsão de demanda de uma série de vendas de cerveja na Austrália. Pesquisa Operacional, v.22, n.3, p.345-358.

[10] Campisi-Pinto S, Adamowski J & Oron G. 2012. Forecasting UrbanWater Demand Via Wavelet-Denoising and Neural Network Models. Case Study: City of Syracuse, Italy. Water Resour Management, 26:3539–3558.

[11] Clemen RT. 1989. Combining forecasts: A review and annotated bibliography. International Journal of Forecasting, 5:559–583.

[12] Dias G J C. 2006. Planejamento estratégico de um centro de distribuição: uma aplicação de redes neurais artificiais de funções de bases radiais para previsão de séries temporais. Dissertação de mestrado. Universidade Federal do Paraná.

[13] Faria F & Mubwandarikwa, E. 2008. The Geometric Combination of Bayesian Forecasting, Journal of Forecasting, 27/6:519-535.

[14] Firat M, Turan ME & Yurdusev MA. 2010. Comparative analysis of neural network techniques for predicting water consumption time series. Journal of Hydrology, 384:46–51.

[15] FIRJAN. Federação das Indústrias do Rio de Janeiro. Índice Firjan de Desenvolvimento Municipal. Available at http://www.firjan.org.br/ifdm/downloads/. Accessed at 12 dec. 2012.

[16] Flores BE & White EM. 1988. A Framework for the Combination of Forecasts. Journal of the Academy of Marketing Science, 16(3-4):95-103.

[17] Ghiassi M, Zimbra D & Saidane H. 2008. Urban Water Demand Forecasting with a Dynamic Artificial Neural Network Model. Journal of Water Resources Planning and Management, 134:138-146.

[18] Goodwin P & Lawton R. 1999. On the asymmetry of the symmetric MAPE. International Journal of Forecasting, 15:405-408.

[19] Granger CWJ & Ramanathan R. 1984. Improved Methods of Forecasting. Journal of Forecasting, 3:197-204.

[20] Herrera M, Torgo L, Izquierdo J & Pérez-García R. 2010. Predictive models for forecasting hourly urban water demand. Journal of Hydrology, 387:141–150.

[21] Hibon M & Evgeniou T. 2005. To combine or not to combine: selecting among the forecasts and their combinations. International Journal of Forecasting, 21(1):15-24.

[22] IBGE. Instituto Brasileiro de Geografia e Estatística. IBGE Cidades @. Available in http://www.ibge.gov.br/cidadesat/topwindow.htm?1. Accessed at 12 dec. 2014.

[23] IPARDES. Instituto Paranaense de Desenvolvimento Econômico e Social. Caderno Estatístico por Município. Available at http://www.ipardes.gov.br/cadernos/. Accessed at 12 dec. 2012.

[24] Kahn KB. 1998. Benchmarking sales forecasting performance measures. The Journal of Business Forecasting Methods & Systems, 17(4):19-23.

[25] Konig AJ, Franses PH, Hibon M & Stekler HO. 2005. The M3-Competition: Statistical tests of the results. International Journal of Forecasting, 21: 397–409.

[26] Makridakis S & Hibon M. 2000. The M3-Competition: results, conclusions and implications. International Journal of Forecasting, 16: 451-476.

[27] Makridakis S, Wheelwright SC & Hyndman RJ. 1998. Forecasting, Methods and Applications. Third Edition. New York: John Wiley & Sons.

[28] Martins V & Werner L. 2012. Forecast combination in industrial series: A comparison between individual forecasts and its combinations with and without correlated errors. Expert Systems with Applications, 39:11479–11486.

[29] MATLAB. 2010. Version 2010b. User's Guide.

[30] Mohamed MM & Al-Mualla AA. 2010. Water demand forecasting in Umm Al-Quwain using the constant rate model. Desalination, 259:161–168.

[31] Morettin P A & Toloi C M C. 2006. Análise de Series Temporais. 2. ed. São Paulo: Editora Edgard Blucher.

[32] Newbold P & Bos T. 1994. Introductory Business & Economic Forecasting. Second edition, South-Western Publishing, Cincinnati. Ohio.

[33] Nikolopoulos K, Makridakis S, Assimakopoulos V & Syntetos A. 2011a. The M4 Competition. The Operational Research Society Conference 2011 – OR53, September 6-8 2011, East Midlands Conference Centre, Nottingham, UK.

[34] Nikolopoulos K, Makridakis S, Assimakopoulos V & Syntetos A. 2011b. The M4 Competition: Facts, Innovations and the way Forward. 31th Annual International Symposium on Forecasting - ISF 2011, June 26-29. 2011, Prague, Check Republic.

[35] Panagiotopoulos A. 2011. Optimising Time Series Forecasts Through Linear Programming. University of Nottingham. Thesis PhD. 205p.

[36] PNUD. Programa das Nações Unidas para o Desenvolvimento. 2010. Ranking do IDH dos Municípios do Brasil. Available at http://atlasbrasil.org.br/2013/ranking. Accessed at 17 may 2014.

[37] Qi C & Chang N. 2011. System dynamics modeling for municipal water demand estimation in an urban region under uncertain economic impacts. Journal of Environmental Management, 92:1628-1641.

[38] Snyder RD & Ord J K. 2009. Exponential Smoothing and the Akaike Information Criterion, Monash Econometrics and Business Statistics. Working Papers, 4/09 working paper, Monash University.

[39] Souza GP. 2008. Método para estruturar a integração de previsões utilizando a técnica Delphi. Tese de doutorado. Universidade Federal de Santa Catarina.

[40] Statgraphics Centurion XVI. 2014. Version 16.1.17. User's Guide.

[41] Stock JH & Watson MW. 1999. A comparison of linear and nonlinear univariate models for forecasting macroeconomic time series. In Engle RF, White H (eds). Cointegration, Causality and Forecasting: A Festschrift for Clive W.J. Granger. Oxford University Press: Oxford; 1–44.

[42] Stock JH & Watson MW. 2004. Combination forecasts of output growth in a seven-country data set. Journal of Forecasting, 23: 405–430.

[43] VWNA. Veolia Water North America. Finding the blue path for a sustainable economy. Chicago, EUA. Available in http:// www.veoliawaterna.com. Accessed at 24 nov. 2012.

[44] Werner L & Ribeiro JLD. 2006. Modelo composto para prever demanda através da integração de previsões. Revista Produção, 16(3):493-509.

[45] Werner L. 2012. Na busca pelo melhor método de combinação de previsões por meio do uso de análise de correlação e análise de cluster. In: II Congresso Brasileiro de Engenharia de Produção. Ponta Grossa, PR, Brasil, 28 a 30 de novembro de 2012.

[46] Wu Y, Liu B, Li M, Tang k & Wub Y. 2013. Prediction of CO2 solubility in polymers by radial basis function artificial neural network based on chaotic self-adaptive particle swarm optimization and fuzzy clustering method. Chinese Journal of Chemistry, 31: 1564-1572.

[47] Yasar A, Bilgili M & Simsek E. 2012. Water Demand Forecasting Based on Stepwise Multiple Nonlinear Regression Analysis. Arabian Journal for Science and Engineering, 37:2333–2341.

[48] Zhai Y, Wang J, Teng Y & Zuo R. 2012. Water demand forecasting of Beijing using the Time Series Forecasting Method. J. Geogr. Sci, 22(5): 919-932.

APPENDIX I - COMPLEMENTARY RESULTS

Table 1:  Performance measured by the validation MSE

Source: Elaborated by the authors. Note: Best mean and standard deviation
for the best TS (ES or SARIMA), according to the MSE criterion.

Table 2: Performance measured by the validation MSE

Source: Elaborated by the authors. Note: best mean, standard deviation, maximum
and minimum for the best TS (ES or SARIMA), according to the MAE criterion.

1. Graduate Program in Industrial & Systems Engineering (PPGEPS/UTFPR-PB). Federal Technological University of Paraná (UTFPR).  Pato Branco Campus-PR. Email: donizetti@utfpr.edu.br
2. Graduate Program in Industrial & Systems Engineering (PPGEPS/UTFPR-PB). Federal Technological University of Paraná (UTFPR).  Pato Branco Campus-PR. Email: gilson@utfpr.edu.br  
3. Graduate Program in Industrial & Systems Engineering (PPGEPS/UTFPR-PB). Federal Technological University of Paraná (UTFPR).  Pato Branco Campus-PR. Email: marcelo@utfpr.edu.br  

4. Graduate Program in Industrial & Systems Engineering (PPGEPS/UTFPR-PB). Federal Technological University of Paraná (UTFPR).  Pato Branco Campus-PR. Email: batistus@utfpr.edu.br  

5. Regional Unit - Pato Branco. Sanitation Company of Paraná (SANEPAR). Pato Branco-PR. Email:bayerl_quimica@hotmail.com