Demand Prediction for Food and Beverage SMEs Using SARIMAX and Weather Data

SARIMAX stands for Seasonal AutoRegressive Integrated Moving Average with eXogenous variables. It's a type of time series model that combines autoregressive, differencing (I), and moving average operations, as well as accounting for seasonality (seasonal ARIMA) and exogenous variables. It's used for forecasting when data has a trend and/or a seasonal pattern, and when external, independent variables might influence the time series variable being forecasted. It expands on simpler models like AR and MA by including trends and seasonality, allowing for more accurate predictions in these complex contexts.

The components of SARIMAX are:

(1) Seasonal (S): SARIMAX accounts for seasonality in the data. Seasonal patterns occur when there are regular fluctuations or patterns in the data that repeat over specific periods.

(2) Autoregressive (AR): This component captures the relationship between a variable and its own lagged values. In other words, it models the relationship between an observation and several lagged observations (from previous time steps).

(3) Integrated (I): The "integrated" part of SARIMAX indicates the differencing of raw observations to make the time series stationary. Stationary time series have constant statistical properties over time, making them easier to model.

(4) Moving Average (MA): This component models the relationship between a variable and a residual error from a moving average model applied to lagged observations.

(5) Exogenous Regressors (X): SARIMAX allows for the inclusion of exogenous variables, which are external factors that can influence the time series being analyzed. These variables are not directly related to the time series but are believed to have an impact on it. Including exogenous regressors in the model can improve its forecasting accuracy, especially when there are known external factors affecting the time series.

And the equation for SARIMAX (p,d,q)(P,D,Q,s) can be seen as follows:

$\Theta(\mathrm{L})^{\boldsymbol{p}} \theta\left(\mathrm{L}^{\boldsymbol{s}}\right)^{\boldsymbol{P}} \Delta^{\boldsymbol{d}} \Delta_{\boldsymbol{s}}^{\boldsymbol{D}} y_t=\Phi(\mathrm{L})^{\boldsymbol{q}} \boldsymbol{\phi}\left(\mathrm{L}^{\boldsymbol{s}}\right)^{\boldsymbol{Q}} \Delta^{\boldsymbol{d}} \Delta_{\boldsymbol{s}}^{\boldsymbol{D}} \epsilon_t+\sum_{i=1}^n \beta_i x_t^i$            (1)

where, $y_t$ is given time series, $p$ is the number of time lags to regress on, $\epsilon_t$ is the noise at time $t, \beta$ is constant, $L$ is regular lag operator, $\Theta(\mathrm{L})^p$ are order $p$ polynomial function of $L$, $q$ is the number of time lags of error to regress on, $\phi$ is defined analogously to $\Theta$, $\Delta^d$ is integration operator where $d$ is the order of differencing used, $\Delta_s^D$ is seasonal differences of the time series where $s$ is the number of time lags comprising one full period of seasonality and $D$ order differencing which applied to seasonal lags. Then, $n$ exogeneous variables defined at each time step $t$, denoted by $x_t^i$ for $i \leq n$.

The reason SARIMAX is chosen due to: 1) SARIMAX models allow for the inclusion of exogenous variables, such as weather features, which can significantly impact demand patterns. By incorporating these variables, SARIMAX models can capture the influence of external factors on demand and improve the accuracy of predictions [27]; 2) SARIMAX models are well-suited for capturing seasonal patterns and autocorrelation in demand data. The seasonal and autoregressive components of SARIMAX models enable them to model and forecast demand accurately, especially when there are recurring patterns and dependencies in the data [28]; 3) SARIMAX models offer flexibility in terms of model specification and customization. They can be tailored to specific demand patterns and adjusted to incorporate different exogenous variables, allowing for a more precise and tailored demand prediction [13]; 4) SARIMAX models provide interpretable results, allowing for a better understanding of the underlying factors driving demand. The coefficients and parameters of SARIMAX models can provide insights into the relationships between demand and exogenous variables, aiding in decision-making and strategy development [13].

One potential weakness is the sensitivity of the model to the selection of its parameters, which may require careful tuning and optimization for different types of demand data [29]. Additionally, SARIMAX may be less effective when dealing with highly volatile or irregular demand patterns, as it relies on capturing and modeling the underlying time series components, which may be challenging in such cases [18].

Vector Autoregression (VAR) is a type of statistical modeling used in econometrics that captures the linear interdependencies among multiple time series. It allows each variable in the system to be a function of the lagged values of all other variables in the system, which makes it crucial in forecasting systems of interrelated time series variables. Essentially, each variable in a VAR is modeled as a linear combination of past values of itself and the past values of all other variables in the system. The primary advantage of VAR is its ability to model the dynamic relationships amongst multiple (typically economic) variables simultaneously.

The equation for vector autoregression can be seen as follows:

$y_t=c+A_1 y_{t-1}+A_2 y_{t-2}+\cdots+A_p y_{t-p}+\varepsilon_t$             (2)

where, $y_t$ is $k \times 1$ vector representing the variable at time $t$. $c$ is $k \times 1$ vector of constants. $A_1$, $A_2$, $\ldots$, $A_p$ are $k \times k$ matrices of coefficients representing the lagged values of the variables. $p$ is the order of the VAR model, indicating how many lagged values are included in the model. Then, $\varepsilon_t$ is $k \times 1$ vector representing the error term at time $t$.

Also, the reason VAR is chosen due to: 1) VAR models are capable of capturing the interdependencies and dynamic relationships among multiple variables simultaneously. This makes them suitable for demand prediction, where multiple factors can influence the demand patterns [19]; 2) VAR models allow for the analysis and forecasting of multiple variables together, considering their mutual interactions. This is particularly useful when there are feedback loops and interrelationships among the variables, which can affect demand dynamics [30]; VAR models can incorporate exogenous variables, such as economic indicators, which can have an impact on demand. By including these variables, VAR models can capture the influence of external factors on demand patterns and improve the accuracy of predictions [31]. One potential weakness of VAR is the challenge of estimating a large number of parameters, especially when dealing with high-dimensional data or a large number of variables. This can lead to increased computational complexity and potential issues related to overfitting, particularly when the number of observations is limited [32].

To obtain optimal parameter for SARIMAX we utilized Akaike Information Criterion (AIC). Then, for VAR, there are several considerations, such as Akaike Information Criterion (AIC), Bayesian Information Criterion (BIC), Final Prediction Error (FPE) and Hannan-Quinn Information Criterion (HQIC).

The Akaike Information Criterion (AIC) is a statistical tool utilized to assess how well a statistical model fits the data. It strikes a balance between the model's fit quality and its complexity, discouraging overly intricate models that might overfit the data. When applied to Vector Autoregression (VAR) models, AIC helps compare various models, aiding in the selection of the most suitable one for forecasting or analyzing complex time series data with multiple variables. A lower AIC value indicates a better trade-off between goodness of fit and model complexity. When comparing different VAR models, you would choose the model with the lowest AIC value, as it suggests that the model provides a good fit to the data without being overly complex. The AIC equation can be seen as follows:

$A I C=-2 \log (L)+2 p$               (3)

where, $L$ is the likelihood of the data given the model. It represents how well the model explains the observed data. Then, $p$ is the number of parameters in the VAR(p) model, which includes the autoregressive coefficients and the error variances.

The Bayesian Information Criterion (BIC) helps balance the goodness of fit of a model with its complexity, penalizing models with more parameters. The BIC penalizes the number of parameters more heavily than the AIC, which means that it tends to favor simpler models. A lower BIC value indicates a better trade-off between goodness of fit and model complexity. The BIC equation can be seen as follows:

$B I C=-2 \log (L)+p \times \log (n)$                    (4)

In here, $n$ is the number of observations in the time series data.

The Final Prediction Error (FPE) is used to estimate the expected mean squared prediction error of a model. FPE provides a measure of forecast accuracy, penalizing complex models to avoid overfitting. Like AIC and BIC, a lower FPE value indicates a better trade-off between goodness of fit and model complexity. The FPE equation can be seen as follows:

$F P E=\left(\frac{n+p+1}{n-p-1}\right) \times\left|\widehat{\sum u}\right|^{-1}$          (5)

In here, $\widehat{\sum u}$ is the estimated residual covariance matrix of the VAR(p) model.

The Hannan-Quinn Information Criterion (HQIC) provides a trade-off between the goodness of fit of the model and its complexity. The HQIC incorporates a logarithmic term related to the sample size $(n)$ that penalizes the number of parameters in the model. Like AIC and BIC, a lower HQIC value indicates a better trade-off between goodness of fit and model complexity. The equation can be seen as follows:

$H Q I C=-2 \log (L)+2 \times$ params $\times \log (\log (n))$            (6)

In here, params is the number of parameters in the VAR model, including the autoregressive coefficients and the error variances.

After gaining the optimal parameter. The performance of algorithm can be measure with several measurements to understand which methods are better than other, such as: MSE, MAE and RMSE. MSE, Mean Squared Error (MSE) calculates the average of the squared variances between our predictions and the real values. This metric gauge the effectiveness of our prediction model. Smaller MSE values suggest that predictions closely align with the real values, whereas larger MSE values indicate greater discrepancies between predictions and reality which we can see in following equation:

$M A E=\frac{\sum_{i=1}^n\left|\widehat{Y}_l-Y_i\right|}{n}$        (7)

Mean Absolute Error (MAE) is the average of the absolute differences between your predictions and the actual values. MAE provides a straightforward and easy-to-understand way of measuring the accuracy of your predictions. Like MSE, lower MAE values indicate that the predictions are closer to the actual values and the MAE equation can be seen as follows:

$M S E=\frac{1}{n} \sum_{i=1}^n\left(Y_i-\widehat{Y}_l\right)^2$                    （8）

RMSE is a way of calculating the average size of the errors in our predictions, with the square root ensuring the scale of the error metric matches the scale of the original values. Like MSE and MAE, lower RMSE values indicate that the predictions are closer to the actual values, signifying a more accurate prediction model. The equation can be seen as follows:

$R M S E=\sqrt{\frac{\sum_i^N\|Y(i)-\hat{Y}(i)\|^2}{n}}$             （9）

For MAE, MSE and RMSE, the symbol can be read as: Where $n$ denotes number of data points, $Y_i$ denotes observed values and $\widehat{Y}_l$ denotes predicted values.

The optimal results for SARIMX can be seen in Tables 4 and 5 as follows:

Table 4. SARIMAX model parameter

The Table 2 contains information about the parameters used and put into model. Thus, the model Number is used as reference to Table 3 as follows:

Table 5. SARIMAX measurement results

From Table 5 above, we chose the lowest AIC (1069.31) model number 2 where the d is 0, D is 1, p is 0, q is 0, P is 2, and Q is 2. Another lowest AIC is model number 4 (1070.11 and obtain better AIC, MSE, MAE and RMSE than before) where d is 1, D is 1, p is 2, q is 3, P is 2 and Q is 2. Although both is not gaining better measurement (MSE, MAE or RMSE) with other model, our primary goal is to understand the underlying patterns in the data, which AIC might be a more appropriate metric.

Thus, rather chose the lowest AIC, we took equilibrium (two lowest value and compare the MSE, MAE, RMSE), we can conclude that our dataset with (d=1, D=1, p=2, q=3, P=2, Q=2) needs differencing (d) and seasonal differencing (D) to become stationary (d), it needs 2 past observations (p), size of the moving average window is set to 2 which signifies the number of lagged forecast errors in the prediction equation (q), seasonal lag order is set to 2 (P), and size of the seasonal moving average window is set to 2.

For VAR, the VAR order can be seen as follows:

Table 6. VAR order selection (*highlights the minimum)

From Table 6 above, the lag number value for AIC, BIC FPE and HQIC is different. Thus, here we want to know, which order is better for our case (the sales qty is taken for best measurement).

Table 7. VAR measurement results

Thus, from Table 7 above, we can conclude that the lag order chosen is 8 due to lowest AIC.

Finally, we compare SARIMAX with VAR to obtain the conclusion which one is better and presented in Table 8.

Table 8. Comparison measurement results

From Table 8, we can conclude that SARIMAX is suitable for demand prediction regarding temp avg and humidity avg due to have lowest value in all measurements.