I. Introduction

The refill drinking water depot industry (DAMIU) has become an essential service in Indonesia, supported by high community dependence on refill drinking water. Based on data from the ministry of industry in collaboration with the Indonesian Refill Water Depot Association (ASDAMINDO), there are more than 78,000 depot units operating nationally, while the Central Statistics Agency (BPS) reports that around 34,49% of households rely on refill water as their main drinking water source [1], [2]. This indicates that depot operations must ensure continuous availability and reliable fulfillment of daily demand.

As part of this ecosystem, Depot XYZ in Balikpapan operates a hybrid B2C and B2B business model with an average daily demand of around 280 gallons. However, demand patterns are highly fluctuating and do not show fully stable seasonal behavior. In operational practice, inventory management is still reactive, where stock replenishment is only performed after depletion occurs during ongoing orders, leading to delays in service fulfillment.

This reactive system has resulted in operational inefficiencies, including 40 recorded out-of-stock incidents over a three-year period. Delivery delays can reach up to several hours depending on production condition while normal lead time is only around 10-50 minutes. With the average yearly out-of-stock (OOS) frequency is 6,67%, this condition is considered significant when OOS reaching more than 4,94% in operational studies [3]. These disruptions contribute to service instability and customer dissatisfaction.

To address this issue, demand forecasting is required as data-driven approach for predicting future demand based on historical patterns. Time series forecasting enables proactive inventory planning so that stock replenishment can be scheduled before demand spikes occur. In this study, SARIMA and Prophet are selected as base models due to their capability in capturing autocorrelation, seasonality, trend, and irregular patterns in demand data [4], [5].

However, No single forecasting model is consistently optimal across all data characteristics. Recent studies show that ensemble approaches can empirically improve prediction stability by combining complementary information from multiple models [6]. Stacking ensemble has also been applied in operational forecasting contexts, such as food-demand forecasting and primary energy consumption, where prediction improvement obtained by learning complementary patterns from different base models[6], [7]. However, its use for daily drinking gallon demand forecasting in refill drinking water depot operations remains underexplored. Therefore, this study proposes a stacking ensemble using SARIMA and Prophet as base learners and XGBoost as a meta-learner, with the aim of improving forecasting accuracy and supporting inventory threshold decisions. The study focuses on daily demand data from October 2023 to Januari 2026 and evaluates performance using RMSE and MAPE [8], [9].

II. Method

A. Research Framework

This study adopts the Cross-Industry Standard Process for Data Mining, also known as CRISP-DM as the main methodological approach. CRISP-DM provides an iterative and structured workflow consisting of business understanding, data understanding, data preparation, modeling, evaluation, and deployment as pictured on the Figure 1 below [10]. Since this article focuses specifically on the forecasting process, the discussion is limited to the first five phases. The deployment phase of any form such as website or inventory implementation is outside the scope of this article.

Figure 1. The flow stages of the CRISP-DM framework [10]

B. Business Understanding

The business understanding phase identifies the operational issue faced by Depot XYZ. The depot still applies a reactive replenishment mechanism, in which stock preparation is conducted after inventory shortages occur during daily operations. This condition may increase the risk of stockout and delay order fulfillment, particularly during sudden increase in demand. Stockout prediction and inventory planning are relevant operational issues because insufficient stock availability may disrupt service performance and reduce the ability of a business to fulfill customer demand [11].

Based on this problem, the analytical objective of this study is to develop a daily gallon demand forecasting model, SARIMA dan Prophet are used as base learners, while XGBoost is applied as a meta-learner to optimize the prediction outputs generated by both models.

C. Data Understanding

The dataset consists of historical daily gallon demand records from Depot XYZ in Indonesia. The observation period covers October 1st, 2023, to January 22nd, 2026. The final dataset contains 845 daily observations, which are sufficient to represent historical demand variation across more than two years of depot operations. The dataset is structured as a univariate time series with two primary variables: the observation date and the total daily demand.

Tie series analysis requires initial understanding of trend, seasonality, and irregular variation because these components influence model selection and forecasting accuracy [9]. Therefore, exploratory analysis is conducted to identify temporal patterns, daily fluctuations, missing dates, and potential weekly seasonality before the modeling process. Thus, the main variables used in this research is as shown on the Table 1 below.

Table 1. Research Variables

D. Data Preparation

The data preparation phase transforms raw operational records into a consistent daily time series format. First, the transaction date is converted into a standardized datetime format. Second, the observations are sorted chronologically top reserve their temporal order. Third, a continuous daily index is constructed to identify dates without recorded demand values. Missing values generated during this process are handled using temporal interpolation to maintain data continuity.

The processed dataset is divided chronologically into training and testing sets using an 80:20 ratio. The first 80% of the observations are used for model training, while the remaining 20% are used for final out-of-sample evaluation. A chronological split is applied instead of random sampling because time-series forecasting must preserve the temporal order between past and future observations. This approach Is consistent with forecasting validation principles, where the testing data should represent future observations that are not used during model training. Previous studies also emphasize that validation design has a strong influence on forecasting performance, particularly when comparing fixed split and walk-forward validation techniques in ARIMA-based forecasting[12], [13]

E. Modelling

1. SARIMA Model

SARIMA (Seasonal Autoregressive Integrated Moving Average) is used as the first base learner. SARIMA extends the ARIMA model by incorporating seasonal components, making it suitable for time series data with recurring patterns. The model is generally represented as::

Where p,d,q represents the non-seasonal components of autoregressive, differencing, and moving-average orders, while P,D,Q represent the seasonal orders. The parameter s represents the seasonal period.

Previous studies indicate that SARIMA remains effective when temporal data can be transformed into a stationary structure and contain identifiable seasonal patterns [4], [5]. In this study, SARIMA parameter selection supported by stationarity testing, ACF and PACF analysis, and information criteria such as AIC and BIC. Manual parameter identification is also compared with the Auto-SARIMA approach to select the most suitable configuration.

2. Prophet Model

Prophet is applied as the second base learner. The model uses a decomposable additive structure consisting of trend, seasonality, holiday effects, and an error term [14]. Prophet is generally expressed as:

Where g(t) represents the trend component, s(t) represents seasonality, h(t) represents holiday effects, and ϵ_t represents the error term.

Prophet is included because its additive structure provides a complementary approach to SARIMA. While SARIMA focuses on autocorrelation and stationarity, Prophet models trend and seasonal behavior more explicitly. This distinction is relevant because empirical studies show that relative performance of SARIMA and Prophet depends on the characteristics of the observed time series [4], [5]..

3. XGBoost as Meta-learner

The prediction outputs from SARIMA and prophet are combined using a stacking ensemble approach. In this framework, both forecasting models act as base learners, while XGBoost acts as the meta-learner. The general formulation is expressed as:

Where represents the final ensemble prediction, while and represent the prediction outputs generated by the two base learners.

Stacking is used because ensemble performance can improve when the meta-learner receives complementary information from both models with different forecasting structures [6], [15]. Jiang et al. also demonstrate that a meta-learning approach can improve prediction stability in short-term forecasting tasks with limited and fluctuating data [16]. In this study, XGBoost is selected because its tree-based boosting structure can model nonlinear relationships between base-model outputs while controlling complexity through regularization [17]

To reduce the risk of data leakage, the meta-learner is trained using out-of-fold predictions generated through time-series validation. In each fold, the base models are trained on earlier observations and used to predict the subsequent validation period. The resulting predictions are collected as meta-features, while the actual demand values are used as the targeted variable. After the meta-learner is trained, SARIMA and Prophet are refitted using the full training set and used to predict the testing period.

This validation strategy is used to ensure that the meta-learner is trained on predictions generated from temporally separated validation periods. Time-series validation is relevant for forecasting tasks because it prevents future observations from being used during model training and provides a more realistic evaluation structure than random cross-validation[13].

The meta-learner was configured using a relatively conservative parameter setting to reduce the risk of overfitting. The model employed 150 estimators with learning rate of 0.05. A maximum tree depth of 1 was selected to limit model complexity, while a subsample ratio of 0.6 was used so that each boosting iteration trained only on a portion of the available observations. The gamma parameter was set to 5 to restrict unnecessary tree splitting. In addition, regularization was strengthened by setting the L2 regularization parameter to 10 and the L1 to 5. These settings were applied to ensure that the meta-learner focused on stable correction patterns from the SARIMA and Prophet outputs rather than fitting short-term noise in the training data.

F. Evaluation

The evaluation phase compares the forecasting performance of SARIMA, Prophet, and the XGBoost-based stacking ensemble. Root Mean Square Error (RMSE) and Mean Absolute Percentage Error (MAPE) are used as the primary evaluation metrics. RMSE measures the magnitude of prediction errors in the original unit, while MAPE expresses prediction errors as percentages [9].

RMSE is calculated as:

MAPE is calculated as:

III. Results a nd Discussion

A. Data Preparation and Stationarity Analysis

The preprocessing stage produced a continuous daily time series containing historical gallon demand data from October 1st, 2023, to January 22nd, 2026. The raw operational records were sorted chronologically, completed into a daily index, and processed using temporal interpolation to handle missing values. The historical visualization showed fluctuating daily demand with recurring weekly pattern as shown in Figure 2. Therefore, a seasonal period or 7 days to represent weekly seasonality in the SARIMA model.

Figure 2.Historical Daily Gallon Demand Time Series.

Before model estimation, stationarity was evaluated using the Augmented Dickey-Fuller (ADF) and Kwiatkowski-Philips-Schimdt-Shin (KPSS) tests. After first-order differencing, the ADF statistic reached -9.7485 with a p-value of 0.0000, while the KPSS reached 0.3337 with a p-value of 0.1000. These results indicate that the differenced series satisfies the stationarity requirement because the ADF test rejects the unit-root hypothesis and the KPSS test does not reject the stationarity hypothesis. Stationarity is important for SARIMA because the model estimates temporal relationships more effectively when the statistical structure of the series is sufficiently stable [9]. Detailed statistical output of this process is as shown in Table 2 below..

Table 2. Stationarity Test After First-Order Differencing

The ACF and PACF plots were subsequently used to identify initial SARIMA parameter candidates. The PACF plot showed significant spikes at early lags, particularly lags 1, 2, 3, and 5, while the ACF plot showed a cut-off pattern after lag 3. These patterns indicate that the demand series contains short-term temporal dependence. However, parameter selection was not based exclusively on visual interpretation because several lags remained significant. Therefore, manual candidates were compared using AIC and BIC, and an automated parameter search was also conducted using Auto-SARIMA.

B. SARIMA Model Results

The manual SARIMA identification process produced several candidate configurations based on the ACF, PACF, and weekly seasonal period. Among the manually tested candidates, the selected configuration was SARIMA(3,1,1)(0,1,1)(7). In parallel, the Auto-SARIMA procedure identified SARIMA(2,1,3)(1,0,0)(7) as the most suitable automated configuration, with an AIC value of 6429.889 and a BIC value of 6461.4982.

The out-of-sample evaluation showed that the automated configuration performed better than the manually selected model. The manual SARIMA model produced an RMSE of 75.9720 and a MAPE of 19.42%, whereas Auto-SARIMA achieved an RMSE of 74.200307 and a MAPE of 18.49%. Since the final model selection was based on forecasting error on the testing period, Auto-SARIMA parameter was selected as the SARIMA base learner. The comparison between manual SARIMA and Auto-SARIMA is presented in Table 3 below.

Table 3. SARIMA Model Evaluation Results

The stronger performance auto-SARIMA suggests that the demand series contains a temporal structure that can be effectively represented through autoregressive, moving-average, and weekly seasonal components. This finding is consistent with Gunawan and Ramadani [5], who reported that SARIMA outperformed Prophet in inventory forecasting when the seasonal structure of the data could be modeled more effectively. It also aligned with Kenyi and Yamamoto [18], who emphasized that SARIMA remains suitable when the series contains identifiable temporal dependence and can be transformed into a stationary structure.

C. Prophet Model Results

Prophet was developed as the second base learner to provide a forecasting structure that differs from SARIMA. The model was configured using weekly and yearly seasonality components to capture recurring temporal patterns. Unlike SARIMA, which relies on autocorrelation and differencing, Prophet decomposes the series into trend, seasonality, holiday effects, and residual components [14].

Table 4. Base Learner Evaluation Results

The performance difference may indicate that the dominant information in the dataset is more closely related to short-term autocorrelation and weekly seasonality than to smoother trend changes. Prophet is designed to flexibly model trend and seasonal components, but this flexibility does not guarantee superior performance when daily fluctuations and short-term variations dominate the observed series. This interpretation is consistent with previous comparative studies showing that Prophet does not automatically outperform statistical models and that forecasting accuracy remains dependent on the characteristics of the dataset [5], [18].

D. XGBoost-Based Stacking Ensemble Results

The SARIMA and Prophet predictions were subsequently used as meta-features for the XGBoost-based stacking ensemble. To avoid data leakage, the meta-training dataset was constructed using out-of-fold predictions generated through five-fold time-series validation. In each fold, the base learners were trained using earlier observations and evaluated on the subsequent validation period. The resulting SARIMA and Prophet predictions were then used as inputs for the XGBoost meta-learner.

The final stacking model achieved and RMSE of 73.648907 and a MAPE of 18.32%. These values are lower than the error values produced by both base learners. Therefore, the XGBoost-based stacking ensemble provides the best numerical performance among the evaluated models.

Table 5. Final Model Evaluation Results

Compared with SARIMA as the strongest base learner, the stacking model reduced RMSE by 0.551400 and MAPE by 0.17%. These results show that the ensemble model improved forecasting accuracy, although the magnitude of improvement was relatively limited.

E. Discussion

The limited improvement obtained from the stacking model can be explained through the principle of base-learner diversity. An ensemble does not automatically produce a substantial performance improvement simply because several models are combined. Its effectiveness depends on whether the base learners provide complementary information and generate sufficiently different error patterns. When the base models capture different aspects of the data, the meta-learner has more useful information for correcting prediction errors [7], [19].

In this study, SARIMA and Prophet were selected because they represent different forecasting structures. SARIMA captures autoregressive and seasonal relationships, while Prophet models trend and seasonality through additive decomposition. However, the evaluation results show a clear performance imbalance: SARIMA produced a MAPE of 18.49%, whereas Prophet produced a MAPE of 25.89%. This difference may indicate that SARIMA already captured most of the dominant temporal information in the dataset, while Prophet provided only limited additional information for the meta-learner.

As a result, XGBoost was still able to learn correction patterns from the two base learners, but the available room for improvement was relatively narrow. The meta-learner reduced prediction error because it did not simply average the two outputs; instead, it learned how to adjust their contributions based on patterns observed during out-of-fold training. However, the weaker contribution from Prophet limited the potential improvement over the strongest base model.

This finding is important because it presents a more realistic interpretation of stacking ensemble performance. Jiang et al. [16] demonstrated that meta-learning can improve forecasting stability under fluctuating data conditions, while Wang and Dong [19] and Divina et al. [7] emphasized the importance of combining diverse learners. The present study supports these findings but also shows that ensemble improvement may remain marginal when the base learners do not contribute equally strong and complementary information.

From a practical perspective, the stacking ensemble is preferable when the primary objective is to minimize forecasting error. However, SARIMA remains a competitive alternative if implementation simplicity and lower computational complexity are prioritized. The main contribution of this study lies in demonstrating that the optimization process must be evaluated empirically rather than assumed to provide a large improvement in every forecasting context.

IV. Conclusions

This study applied SARIMA and Prophet as base learners and XGBoost as a meta-learner to forecast daily gallon demand at Depot XYZ. The results show that SARIMA achieved an RMSE of 74.200307 and a MAPE of 18.49%, while Prophet produced an RMSE of 92.954735 and a MAPE of 25.89%. The XGBoost-based stacking ensemble generated the best forecasting performance, with an RMSE of 73.648907 and a MAPE of 18.32%. Although the stacking ensemble reduced the prediction error compared with the strongest base learner, the improvement remained limited, with a reduction of 0.551400 in RMSE and 0.17 percentage points in MAPE. This finding indicates that the effectiveness of stacking depends not only on the use of multiple models but also on the diversity and complementary information provided by the base learners. In this case, SARIMA already captured most of the dominant temporal patterns in the dataset, while Prophet provided a more limited additional contribution to the meta-learner. Therefore, the XGBoost-based stacking ensemble can be selected as the final forecasting model when the primary objective is to minimize prediction error. However, SARIMA remains a competitive alternative when implementation simplicity and computational efficiency are prioritized. Future studies may incorporate external variables, such as holidays, weather conditions, or customer-related characteristics, to provide more diverse information and further improve forecasting accuracy.