## 2 Methodology
### SARIMA Module  
In the context of our proposed method, the SARIMA module can be defined by a mathematical formula that represents its components: autoregressive (AR), differencing (I), moving average (MA), and the seasonal components of each. The generic form of a SARIMA model is typically denoted as follows:  
\[SARIMA(p,d,q)\times(P,D,Q)_{s} \tag{30}\]  
where p is the number of autoregressive terms, d is the number of non-seasonal differences needed for stationarity, q is the number of lagged forecast errors in the prediction equation, P is the number of seasonal autoregressive terms, D is the number of seasonal differences, Q is the number of seasonal moving average terms, and s is the number of periods in a season. Even in seasonally inactive scenarios, Equation (30) provides a general framework that ensures the model is robust across various datasets, including those where seasonality is not a major factor. It serves as a baseline for comparing more complex, seasonally active systems. Specifically, within our model, the formula for the SARIMA module can be expressed as follows:  
\[\phi(B)(1-B)^{d}(1-B^{s})^{D}X_{t}=\theta(B)\theta(B^{s})Z_{t} \tag{31}\]  
where \(\phi(B)\) represents the non-seasonal autoregressive operator, \(\theta(B)\) presents the non-seasonal moving average operator, \(\Theta(B^{s})\) is the seasonal moving average operator, \(X_{t}\) is the time series being modeled, and \(Z_{t}\) is the white noise series. The above equation can be further decomposed into the following:  
\[\phi(B)(1-B)^{d}(1-B^{s})^{D}X_{t}=\theta(B)\epsilon_{t} \tag{32}\]  
\[\theta(B^{s})\epsilon_{t}=Z_{t} \tag{33}\]  
For the SARIMA module proposed in this paper, its establishment mainly includes the following steps:  
(1) Model identification: According to the training dataset, the appropriate parameters for the corresponding IMFs are automatically determined. This process is automatically completed by fitting various combinations of ARIMA models for different orders \((p,d,q)\) and seasonal orders \((P,D,Q)_{s}\). The models are evaluated based on a given information criterion (the Akaike Information Criterion), and the model with the lowest criterion value is selected. (2) Using the parameters determined in the previous step, future values of the time series X are predicted with the SARIMA model. The prediction is conducted in a rolling forecast manner, that is, predicting a certain number of future steps, with the model being updated with new data at each step. (3) The model is continuously predicted and updated upon the receipt of new data. This iterative process continues until predictions have been made for the entire forecast range. We also considered several boundary conditions during the application of the SARIMA model to ensure the accuracy and relevance of the predictions. These boundary conditions include the following:  
1. Stationarity: The SARIMA model assumes that the time series data are stationary. To satisfy this condition, we employed the auto arima function, which automatically determines the differencing required to achieve stationarity if the raw data exhibit non-stationary behavior. 2. Autocorrelation: The autoregressive component of the SARIMA model requires that the time series data exhibit significant autocorrelation. Through the model selection process, auto arima tests various parameter combinations to ensure that the chosen model has sufficient autocorrelation for effective forecasting. 3. Short-term seasonality: We explicitly modeled short-term seasonality by specifying a weekly cycle (m = 7). This reflects the recurring patterns in subway passenger flow that occur within a week, such as weekday versus weekend differences. The seasonal component of the SARIMA model captures these variations to improve forecast accuracy. In summary, in this paper, the original subway passenger flow data are first decomposed into multiple Intrinsic Mode Functions (IMFs) using the VMD method.