Water Consumption Forecasting Using ARIMA and Holt-Winters Methods: A Case Study in Vouzela, Portugal

1. Introduction

The efficient management of water resources is a growing challenge worldwide, especially in the face of climate change and increased demand for water. In this context, accurate forecasting of water consumption becomes a crucial tool for planning and decision-making by water managers. Water consumption in residential areas is influenced by several factors, including seasonal variations, household usage patterns or demographic trends. The complexity of these factors makes forecasting water consumption a significant challenge for local councils as water managers. Advanced statistical methods, such as ARIMA (Autoregressive Integrated Moving Average) or Holt-Winters, have proven effective in forecasting time series in various fields, including water consumption, as they are particularly useful for capturing complex patterns in historical data and projecting future trends. This study focuses on the application of the ARIMA and Holt-Winters methods, using python technology, to forecast monthly water consumption in a specific area of Cambra, in the municipality of Vouzela, as illustrated in Figure 1. The choice of this locality as a case study offers a unique opportunity to examine consumption patterns in a medium-sized Portuguese community, potentially providing knowledge applicable to other nearby regions with similar characteristics. The main objective of this work is to develop and evaluate robust forecasting models that assist in water supply management and planning. By using five years worth of historical data (2018-2022), in addition to forecasting consumption for the following year, a better understanding of the underlying patterns of water consumption in the region is obtained. In addition to the application of statistical methods, this study addresses common challenges in time series analysis, such as data normalisation, stationarity checking and the selection of optimal parameters, specifically for the ARIMA method. The methodology applied aims to provide a replicable approach that can be adapted and applied in other water resource management contexts. The following sections present a review of the relevant literature, details of the methodologies employed for forecasting the data through experimental work, a discussion of the results obtained and a brief approach to the implications of this study for water resources management. Finally, directions for future work will be indicated, including the recommendation to implement telemetry with LoRaWAN technology [13]. The rest of the article is organised as follows. Section 2 presents the related work. Section 3 describes the experimental work. Section 4 deals with the results observed during the application of the methods and, finally, the conclusions and future work are presented in Section 5.

2. Related Work

The ARIMA (Autoregressive Integrated Moving Average) and Holt-Winters methods have been widely used in various fields to forecast time series, demonstrating their versatility and effectiveness. This section presents an overview of some studies that have applied these statistical methods in different contexts, providing a comparative basis for the present work. In the field of energy, Ediger and Akar (2007) applied ARIMA methods and SARIMA (Seasonal AutoRegressive Integrated Moving Average), an extension of ARIMA that specifically deals with seasonal components in time series, to forecast primary energy consumption in Turkey [14]. The results of the forecast indicated that the growth rate of energy needs will decrease, with the exception of wood and vegetable waste, emphasising the relevance of renewable sources to guarantee energy security. In the electricity sector, Contreras et al. (2003) used ARIMA methods to forecast electricity prices in competitive markets. Tests carried out in the Spanish and California markets revealed average daily errors of 5-10%, demonstrating the applicability of ARIMA in volatile markets [15]. For natural gas consumption, Taspinar et al. (2013) compared ARIMA methods with approaches based on Fourier series expansion. Although ARIMA was outperformed in this case, the study emphasises the importance of considering seasonal variations and temperature feedback in energy consumption forecasts [16]. The COVID-19 pandemic has prompted several studies using ARIMA to forecast cases. Chakraborty and Ghosh (2020) applied the method to forecast the number of cases in five severely affected countries. Similarly, Zeroual et al. (2020) compared ARIMA with Machine Learning (NARNN - Nonlinear Autoregression Neural Network) and (LSTM - Long-short term memory) approaches to predict cases in European countries [17,18]. Maleki et al. (2020) proposed an algorithm to calculate the best ARIMA parameters per country, analysing 145 countries in six geographical regions. This study highlighted the relationship between the behaviour of the virus and population variables, suggesting the incorporation of factors such as humidity and culture in future models [19]. In the financial field, Mondal et al. (2014) evaluated the effectiveness of ARIMA in forecasting Indian stock prices, achieving an accuracy of over 85 per cent in various sectors. This study demonstrates the applicability of ARIMA in financial markets, although it emphasises the need for more specific methods for certain sectors [20]. In the field of aviation, Alhindi et al. (2021) used ARIMA to predict visibility based on meteorological data, a crucial factor for air safety. The promising results suggest exploring Machine Learning techniques to further improve forecast accuracy [21]. In the context of agricultural price forecasting, the study applied the ARIMA method to forecast potato prices in Turkey, demonstrating that ARIMA outperformed other methods, such as Holt-Winters, in accuracy. The results indicate that ARIMA is effective in modelling time series, in this case prices, highlighting its relevance for formulating economic policies. This work contributes to the literature by highlighting the importance of price stability in agriculture [22]. Still in the area of time series forecasting, Koehler et al. (2001) explored new methods based on the multiplicative Holt-Winters method, which is fundamental for dealing with data that has a trend and seasonality. The research demonstrated that the forecast error variance can be adjusted to reflect these characteristics, improving forecast accuracy. The results indicate that the correct choice of model is crucial, especially in series with an upward trend, which can have significant implications in a number of practical applications [23]. These studies demonstrate the versatility and effectiveness of statistical methods in various fields, from energy and public health to finance and aviation. The successful application of ARIMA in forecasting energy consumption, as demonstrated by Ediger and Akar (2007), provides a solid basis for our approach to forecasting water consumption. Furthermore, the lessons learnt from the COVID-19 case study, particularly the importance of considering regional and seasonal factors, can be applied to our water consumption analysis. Maleki et al.’s (2020) methodology for optimising ARIMA parameters by region offers important insights for the localised approach in Cambra, Portugal. This work builds on these successful applications of ARIMA and Holt-Winters, adapting them to the specific context of water consumption forecasting, with the aim of contributing to more efficient management of water resources.

3. Experimental Work

The study is based on a set of water consumption data, in cubic meters (m³), from 2018 to 2022, from 41 meters in a specific area of Confulcos, in the parish of Cambra, in the municipality of Vouzela. The first step in processing the data was to build a dataset, as shown in Table 1, and then a time series.

Using a list previously generated with the monthly data extracted from the dataset, and the indices generated for this purpose, we can see in Figure 2 the series in the form of a line graph, where we can see an upward trend and a pattern of higher consumption at certain periods of each year.

In Figure 3, after decomposing the series, we can clearly see the trend and seasonality, as well as the residuals. These residuals correspond to the part of the series that is not explained by the trend and seasonality. The analysis reveals that the series shows a seasonal pattern, with variations over time, and an upward trend, showing an increase in water consumption, particularly in the summer months.

These steps proved crucial to understanding the time series. The visual analysis allowed us to identify a clear trend in the series, indicating that the average does not oscillate around a constant value. It was also possible to detect a recurring pattern of consumption over the period analysed. With the aim of drawing up a forecast using the ARIMA method in mind, the next step was to assess whether the data had an approximately normal distribution. Although not strictly necessary for the application of the model, this check is advisable because it can help validate the suitability of the model, ensuring that the residuals are random and do not show systematic patterns, which is fundamental for the reliability of the forecasts. However, looking at the graph shown in Figure 4 raised some uncertainties about this normality.

One of the techniques for checking the normality of data is to perform the Shapiro-Wilk test, with a significance level of 5% (commonly used). This test is recommended for a set of small and moderate samples, but it maintains high precision even at higher numbers (n < 100), compared to other tests such as Kolmogorov-Smirnov with Lilliefors correction [24,25]. The Shapiro–Wilk (SW) test is a well-known method to assess whether a sample comes from a normally distributed population. It is based on the following statistic:

$$W={\displaystyle \frac{{\left({\sum }_{i=1}^n{a}_i{X}_{\left(i\right)}\right)}^2}{{\sum }_{i=1}^n{(X_i-\overline{X})}^2}}$$

(1)

where:

• W is the Shapiro–Wilk test statistic,
• n is the number of observations in the sample,
• $X_{\left(i\right)}$ are the ordered sample values such that $X_{\left(1\right)}\le X_{\left(2\right)}\le \dots \le X_{\left(n\right)}$,
• $a_i$ are constants derived from the means, variances, and covariances of the order statistics of a sample from the standard normal distribution,
• $X_i$ are the original (unordered) sample values,
• $\overline{X}$ is the mean of the sample.

The test provides a p-value, which is the probability of obtaining a result as extreme as the one observed, assuming that the null hypothesis (of normality) is true [26]. In this case, the p-value was found to be less than the significance level (0.05), which leads to the rejection of the null hypothesis and the acceptance of the alternative hypothesis, indicating that the distribution is not normal. In order to normalise the distribution of the data, we chose to apply a natural logarithmic transformation (base ‘e’). This choice was based on the fact that the series contains only positive integers, and this transformation is effective in reducing variance and improving normality. After applying the logarithmic transformation, the Shapiro-Wilk test was performed again. The new result ( 0.051) indicated that the distribution of the transformed data is now approximately normal, which is reflected in Figure 5.

In addition to normalising the data, applying the ARIMA method requires the series to be stationary. This means that its statistical properties, such as mean and variance, must remain constant over time, without exhibiting trends or seasonal patterns. Although the initial decomposition revealed an increasing trend and a periodic pattern in the series, the Kwiatkowski-Phillips-Schmidt-Shin (KPSS) test was carried out to confirm stationarity. The result of the test indicated that the series was not stationary, as the test statistic exceeded the critical value, leading to the rejection of the null hypothesis of stationarity. To make the series stationary, the differentiation technique was applied. This method creates a series by calculating the difference between consecutive values, subtracting each value from its successor [27]. This approach helps to eliminate trends and stabilise the average. Figure 6 shows the graph of the series after differentiation, where you can see that the trend has been removed and the average oscillates around a constant value, indicating that the series has become stationary.

After processing the series, ensuring that it met the necessary requirements for applying the ARIMA method, we proceeded to implement the method itself on the differentiated series.

3.1. ARIMA and SARIMA Method

Once the series had been processed, ensuring that it met the necessary requirements for the application of the ARIMA method, the method itself was implemented on the differentiated series. The ARIMA method has three parameters (p, d, q) and is often written as ARIMA(p, d, q), where p represents the order of the autoregressive component, d is the number of differentiations required to obtain a stationary ARMA(p, q) model, and q corresponds to the order of the moving average component. The general ARIMA(p, d, q) formula is expressed by a lag operator [28], represented as:

$$\varphi \left(B\right){(1-B)}^d{Y}_t=\theta \left(B\right)e_t$$

(2)

where the characteristic operators are:

$$\varphi \left(B\right)=1-{\varphi }_{1}B-{\varphi }_2{B}^2-\dots -{\varphi }_p{B}^p$$

(3)

$$\theta \left(B\right)=1-{\theta }_{1}B-{\theta }_2{B}^2-\dots -{\theta }_q{B}^q$$

(4)

and

$${(1-B)}^d{Y}_t={\nabla }^d{Y}_t$$

(5)

The definitions of the symbols are:

• $\varphi$ is the estimate of the autoregressive component parameter,
• $\theta$ is the estimate of the moving average component parameter,
• ${\nabla }^d$ is the differencing operator,
• B is the lag operator,
• $e_t$ is a purely random process with zero mean and variance ${\sigma }_e^2$ (white noise).

Crucial parameters were considered, such as the number of combinations for each parameter and the monthly seasonal period. We chose to use Auto-ARIMA, a method that automatically combines various regression orders (p), moving averages (q) and degrees of differentiation (d), with the aim of identifying the most suitable model. It is important to emphasise that Auto-ARIMA is an algorithm that can automatically select the best model generated by the ARIMA or SARIMA methods for a time series. When selecting a SARIMA method, it incorporates seasonality, resulting in a combination of the parameters of the standard method (p,d,q) with the seasonal parameters (P,D,Q)m, where m represents the seasonal period [29]. The SARIMA(p,d,q)(P,D,Q)m method represented as:

$${\Phi }_P\left(B^m\right)\varphi \left(B\right){(1-B)}^d{(1-B^m)}^D{Y}_t={\Theta }_Q\left(B^m\right)\theta \left(B\right)e_t$$

(6)

where the seasonal operators are:

$${\Phi }_P\left(B^m\right)=1-{\Phi }_1{B}^m-{\Phi }_2{B}^{2m}-\dots -{\Phi }_P{B}^{Pm}$$

(7)

$${\Theta }_Q\left(B^m\right)=1+{\Theta }_1{B}^m+{\Theta }_2{B}^{2m}+\dots +{\Theta }_Q{B}^{Qm}$$

(8)

and

$${(1-B^m)}^D{Y}_t={\nabla }_m^D{Y}_t$$

(9)

The definitions of the symbols are:

• $\varphi$ is the estimate of the non-seasonal autoregressive component parameter,
• $\Phi$ is the estimate of the seasonal autoregressive component parameter,
• $\theta$ is the estimate of the non-seasonal moving average component parameter,
• $\Theta$ is the estimate of the seasonal moving average component parameter,
• B is the lag operator,
• d and D are the orders of non-seasonal and seasonal differencing, respectively,
• m is the number of periods per season,
• $e_t$ is a purely random process with zero mean and variance ${\sigma }_e^2$ (white noise).

Since the method was applied to an already differentiated series, Auto-ARIMA identified the optimal model without the need for additional degrees of differentiation, as can be seen in Figure 7.

The report produced by the Auto-ARIMA method identified the most appropriate model using the Akaike Information Criterion (AIC). The AIC, is a statistical tool that assesses the relative quality of statistical models for a specific data set, balancing the accuracy of the fit with the complexity of the model [29,30]. The AIC represented as:

$$AIC=-2log\left(L\right)+2k$$

(10)

where:

• L is the maximized likelihood of the model,
• k is the number of estimated parameters in the model.

Figure 8 shows the detailed results of the SARIMA method chosen, which includes seasonality. This method integrates autoregressive and moving average elements, with an order (1,0,1) for the non-seasonal component and an order (3,0,0) for the seasonal component, considering a 12-month seasonal cycle. Specifically, for the non-seasonal component, we have an autoregressive order (p) of 1, a degree of differentiation (d) of 0, and a moving average order (q) of 1. For the seasonal component, the seasonal autoregressive order (P) is 3, the seasonal degree of differentiation (D) is 0, and the seasonal moving average order (Q) is 0. The seasonal period (s) is 12 months. The choice of this particular model suggests that it provides the best balance between fit to the data and complexity, as assessed by the AIC. This indicates that the model efficiently captures the fundamental characteristics of the time series, including trends, seasonal patterns and other fluctuations, while avoiding overfitting the data.

After modeling the series, the next step was to extract the residuals and analyze the autocorrelation between these data. In the context of time series, residuals are the differences between the observed values and the values predicted by the model. The analysis of these residuals is crucial, as it examines the relationship of the residuals within the same set, all associated with the water consumption variable. If the model is suitable, the residuals should be random and show no systematic patterns. It is important to note that the forecast is generated from these residuals, as they are derived from a series that has been properly treated and modeled. Analysis of the residuals made it possible to validate the model produced by the Auto-ARIMA method. This validation was possible because the model passed the normality test and, as mentioned above, the evaluation of autocorrelation and partial autocorrelation. The graphs presented illustrate that the residuals in both the Autocorrelation Function (Figure 9) and the Partial Autocorrelation Function (Figure 10) do not exhibit significant autocorrelation. This is evidenced by the fact that the bars, which represent periods in the time domain, remain within the 95% confidence interval.

This observation suggests that the model was able to efficiently capture the temporal structure of the data.

Once the series of residuals had passed the tests, we proceeded to forecast the data for the next 12 months, corresponding to the year 2023. It is important to note that the residuals are on a reduced scale, resulting from the transformations applied to the original series - first the logarithmic transformation and then the differentiation. Consequently, the forecast is also on this reduced scale, making it necessary to revert all the data to the original scale. The process of reversing the differentiation involves the reverse procedure to the one applied initially. Instead of subtracting consecutive values, each value of the differentiated forecast is added to the corresponding value of the non-differentiated series. This process is carried out iteratively, starting with the addition of the first value of the differentiated series to the last known value of the non-differentiated series. After reversing the differentiation, it is still necessary to apply the exponential function to the series, since the logarithmic function was initially applied to normalize the distribution, as indicated by the Shapiro-Wilk test. Figure 11 shows a graphical view of the consumption forecast for 2023, highlighted in red, after the data has been processed to return to the initial scale.

The numerical results of the forecast, obtained by the ARIMA method considered optimal according to the Akaike Information Criterion (AIC), are presented in tabular format together with the results of the Holt-Winters method for comparative purposes.

3.2. Holt-Winters

After forecasting water consumption for 2023 using the ARIMA method, the Holt-Winters method was applied to the original series, before the transformations carried out for the ARIMA method. Forecasting using the Holt-Winters additive method, also known as triple exponential smoothing, was chosen due to the observation of a constant seasonality in the time series. This method uses three fundamental equations to model the level, trend and seasonality of the series [28,31], represented as:

Level:

$$l_t=\alpha (y_t-s_{t-m})+(1-\alpha )(l_{t-1}+b_{t-1})$$

(11)

Trend:

$$b_t=\beta (l_t-l_{t-1})+(1-\beta )b_{t-1}$$

(12)

Seasonality:

$$s_t=\gamma (y_t-l_t)+(1-\gamma )s_{t-m}$$

(13)

Forecast:

$${\widehat{y}}_{t+h}=l_t+h{b}_t+s_{t+h-m(k+1)}$$

(14)

The definitions of the symbols are:

• $l_t$ represents the level of the series,
• $b_t$ represents the trend,
• $s_t$ is the seasonal component,
• $y_t$ is the observed value at time t,
• $\alpha$, $\beta$, $\gamma$ are the smoothing parameters,
• m is the number of observations per seasonal cycle,
• h is the forecast horizon,
• $k=\lfloor (h-1)/m\rfloor$.

The additive method is more appropriate when the amplitude of seasonal fluctuations remains relatively constant over time, regardless of the overall level of the series. In contrast, the multiplicative method would be preferable if the magnitude of seasonal variations increased proportionally over the series. While trend and seasonality follow similar principles to those discussed in the application of the ARIMA method, the Holt-Winters method introduces the concept of “level”. This is characterized as a weighted average that gives greater importance to the most recent data. The “level” is updated at each iteration, taking into account the trend and seasonality of the series. Although this method does not build a model with parameters (p,d,q)(P,D,Q) like SARIMA, it does generate a model with specific parameters: level smoothing ($\alpha$), trend smoothing ($\beta$) and seasonal smoothing ($\gamma$). These parameters control the rate at which the model adapts to changes in the level, trend and seasonal pattern of the series, respectively. Figure 12 illustrates that the forecast (represented by the black line) shows greater similarity to the immediately preceding period than to any other observed period, due to the influence of the “level” variable.

In the following Table 2 shows the actual values for the year 2023, the values predicted by the ARIMA method for the same year, as well as the values predicted by the Holt-Winters method for the same period. The comparison between three sets of data for the year 2023, i.e. the actual values observed, the forecasts generated by the ARIMA method and the forecasts produced by the Holt-Winters method. This Table 2 allows for a direct comparative analysis between the data actually recorded and the projections made by both forecasting methods, making it easier to assess the accuracy and effectiveness of each modeling approach.

4. Results

The section discussing the results focuses on the comparative analysis of the water consumption forecasts obtained using the Holt-Winters and ARIMA statistical methods. Figure 13 shows a graph overlaying the forecasts generated by both methods with the actual values observed. This visualization reveals that the Holt-Winters method (represented by the blue line) shows a lower error in its forecasts compared to the ARIMA method (red line). The superior effectiveness of the Holt-Winters method can be attributed to the inconsistency in the seasonality and trend patterns observed in previous years, possibly influenced by the global COVID-19 pandemic, which significantly altered consumption habits during 2020 and 2021. The Holt-Winters method manages to produce forecasts that are closer to the real values (green line) due to its characteristic of giving greater weight to the most recent data in the forecast equation. This approach allows for faster adaptation to recent changes in consumption patterns.

By graphing the actual time series (2018-2022), together with the actual consumption figures for the year 2023 and the forecasts generated by both statistical methods, in Figure 14, it becomes clear that the Holt-Winters method produced results closer to the observed reality. The line graph shows a clear visual comparison between the original time series from 2018 to 2022, actual consumption figures recorded in 2023, the forecasts generated by the ARIMA method and the forecasts produced by the Holt-Winters method. This graphical representation allows for an immediate comparative analysis, highlighting the superior accuracy of the Holt-Winters method forecasts compared to the ARIMA method when compared to the actual consumption data in 2023.

To conclude the analysis, we evaluated the accuracy of both models generated by the statistical methods SARIMA and Holt-Winters, using two complementary metrics: the Mean Absolute Error (MAE) and the Mean Squared Error (MSE). These evaluation metrics provide different perspectives on the quality of the predictions, there by enabling a more robust analysis of the models performance [32,33,34].

5. Conclusions and Future Work

This study demonstrated the successful application of ARIMA and Holt-Winters methods in forecasting monthly water consumption in Cambra, Vouzela. The detailed methodology for ARIMA, including logarithmic transformation, KPSS stationarity testing, differentiation, and (automatic) parameter selection via Auto-ARIMA, provided a robust approach to handling the complexities of the water consumption time series. The comparison between ARIMA and Holt-Winters yielded valuable conclusions regarding the relative effectiveness of these methods for this specific dataset. The use of the Mean Absolute Error (MAE) and Mean Squared Error (MSE) as an evaluation metrics allowed for an objective comparison, highlighting the strengths and limitations of each approach. Likely due to changes in water consumption habits in 2020 and 2021, influenced by the COVID-19 pandemic, the Holt-Winters method proved more effective, as its characteristic of assigning greater weight to more recent data made it better suited for forecasting. The results obtained have significant implications for water resource management in the parish of Cambra and potentially in other similar regions. However, Auto-ARIMA may prove to be the more effective method in other nearby regions, as the choice will always depend on the constructed time series, with seasonality playing a particularly relevant role. The generated forecasts can assist municipal managers in more efficiently planning water supply, providing advance information on periods of higher consumption. Based on these insights, it is possible to adjust distribution capacity, optimise resource allocation—such as personnel and infrastructure—and make informed decisions regarding maintenance or future investments. Additionally, forecast-based planning can help identify and mitigate potential water scarcity risks, preventing wastage and promoting the sustainability of water resources. For future work, it is recommended to implement telemetry using LPWAN (Low Power Wide Area Network) technology for real-time data collection, enabling more frequent and accurate updates to statistical methods, as well as developing an early warning system based on forecasts to identify potential anomalies in water consumption. This research makes a significant contribution to the field of water resource management, providing a replicable methodology and valuable insights into water consumption patterns. The continuous refinement and expansion of these models have the potential to substantially enhance the efficiency and sustainability of water usage in both urban and rural contexts.