Hierarchical Spatio-Temporal Transformer for Coherent Emergency Department Forecasting Filipa Lino1 ( ), Bárbara Tavares1 , Carlos Santiago1 , Cláudia Soares2 , and Manuel Marques1
arXiv:2607.27106v1 [cs.LG] 29 Jul 2026
1
Institute for Systems and Robotics, LARSyS, Instituto Superior Técnico, Portugal 2 NOVA School of Science and Technology, Caparica, Portugal [email protected]
Abstract. Emergency Departments (EDs) are critical access points in healthcare systems, yet they face persistent pressure from unpredictable patient demand, seasonal surges, and non-urgent visits. Effective ED planning requires forecasts at multiple decision-making levels: hospitals need local demand estimates for staffing and bed management, regions require forecasts to coordinate healthcare units, and national authorities need system-wide projections for capacity planning. However, most existing approaches forecast ED demand independently at a single level, ignoring the hierarchy linking hospitals, regions, and national systems. This can produce incoherent predictions, where hospital-level forecasts do not aggregate consistently to regional or national demand. We propose HierSTT, a hierarchical Transformer-based framework for coherent multi-level ED forecasting. HierSTT jointly predicts hospital, regional, and national level demand in a single end-to-end model. A Temporal Fusion Transformer captures national dynamics, while spatio-temporal Transformer encoder-decoder modules model regional and hospital demand conditioned on higher-level forecasts. A coherence-aware loss penalizes cross-level inconsistencies during training. We further introduce a nationwide Portuguese ED dataset covering 81 hospitals across 5 regional health administrations, with heterogeneous covariates at each level. Experiments show that HierSTT reduces average WAPE by 32% relative to the best non-hierarchical deep learning baseline and outperforms all classical hierarchical reconciliation methods, while producing near-coherent predictions across levels. Additional resources associated with this work are available at https://github.com/FilipaLino/HierSTT. Keywords: Emergency Department Forecasting · Hierarchical Time Series · Transformers · Spatio-Temporal Modeling · Forecast Coherence.
1
Introduction
Emergency Departments (EDs) play a critical role in the healthcare system, serving as the largest source of hospital admissions. These departments operate under persistent pressure due to unpredictable fluctuations in patient demand [27],
2
F. Lino et al.
(a)
(b)
Fig. 1: Geographical and hierarchical organization of the ED forecasting dataset. (a) Spatial distribution of EDs in Portugal, grouped by RHAs. Each point represents a hospital. (b) Hierarchical structure of the dataset.
which when exceeding operational capacity leads to overcrowding, a significant challenge [20] causing negative patient outcomes, longer waiting times, increased costs, and staff burnout. Accurate forecasting of ED demand is therefore crucial for improving patient care, optimizing resource allocation, and supporting policy decision-making [25]. Forecasting ED visits is challenging due to seasonality, weekday effects, holidays, epidemics, weather patterns, and regional mobility dynamics [26]. Traditional approaches such as autoregressive integrated moving average (ARIMA), exponential smoothing (ETS), and other statistical methods have been widely applied [26,32], while more recently deep learning models, including long shortterm memory (LSTM), and transformer-based architectures, have demonstrated improved capacity to capture nonlinear temporal dependencies [17,19,18]. Despite these advances, most existing studies formulate ED forecasting as an isolated single-level problem, focusing independently on one hospital or one region. However, healthcare systems naturally exhibit a hierarchical structure: individual hospitals belong to regional administrative units, which in turn compose the national healthcare network (Fig. 1). Ignoring this structure may lead to incoherent forecasts, where prediction generated independently at different levels are mutually inconsistent. Furthermore, data across levels is often heterogeneous. While hospital data may capture operational and demographic factors, regional and national levels provide broader environmental and population-level signals, making joint modeling desirable. Hierarchical time-series forecasting addresses this challenge by producing consistent forecasts across aggregation levels [14,29]. Classical reconciliation meth-
Hierarchical Spatio-Temporal Transformer for Coherent ED Forecasting
3
ods adjust independently generated predictions post-hoc [22], but fail to exploit cross-level dependencies during learning. To address these limitations, we propose HierSTT, a hierarchical Transformerbased deep learning framework for coherent multi-level ED demand forecasting. The model simultaneously predicts patient visits over a 4-week horizon at hospital (81), Regional Health Areas (RHAs) (5 regions), and national levels of the Portuguese National Health Service. It combines a Temporal Fusion Transformer [6] at the national level with spatio-temporal Transformer encoders and decoders [10] at regional and hospital levels. Predictions are generated in a topdown manner, enabling higher-level forecasts to inform lower-level estimates, simultaneously integrating historical demand with heterogeneous level-specific covariates. A coherence-aware loss penalizes inconsistencies between aggregated lower-level predictions and upper-level forecasts during training. Alongside the model, we introduce a nationwide Portuguese ED dataset reflecting real-world heterogeneous data availability across administrative levels. Extensive experiments demonstrate that HierSTT achieves superior balance between predictive accuracy and hierarchical coherence compared to statistical baselines, non-hierarchical deep learning models, and reconciliation methods. In summary, our main contributions are the following: – A hierarchical Transformer-based framework for simultaneous ED demand forecasting across multiple levels; – A coherence-aware training objective enforcing cross-level consistency; – A new nationwide Portuguese ED dataset with heterogeneous, level-specific covariates.
2
Related Work
2.1
Classical Approaches for ED Forecasting
Classical time-series forecasting methods have long been used for demand prediction. Naïve forecasting assumes that future values follow the most recent observations, providing a strong baseline in persistent series [8]. Yule et al. [31] introduced auto-regressive (AR) and moving-average (MA) processes, later unified under the ARIMA framework [4], which models temporal dependencies through linear relationships in past observations and residuals. Exponential Smoothing methods, such as ETS by Hyndman et al. [9], decompose time series into Error, Trend and Seasonality components, automatically selecting most suitable additive formulation for different temporal patterns [7]. In ED patient volume forecasting, these methods have been adopted due to its strong temporal structure, which follows calendar-driven patterns such as weekdays, holidays, and seasonal trends [3]. Classical statistical approaches, such as moving averages and ARIMA, have therefore been used to model these dynamics [16]. Their interpretability and computational efficiency make them attractive for operational settings.
4
F. Lino et al.
However, ED demand is influenced by complex and interactive factors, including patient demographics, epidemiological trends, environmental conditions, and operational characteristics. The nonlinear and dynamic nature of these relationships challenges the assumptions of linearity and stationarity of classical models, often limiting their predictive performance in real-world scenarios [12]. 2.2
Deep Learning for Time Series
Deep learning (DL) has become a major research direction in time series forecasting due to its ability to learn complex nonlinear patterns and multivariate interactions [13]. In ED forecasting, Kadri et al. [11] applied Recurrent Neural Networks (RNNs) to patient-flow prediction, but RNNs often suffer from training instability and limited ability to capture long-range temporal dependencies. To address these limitations, more advanced recurrent architectures have been proposed. Kashani et al. [19] showed that Long Short-Term Memory (LSTM) networks improve performance in highly correlated temporal settings. Building on recurrent models, sequence-to-sequence (Seq2seq) architectures became relevant for forecasting by encoding historical information into latent representations and decoding future horizons. Seq2Seq models based on Gated Recurrent Units (GRUs) and LSTMs have shown strong performance in capturing long-range temporal dependencies and multistep forecasting patterns [30,15]. A key advantage of DL models is their ability to integrate multiple exogenous variables. Several studies have shown that calendar-related features, weather conditions, and air quality influence ED demand variability [6,5,2]. These findings have motivated the use of attention mechanisms to dynamically weight relevant inputs. Chaves et al. [2] demonstrated the effectiveness of attention-based architectures for forecasting ED patient admissions. Transformer-based architectures have advanced time series forecasting through self-attention. In the ED domain, Ming et al. [18] demonstrated the strong predictive performance of Transformers, while Ji et al. [10] proposed a spatio-temporal variant for jointly modeling temporal and spatial dependencies. Temporal Fusion Transformers (TFT) have shown promising results in healthcare, with Pulkkinen et al. [24] applying TFT to short-term hospital demand prediction and Caldas et al. [6] extending it to longer-term forecasting across Portugal’s RHAs. Alongside Transformer-based approaches, alternative deep forecasting architectures have also emerged. Oreshkin et al. [23] introduced N-BEATS, a deep residual architecture that achieved strong forecasting performance without recurrent or attention mechanisms. However, despite these advances, most DL-based approaches remain designed for single-level forecasting and do not explicitly model hierarchical dependencies across healthcare system levels. 2.3
Hierarchical Time Series Forecasting
Hierarchical Time Series (HTS) forecasting generates predictions across multiple aggregation levels while ensuring coherence, i.e., that lower-level forecasts aggregate consistently to higher-level ones. HTS has been applied in domains
Hierarchical Spatio-Temporal Transformer for Coherent ED Forecasting
5
such as tourism demand forecasting across national and regional levels [22] and retail sales forecasting across items, departments, and product categories [21]. Classical HTS approaches include bottom-up, top-down, middle-out, and reconciliation-based strategies, which usually generate forecasts independently before enforcing coherence through post-hoc aggregation constraints [22]. However, post-hoc reconciliation limits the ability to learn cross-level dependencies and may propagate errors across the hierarchy. Recent work have sought to integrate reconciliation into model training. Mancuso et al. [14] proposed a deep neural network with a reconciliation layer and a coherence-aware loss, reconciling lower-level predictions using outputs from immediately higher level, through sequential training across hierarchical level pairs. Wang et al. [29] introduced the Neural Reconciler, combining encoderdecoder architectures with Normalizing Flows to model hierarchical forecasts. However, this method still follows a multi-stage pipeline, where base forecasts are first generated and later refined through reconciliation. Thus, hierarchical constraints are not learned jointly with temporal representations, potentially leading to suboptimal performance. Overall, while ED forecasting has advanced through statistical and deep learning approaches, the integration of hierarchical structure and heterogeneous covariates remains underexplored. To the best of our knowledge, no prior work jointly addresses hierarchical coherence and heterogeneous multi-level covariates in ED forecasting.
3
Dataset
3.1
Data Sources and Collection
This work is based on a nationwide dataset of ED activity in Portugal, comprising historical records collected across multiple administrative levels of the healthcare system. Data were initially available for 91 hospitals. Ten recently established hospitals were excluded due to insufficient data, resulting in a final dataset of 81 hospitals, grouped into 5 RHAs, as well as aggregated nationallevel indicators. This hierarchical structure is detailed in Fig. 1b. The data span the period from January 1, 2021 to April 20, 2024. The Portuguese Ministry of Health provides open-source data on the National Health Service (SNS) activities in the Transparência Website 3 and daily ED monitoring platform 4 . These sources were combined and cross-referenced with external data, including Portugal air quality data 5 , Temperature 6 and Mortality 7 to construct a unified dataset. 3
https://transparencia.sns.gov.pt/explore/?sort=modified https://www.sns.gov.pt/monitorizacao-do-sns/servicos-de-urgencia/ 5 https://aqicn.org/map/portugal/pt/ 6 https://mesonet.agron.iastate.edu/request/download.phtml?network=PT__A SOS 7 https://evm.min-saude.pt/#shiny-tab-a_total 4
6
3.2
F. Lino et al.
Hierarchical Structure and Variables
Hospital Level At the hospital level, each record includes information such as the hospital unit (D0), locality, the total number of urgent episodes (M1), and detailed triage information according to the Manchester Triage System (M2 blue, M3 - white, M4 - green, M5 - without triage, M6 - yellow, M7 - orange, M8 - red). Additional variables include average waiting time between triage and medical evaluation, type of emergency service (Basic Emergency Service - SUB, Medical-Surgical Emergency Service - SUMC, Polyvalent Emergency Service SUP and Trauma Center - SUPCT), patient access pathways (self-admission, SNS24, through primary health care warrant, over emergency doctors or hospital doctors referrals), and admitted age ranges. Calendar-related variables, such as day, weekday, month, and holiday indicators, are also included. Furthermore, an operational variable was introduced to indicate whether each ED was open on a given day, as some units occasionally close, resulting in zero recorded visits. Regional Level At the regional level, data include aggregated demand indicators, average waiting times, and environmental variables such as air quality index (AQI) and temperature. Due to limited availability of localized environmental measurements, AQI and temperature were aggregated at the RHA level. The AQI reflects concentrations of major pollutants, including ozone (O3 ), particulate matter (PM2.5 and PM10), nitrogen dioxide (NO2 ), and sulfur dioxide (SO2 ), and was computed according to established guidelines [1]. National Level At the national level, the dataset includes system-wide indicators such as aggregated ED demand, average waiting times, and mortality rates, providing a macro-level view of the healthcare system pressure. The availability of covariates varies across hierarchical levels. Hospital-level data capture local operational and demographic characteristics, whereas regional and national levels provide complementary environmental and health indicators, such as temperature and mortality rates. This creates a heterogeneous feature space across the hierarchy, reflecting real-world data availability and motivating models capable of integrating level-specific information. To ensure dataset completeness over the study period, missing values at specific timestamps were filled using linear interpolation.
4
Methodology
This section presents our Hierarchical Spatio-Temporal Transformer (HierSTT) for multi-step ED demand forecasting across hospital, regional, and national levels. The task consists of jointly predicting patient visits over a 28-day horizon using the previous 42 days, while ensuring hierarchical coherence and effectively leveraging heterogeneous covariates. HierSTT jointly models temporal dynamics and cross-level dependencies through a top-down architecture in which higher-level representations guide lower-level forecasts. A coherence-aware loss further enforces consistency across the hierarchy, enabling the model to capture both global demand trends and local variations.
Hierarchical Spatio-Temporal Transformer for Coherent ED Forecasting
4.1
7
Problem Formulation
We consider a hierarchical forecasting problem defined over a healthcare system organized into three aggregation levels: hospital, regional, and national. At the lowest level, let H = {1, . . . , H} denote the set of hospitals, grouped into a set of regions R = {1, . . . , R}, which together form the national level. For each hospital h ∈ H, we observe a univariate time series of daily ED (h) visits, denoted by yt ∈ R0+ at time t. Aggregated time series at the regional and national levels are defined as the sum of their corresponding lower-level series. Specifically, for each region r ∈ R, and at the national level, ∑ (h) ∑ (h) ∑ (r) (r) (n) (n) yt = yt , y t = yt , and yt = yt . (1) h∈Hr
h∈H
r∈R
In addition to the target series, we consider exogenous covariates that vary across hierarchy levels, with some variables only available at specific levels. Given a fixed-length input sequence of Tin = 42 past observations, the goal is to predict the future values of the time series over a horizon of Tout = 28 days for all hierarchy levels. Formally, for each time t, the model receives as input: (r) (n) t Xt−Tin :t = {x(h) τ , xτ , xτ }τ =t−Tin ,
(2)
(·)
where xτ represents the observed features (including past visits) at each level. The objective is to jointly predict hospital, regional, and national level demand over the forecasting horizon, while ensuring hierarchical coherence, i.e., lower-level forecasts must aggregate consistently to higher-level predictions according to Eq. (1). To address this problem, we formulate hierarchical forecasting as a joint mapping from multi-level historical observations to future predictions across all levels: ( ) (n) (r) (h) ŷt+1:t+Tout , ŷt+1:t+Tout , ŷt+1:t+Tout = fθ (Xt−Tin :t ) , (3) for all hospitals h ∈ H, regions r ∈ R, and the national level, where fθ models temporal dependencies and cross-level interactions. To explicitly model hierarchical dependencies, forecasting is performed in a top-down manner, where national forecasts are first generated, followed by regional forecasts conditioned on national predictions, and hospital forecasts conditioned on regional predictions: ( ) (n) (n) (n) ŷt+1:t+Tout = fθ Xt−Tin :t , (4) ( ) (r) (r) (r) (n) ŷt+1:t+Tout = fθ Xt−Tin :t , ŷt+1:t+Tout , (5) ( ) (h) (h) (h) (r) ŷt+1:t+Tout = fθ Xt−Tin :t , ŷt+1:t+Tout . (6) (n)
The national component fθ is implemented using a TFT [6] to capture global temporal patterns and heterogeneous covariates. Regional and hospital
8
F. Lino et al.
Fig. 2: Overview of the HierSTT. A TFT models national-level ED demand, while spatio-temporal Transformer encoder–decoder modules forecast regional and hospital-level demand. Predictions are generated in a top-down manner: national forecasts condition regional decoding, and regional forecasts condition hospital decoding. (r)
(h)
components, fθ and fθ , use spatio-temporal Transformer encoder–decoder architectures [10] to jointly model temporal dynamics and cross-entity interactions. Higher-level predictions are incorporated into the decoding process, enabling lower-level forecasts to be conditioned on national and regional signals through attention mechanisms. This formulation allows the model to jointly integrate heterogeneous covariates and hierarchical dependencies through top-down conditioning. A schematic overview of the proposed architecture is shown in Fig. 2. 4.2
Temporal Fusion Transformer for National Forecasting (n)
At the national level, the model fθ is implemented as a TFT [6], which is designed for multi-horizon forecasting with heterogeneous covariates. The model processes a window of Tin + Tout days, encoding the past Tin = 42 days through an LSTM and decoding the future Tout = 28 days using a second LSTM conditioned on static context, followed by multi-head self-attention over the decoded sequence. Input covariates include static variables, time-varying real-valued variables observed in the past (e.g., waiting times, ED demands, and mortality counts), and time-varying categorical variables that are known at future time steps (e.g., month, weekday, and holiday indicators). Input variables are processed through Variable Selection Networks (VSNs), which learn time-dependent weights over covariates and suppress less relevant signals. Static covariates are encoded into context vectors that condition both variable selection and subsequent recurrent layers. However, since no natural
Hierarchical Spatio-Temporal Transformer for Coherent ED Forecasting
9
Fig. 3: Spatio-temporal Transformer encoder–decoder used for regional and hospital forecasting. The encoder alternates temporal self-attention within entities and spatial self-attention across entities. The decoder incorporates future-known covariates and higher-level forecasts through temporal, cross-, and spatialattention mechanisms to generate predictions ŷ(ℓ) .
static covariates exist at the national level, a dummy static feature is used for compatibility with the TFT architecture. The final output is a national forecast (n) ŷt+1:t+Tout ∈ RTout , which is then used to condition regional forecasting. 4.3
Spatio-Temporal Transformer for Regional and Hospital Forecasting (r)
(h)
At regional and hospital levels, the functions fθ and fθ are implemented as spatio-temporal Transformer encoderdecoder modules [10]. Unlike the national model, these levels include multiple interacting entities (5 regions and 81 hospitals). The architecture therefore models both temporal dependencies within each entity and spatial dependencies across entities, as shown in Fig. 3. Spatio-Temporal Encoder For each level ℓ ∈ {r, h}, the encoder receives historical features X(ℓ) ∈ RB×Eℓ ×Tin ×Fℓ , where B is the batch size, Eℓ the number of entities, and Fℓ the number of features. Inputs are projected to dimension dmodel , enriched with temporal positional encodings, and processed by L spatiotemporal blocks. Each block alternates between two types of attention. Temporal self-attention is applied independently to each entity over the 42-day history. Spatial selfattention is then applied across all entities at each time step. This allows the encoder to learn both intra-entity temporal patterns and inter-entity correlations. Residual connections, layer normalization, and feed-forward layers are applied throughout.
10
F. Lino et al.
Spatio-Temporal Decoder The decoder predicts the next Tout = 28 days using the encoded memory, future-known covariates, and higher-level forecasts. Decoder queries are built from embeddings of month, weekday, holiday, and open status indicators, together with positional encodings. This yields a query tensor Q ∈ RB×Eℓ ×Tout ×dmodel . Higher-level predictions are projected and added to the decoder queries: national forecasts ŷ(n) ∈ RB×Tout condition regional decoding, while the corresponding regional forecasts ŷ(rh ) ∈ RB×Tout , where rh denotes the region which hospital h belongs, condition hospital decoding. Each decoder layer applies temporal self-attention over the prediction horizon within each entity, cross-attention to the encoder memory, and spatial self-attention across entities. A final linear projections yields the regional and hospital forecasts. 4.4
Hierarchical Coherence-Aware Loss
To encourage hierarchical consistency without hard-constraining the outputs, we add a soft coherence term to the training objective: ( ) L = (1 − α) L(n) + L(r) + L(h) + αLcoh ,
(7)
where, α ∈ [0, 1] balances forecasting accuracy and hierarchical coherence, and each direct loss L(.) is computed using Smooth L1 loss. The coherence loss Lcoh penalizes inconsistencies between predictions across hierarchy levels. Specifically, we impose three aggregation constraints: ( Lcoh = L ŷ
(n)
,
∑
) ŷ
(r)
( + L ŷ
r∈R
(n)
,
∑ h∈H
) ŷ
(h)
( ) R ∑ (h) 1 ∑ (r) L ŷ , + ŷ , R r=1
(8)
h∈Hr
for R = 5 regions. The first two terms enforce national-level consistency with aggregated regional and hospital forecasts, while the third enforces regional consistency with the corresponding hospital aggregates. Importantly, Lcoh compares predictions across levels rather than comparing aggregated predictions with ground truth, avoiding redundancy with the direct forecasting losses.
5
Experimental Setup
The proposed framework is evaluated against statistical, deep learning, and hierarchical reconciliation baselines across hospital, regional, and national forecasting levels. Experiments assess both forecasting accuracy and hierarchical coherence using real-world ED demand data from the Portuguese National Health Service. Results demonstrate that the proposed model consistently outperforms competing approaches on both dimensions. All experiments were implemented in PyTorch and executed on an Intel Xeon Platinum 8260 CPU (2.40 GHz) and an NVIDIA RTX A6000 GPU.
Hierarchical Spatio-Temporal Transformer for Coherent ED Forecasting
11
Dataset and Splits Experiments were conducted on the Portuguese National Health Service ED dataset described in Section 3. Data prior to August 2021 were excluded to reduce the impact of COVID-19-related demand shifts. The remaining data were chronologically split into training, validation, and test sets using cutoff dates of 15 July 2023 and 2 December 2023, yielding 713, 140, and 140 days, respectively. No shuffling was applied to avoid temporal leakage. All models used a 42-day encoder window (Tin ) and a 28-day forecasting horizon (Tout ). This configuration was chosen to capture seasonal trends and support effective resource management in hospital settings. Samples were generated using a sliding-window strategy with one day stride. Data Preprocessing All targets were log-transformed using ỹ = log(1 + y) to stabilize variance across hierarchical levels. Input features were level-specific, comprising 16, 17, and 26 features at the national, regional, and hospital levels, respectively. Features and targets were independently normalized using RobustScaler fitted on the training set. To compute coherence penalties Lcoh in the original count domain, predictions were inverse-scaled and mapped back via ŷ = expm1(·) before evaluating aggregation constraints. 5.1
Implementation Details
The proposed model was instantiated with the following hyperparameters. At the national level, the TFT used a hidden dimension of 128, two LSTM layers, four attention heads, an embedding dimension of 8, and dropout rate of 0.1. The regional and hospital spatio-temporal Transformers used L = 2 blocks, dmodel = 128, nhead = 4, and a dropout rate of 0.2. The coherence weight was set to α = 0.3, which provided the best accuracycoherence trade-off among tested values (detailed results are reported in the supplementary material B). Models were trained for up to 200 epochs using AdamW with 10−4 weight decay, OneCycleLR scheduling with peak learning rate of 3 × 10−4 , gradient clipping at 1.0, and early stopping with 50-epoch patience based on validation loss. The best validation checkpoint was used for testing. For fair comparison, all deep learning experiments were repeated across six fixed random seeds, with results reported as mean ± standard deviation. 5.2
Baselines
We compare the proposed model against statistical, deep learning, and hierarchical reconciliation baselines. Statistical baselines include Naïve [8], ARIMA [4], and ETS [9], applied independently to each series. Non-hierarchical deep learning baselines include LSTM [17], Seq2Seq-LSTM [15], Seq2Seq-GRU [30], Transformer [28], TFT [6], and N-BEATS [23], all trained independently at each hierarchy level without coherence constraints. Finally, using the Nixtla framework [22], we evaluate Bottom-Up, Top-Down, and Middle-Out reconciliation strategies applied to the statistical forecasts. These baselines assess the benefit of jointly learning hierarchical dependencies.
12
F. Lino et al. Regional
National
Lower-left is better ↙
Lower-left is better ↙
Hospital Lower-left is better ↙
700
100
5000 600 4000
60
500
RMSE
RMSE
RMSE
80
400
3000
2000 HierSTT
40
300
↙
HierSTT #
20
40
60
80
100
200 100
1000 200
300
MAE Naive T.D. Naive
ARIMA T.D. ARIMA
400
500
600
↖
HierSTT
1000
2000
MAE ETS T.D. ETS
LSTM M.O. Naive
Seq2seq-LSTM M.O. ARIMA
Seq2seq-GRU M.O. ETS
3000
4000
5000
MAE Transformer B.U. Naive
FA-Transformer B.U. ARIMA
TFT B.U. ETS
N-BEATS HierSTT (Ours)
Fig. 4: Mean MAE vs. mean RMSE for the hospital, regional, and national levels. Error bars denote standard deviation across seeds for deep learning models. T.D. = Top Down; M.O. = Middle Out; and B.U. = Bottom Up.
5.3
Evaluation Metrics
Performance was evaluated in count space after inverting all preprocessing transformations. We report Mean Absolute Error (MAE), Root Mean Squared Error (RMSE), and Weighted Absolute Percentage Error (WAPE): ∑ |ŷt − yt | WAPE = t∑ × 100. (9) t yt MAE and RMSE capture the average magnitude of prediction errors, with RMSE penalizing large errors more strongly. WAPE provides a scale-normalized measure of relative error and is preferred over Mean Absolute Percentage Error (MAPE) because some EDs may close, resulting in zero-demand days for which MAPE is undefined. Hierarchical coherence was evaluated using the Hierarchical Aggregation Error (HAgE), computed as the WAPE between aggregated lower-level forecasts and higher-level predictions (Pred. HAgE). We additionally report HAgE against the observed higher-level ground truth (G.T. HAgE) to assess error propagation across hierarchy levels. Three aggregation settings were considered: hospital-toregional, hospital-to-national, and regional-to-national. Perfectly coherent forecasts achieve Pred. HAgE = 0.
6
Results
We evaluate HierSTT against all baselines across two complementary dimensions: forecasting accuracy and hierarchical coherence. Fig. 4 reports absolute forecasting errors separately for the hospital, regional, and national levels, while Fig. 5 summarizes the trade-off between accuracy and hierarchical coherence. Detailed numerical results for each level are reported in the supplementary material C.
Hierarchical Spatio-Temporal Transformer for Coherent ED Forecasting Predicted Aggregation Coherence
Ground-Truth Aggregation Error Naive ARIMA ETS LSTM Seq2seq-LSTM Seq2seq-GRU Transformer FA-Transformer TFT N-BEATS T.D. Naive T.D. ARIMA T.D. ETS M.O. Naive M.O. ARIMA M.O. ETS B.U. Naive B.U. ARIMA B.U. ETS HierSTT (Ours)
20
Mean G.T. HAgE (%)
Mean Pred. HAgE (%)
20
15
10
5
15
10
5 ↑ HierSTT HierSTT
5
10
15
20
13
25
5
Mean WAPE (%)
10
15
20
25
Mean WAPE (%)
Fig. 5: Trade-off between forecasting accuracy (mean WAPE) and hierarchical coherence (mean HAgE) across all models. The left plot measures coherence between aggregated lower-level and direct higher-level forecasts, while the right compares aggregated predictions against ground-truth aggregates. Results are averaged over hospital→regional, hospital→national, and regional→national transitions.
6.1
Overall Forecasting Performance
Fig. 4 reports MAE and RMSE for each hierarchy level. At the hospital level, HierSTT achieves the lowest MAE, improving by approximately 8% over NBEATS, the second best model. However, it incurs a slightly higher RMSE, suggesting that HierSTT reduces average error at the cost of sensitivity to peak deviations, consistent with the higher variability of hospital demand. Although N-BEATS obtains the lowest hospital-level RMSE, this advantage does not transfer to higher levels. At the regional level, the gains are most pronounced: HierSTT reduces MAE and RMSE by approximately 38% and 27%, respectively, over the Naïve baseline, which is the second best at this level. At the national level, HierSTT employs a TFT as its top-level forecaster, and the results confirm that this component retains its standalone performance when trained within the full hierarchical architecture. Despite the added complexity of jointly modeling two additional hierarchy levels, HierSTT achieves the best RMSE and is essentially tied with TFT on MAE. Across all three levels, HierSTT is the only model that consistently occupies the lower-left region of each plot, while competing models tend to excel at one level at the expense of others. 6.2
Hierarchical Coherence
Coherence arises differently depending on the forecasting strategy. Naïve forecasts are implicitly coherent because they preserve the aggregation structure of historical observations, while post-hoc reconciliation methods enforce coherence
14
F. Lino et al.
through aggregation or disaggregation operations. Independently trained models ignore hierarchical consistency, whereas HierSTT learns coherence directly through its training objective. Fig. 5 reports prediction-side HAgE against mean WAPE (left) and groundtruth HAgE against mean WAPE (right). Naïve and reconciliation-based models achieve zero prediction-side HAgE by construction and are therefore excluded from the left plot. Among unconstrained models, the standard Transformer exhibits the worst prediction-side coherence, while ARIMA achieves the lowest values among statistical models. HierSTT achieves the best prediction-side HAgE overall, improving by approximately 13% over ARIMA. On the right plot, which measures how well aggregated lower-level predictions match ground-truth higherlevel values, HierSTT again outperforms all models, including Naïve and reconciliation baselines, reducing ground-truth HAgE by approximately 18% relative to Seq2seq-LSTM, the second best. This indicates that its lower-level forecasts remain both coherent and accurate when aggregated, a property reconciliation baselines cannot guarantee given the quality of their base forecasts. HierSTT also reduces mean WAPE by over 30% relative to the strongest competing models (N-BEATS, Naïve, and FA-Transformer), placing it clearly in the lower-left corner of both plots. A model can appear competitive in Fig.4 while performing poorly in Fig. 5. TFT is a clear example: strong national-level accuracy but poor hospital-level performance leads to high mean WAPE and weak coherence. This highlights the importance of jointly evaluating accuracy and coherence across all hierarchy levels. HierSTT is the only model that consistently performs strongly on both dimensions, combining the best mean WAPE with learned hierarchical coherence. Error bars further show that it is also the most stable deep learning model across random seeds.
7
Conclusion
In this work, we proposed HierSTT, a hierarchical Transformer-based framework for coherent multi-level forecasting of Emergency Department demand. HierSTT jointly models temporal dynamics, heterogeneous covariates, and cross-level dependencies through top-down conditioning, that allows higher-level forecasts to inform lower-level estimates, and a coherence-aware training objective. We additionally introduced a nationwide Portuguese ED dataset covering 81 hospitals across 5 regional health administrations. Experiments show that HierSTT consistently outperforms statistical, nonhierarchical deep learning, and reconciliation baselines across all hierarchy levels, achieving over 30% average WAPE improvement over the strongest competitor. Unlike post-hoc reconciliation methods or independently trained models, HierSTT jointly achieves high accuracy and near-coherent forecasts. Results also highlight the importance of evaluating hierarchical forecasting using both accuracy and coherence metrics such as HAgE.
Hierarchical Spatio-Temporal Transformer for Coherent ED Forecasting
8
15
Acknowledgments
This work was supported by Fundação para a Ciência e a Tecnologia (FCT) through LARSyS funding (DOIs: 10.54499/LA/P/0083/2020, 10.54499/UIDP/ 50009/2020, and 10.54499/UIDB/50009/2020), and PhD grant 2025.03757.BD (DOI: 10.54499/2025.03757.BD).
References 1. Agency, U.S.E.P.: Technical assistance document for the reporting of daily air quality the air quality index (aqi) 2. Álvarez-Chaves, H., Maseda-Zurdo, I., Muñoz, P., R-Moreno, M.D.: Evaluating the impact of exogenous variables for patients forecasting in an emergency department using attention neural networks. Expert Systems with Applications 240, 122496 (2024) 3. Batal, H., Tench, J., McMillan, S., Adams, J., Mehler, P.S.: Predicting patient visits to an urgent care clinic using calendar variables. Academic Emergency Medicine 8(1), 48–53 (2001) 4. Box, G.E., Jenkins, G.M., Reinsel, G.C., Ljung, G.M.: Time series analysis: forecasting and control. John Wiley & Sons (2015) 5. Boyle, J., Jessup, M., Crilly, J., Green, D., Lind, J., Wallis, M., Miller, P., Fitzgerald, G.: Predicting emergency department admissions. Emergency Medicine Journal 29(5), 358–365 (2012) 6. Caldas, F.M., Soares, C.: A temporal fusion transformer for long-term explainable prediction of emergency department overcrowding. In: Joint European Conference on Machine Learning and Knowledge Discovery in Databases. pp. 71–88. Springer (2022) 7. Findley, D.F., Monsell, B.C., Bell, W.R., Otto, M.C., Chen, B.C.: New capabilities and methods of the x-12-arima seasonal-adjustment program. Journal of Business & Economic Statistics 16(2), 127–152 (1998) 8. Hyndman, R.J., Athanasopoulos, G.: Forecasting: principles and practice. OTexts (2018) 9. Hyndman, R.J., Koehler, A.B., Snyder, R.D., Grose, S.: A state space framework for automatic forecasting using exponential smoothing methods. International Journal of forecasting 18(3), 439–454 (2002) 10. Ji, J., He, J., Lei, M., Wang, M., Tang, W.: Spatio-temporal transformer network for weather forecasting. IEEE Transactions on Big Data 11(2), 372–387 (2024) 11. Kadri, F., Abdennbi, K.: Rnn-based deep-learning approach to forecasting hospital system demands: application to an emergency department. International Journal of Data Science 5(1), 1–25 (2020) 12. Kayacan, E., Ulutas, B., Kaynak, O.: Grey system theory-based models in time series prediction. Expert Systems with Applications 37(2), 1784–1789 (2010) 13. Makridakis, S., Spiliotis, E., Assimakopoulos, V., Semenoglou, A.A., Mulder, G., Nikolopoulos, K.: Statistical, machine learning and deep learning forecasting methods: Comparisons and ways forward. Journal of the Operational Research Society 74, 1–20 (09 2022) 14. Mancuso, P., Piccialli, V., Sudoso, A.M.: A machine learning approach for forecasting hierarchical time series. Expert Systems with Applications 182, 115102 (2021)
16
F. Lino et al.
15. Masood, Z., Gantassi, R., Choi, Y.: A multi-step time-series clustering-based seq2seq lstm learning for a single household electricity load forecasting. Energies 15(7), 2623 (2022) 16. Milner, P.C.: Forecasting the demand on accident and emergency departments in health districts in the trent region. Statistics in Medicine 7(10), 1061–1072 (1988) 17. Ming, C., Lee, G.J., Teo, Y.N., Teo, Y.H., Zhou, X., Ho, E.S., Toh, E.M., Ong, M.E.H., Tan, B.Y., Ho, A.F.: Deep learning modelling to forecast emergency department visits using calendar, meteorological, internet search data and stock market price. Computer methods and programs in biomedicine 267, 108808 (2025) 18. Ming, C., Leung, K.B., Shen, Y., Ho, A.F.: Transformers outperform traditional forecasting models and perform comparably to recurrent neural networks in the prediction of emergency department visits using calendar and meteorological data. Artificial Intelligence in Emergency Medicine 1, 100006 (2026) 19. Moosavi Kashani, S., Zargar Balaye Jame, S., Markazi-Moghaddam, N., Omrani Nava, A.: Comparison of long short-term memory and convolutional neural network models for emergency department patients arrival daily forecasting. Journal of Archives in Military Medicine 12(1), e140888 (2024) 20. Moreno, A., Muñoz, O.: Application of queuing theory to optimize the triage process in a tertiary emergency care (er) department. Journal of Emergencies Trauma and Shock 12, 268–273 (11 2019) 21. Nasios, I., Vogklis, K.: Blending gradient boosted trees and neural networks for point and probabilistic forecasting of hierarchical time series. International Journal of Forecasting 38(4), 14481459 (Oct 2022) 22. Olivares, K.G., Garza, F., Luo, D., Challú, C., Mergenthaler, M., Taieb, S.B., Wickramasuriya, S.L., Dubrawski, A.: Hierarchicalforecast: A reference framework for hierarchical forecasting in python (2023) 23. Oreshkin, B.N., Carpov, D., Chapados, N., Bengio, Y.: N-beats: Neural basis expansion analysis for interpretable time series forecasting. arxiv 2019. arXiv preprint arXiv:1905.10437 (2019) 24. Pulkkinen, E.: Forecasting emergency department arrivals with neural networks. B.S. thesis (2020) 25. Savioli, G., Ceresa, I.F., Gri, N., Bavestrello Piccini, G., Longhitano, Y., Zanza, C., Piccioni, A., Esposito, C., Ricevuti, G., Bressan, M.A.: Emergency department overcrowding: understanding the factors to find corresponding solutions. Journal of personalized medicine 12(2), 279 (2022) 26. Silva, E., Pereira, M.F., Vieira, J.T., Ferreira-Coimbra, J., Henriques, M., Rodrigues, N.F.: Predicting hospital emergency department visits accurately: A systematic review. The International Journal of Health Planning and Management 38(4), 904–917 (2023) 27. Turner, A.J., Anselmi, L., Lau, Y.S., Sutton, M.: The effects of unexpected changes in demand on the performance of emergency departments. Health Economics 29(12), 1744–1763 (2020) 28. Vaswani, A., Shazeer, N., Parmar, N., Uszkoreit, J., Jones, L., Gomez, A.N., Kaiser, Ł., Polosukhin, I.: Attention is all you need. Advances in neural information processing systems 30 (2017) 29. Wang, S.: Neuralreconciler for hierarchical time series forecasting. In: Proceedings of the 17th ACM International Conference on Web Search and Data Mining. p. 731739. WSDM ’24, Association for Computing Machinery, New York, NY, USA (2024)
Hierarchical Spatio-Temporal Transformer for Coherent ED Forecasting
17
30. Xu, J., Wang, K., Lin, C., Xiao, L., Huang, X., Zhang, Y.: Fm-gru: A time series prediction method for water quality based on seq2seq framework. Water 13(8), 1031 (2021) 31. Yule, G.U.: On a method of investigating periodicities in disturbed series with special reference to wolfers sunspot numbers. Statistical Papers of George Udny Yule pp. 389–420 (1971) 32. Zhang, Q., Wang, K., Guo, H., Yang, S., Li, C.: Base on ets model for forcasting emergency department visits. In: 2018 IEEE 3rd Advanced Information Technology, Electronic and Automation Control Conference (IAEAC). pp. 2148–2151. IEEE (2018)
18
F. Lino et al.
A
Analysis of the Hierarchical Portuguese ED Dataset
A.1
Variable Availability Across Hierarchical Levels
Table 1: Available variables at each hierarchical level of the dataset of ED activity in Portugal. Variable Type Hospital ID RHA Locality ED Visits (M1) Triage (M2-M8) Waiting Time Calendar (weekday/month) Holidays Air Quality (AQI) Temperature Mortality Open ED Indicator ED Type / Access / Age
Hospital Level Regional Level National Level 3 3 3 3 3 7 3 3 7 7 7 3 3
7 3 7 3 3 3 3 3 3 3 7 3 7
7 7 7 3 3 3 3 3 7 7 3 3 7
The forecasting dataset integrates information from the hospital, regional, and national levels of the Portuguese healthcare system. Since variable availability differs across levels, the resulting feature space is heterogeneous across the hierarchy. Table 1 summarizes the variables available at each level. Hospital level variables includes local operational and demographic information, such as hospital identifier, Regional Health Administration (RHA), locality, emergency department (ED) visits, triage indicators, calendar variables, holidays, open/closed status, ED type, patient access pathways, and admitted age ranges. Regional level variables include aggregated ED demand and triage indicators, waiting times, calendar variables, holidays, open/closed status, and environmental variables, including air quality index (AQI) and temperature. Environmental variables are included only at this level due to the limited availability of reliable hospital-specific measurements. National level variables include aggregated ED visits, triage indicators, waiting time, calendar variables, holidays, open/closed status, and mortality, providing a system-wide view of healthcare demand. This level-dependent variable availability reflects real-world data collection constraints and motivates a hierarchical forecasting framework capable of combining heterogeneous level-specific information while producing coherent predictions across the hierarchy.
Hierarchical Spatio-Temporal Transformer for Coherent ED Forecasting
Temporal Patterns in Hospital-Level ED Visits
COVID-19
450
Val split
Region 1: Hospital-level M1 time series - full timeframe 500
Val split
COVID-19
Train split
Region 0: Hospital-level M1 time series - full timeframe 300
Train split
A.2
19
250 400 350 Region 0 hospitals H2 H5 H6 H14 H15 H20 H35 H57 H64 H77
150
100
300 M1
M1
200
Region 1 hospitals H0 H9 H29 H45 H58
250 200 150 100
50 50
2021-07-01
2022-01-01
2022-07-01
2023-01-01 Date
2023-07-01
2024-01-01
0 2021-01-01
2024-07-01
2021-07-01
2022-01-01
2022-07-01
(a)
2024-01-01
2024-07-01
800
M1
400
300
200
Region 3 hospitals H3 H12 H17 H25 H27 H32 H34 H38 H39 H42 H43 H44 H46 H49 H61 H63 H66 H67 H68 H73 H80
700 600 500 M1
Region 2 hospitals H1 H18 H23 H26 H28 H30 H33 H40 H47 H50 H51 H53 H55 H59 H70 H76 H78 H79
500
Train split
COVID-19
Val split
Train split
COVID-19
Region 3: Hospital-level M1 time series - full timeframe 900
600
400 300 200
100
0 2021-01-01
2023-07-01
(b)
Region 2: Hospital-level M1 time series - full timeframe 700
2023-01-01 Date
Val split
0 2021-01-01
100
2021-07-01
2022-01-01
2022-07-01
2023-01-01 Date
2023-07-01
2024-01-01
2024-07-01
0 2021-01-01
2021-07-01
2022-01-01
2022-07-01
(c)
2023-01-01 Date
2023-07-01
2024-01-01
2024-07-01
(d) Val split
COVID-19
900
Train split
Region 4: Hospital-level M1 time series - full timeframe 1000
Region 4 hospitals H4 H7 H8 H10 H11 H13 H16 H19 H21 H22 H24 H31 H36 H37 H41 H48 H52 H54 H56 H60 H62 H65 H69 H71 H72 H74 H75
800 700
M1
600 500 400 300 200 100 0 2021-01-01
2021-07-01
2022-01-01
2022-07-01
2023-01-01 Date
2023-07-01
2024-01-01
2024-07-01
(e)
Fig. 6: Daily hospital-level M1 time series within each region acrosss the full available time frame. Each line represents one hospital. The plots include the COVID-19 affected period and the subsequent modeling period.
To contextualize the forecasting task, Fig. 6 shows the hospital-level M1 time series across the five regions over the full available time frame, including the COVID-19 period and the subsequent modeling period. The plots highlight the variability across hospitals, both in demand magnitude and temporal dynamics. This heterogeneity is particularly evident in Regions 2–4, where hospitals operate at different demand scales within the same regional unit. The figures also show atypical demand patterns during the COVID-19 period, motivating the use of data from 1 August 2021 onward for forecasting experiments. The modeling start date and train/validation splits are indicated to illustrate the temporal evaluation protocol.
F. Lino et al. Region 0: Two-weeks M1 zoom (2023-12-03 to 2023-12-16) 240
Val split
Month start
Train split
Region 0: Hospital-level M1 time series - modelling period 300
Saturday
20
220
250 200 180 Region 0 hospitals H2 H5 H6 H14 H15 H20 H35 H57 H64 H77
150
100
Region 0 hospitals H2 H5 H6 H14 H15 H20 H35 H57 H64 H77
160 M1
M1
200
140 120 100 80
50 60 0 2021-07-01
2022-01-01
2022-07-01
2023-01-01 Date
2023-07-01
2024-01-01
40 2023-12-03
2024-07-01
2023-12-06
2023-12-09
2023-12-12
(a)
Region 1: Two-weeks M1 zoom (2023-12-03 to 2023-12-16) Saturday
450
Val split
Train split
Month start
450
2023-12-18
(b)
Region 1: Hospital-level M1 time series - modelling period 500
2023-12-15
Date
400
400 350 350
M1
Region 1 hospitals H0 H9 H29 H45 H58
250 200
M1
300
300
Region 1 hospitals H0 H9 H29 H45 H58
250
200
150 150 100 100
50 0 2021-07-01
2022-01-01
2022-07-01
2023-01-01 Date
2023-07-01
2024-01-01
2024-07-01
50 2023-12-03
2023-12-06
(c)
2023-12-09
2023-12-12
2023-12-15
2023-12-18
Date
(d)
Fig. 7: Hospital-level M1 series over the modeling period for Region 0 (a) and Region 1 (c), selected as illustrative examples due to their smaller number of hospitals. Panels (b) and (d) show two-week zoomed examples of hospital-level trajectories for the same regions. Vertical dashed lines indicate Saturdays.
Fig. 7 presents modeling-period and short-term examples for Regions 0 and 1, selected because their smaller number of hospitals facilitates visualization. The modeling-period plots reveal persistent differences between hospitals within the same region, while two-week zoom plots (Fig. 7b and 7d) illustrate the short-term variability and recurring weekly patterns. The Saturday markers further highlight the importance of calendar structure, supporting the inclusion of futureknown temporal covariates. A.3
Sliding Window Sample Construction
Forecasting samples were generated using a sliding-window strategy with stride one. Each sample consists of a 42-day encoder window, Tin = 42, followed by a 28-day forecasting horizon, Tout = 28, as illustrated in Fig. 8. This input–output configuration was selected to capture short-term and weekly seasonal patterns while producing forecasts over a horizon relevant for operational planning in hospital emergency departments. The chronologically ordered time series were divided into training, validation, and test sets using cutoff dates of 15 July 2023 and 2 December 2023, resulting in 713 training days, 140 validation days, and 140 test days. No shuffling was applied, preserving the temporal order of the observations and preventing information leakage between splits.
Hierarchical Spatio-Temporal Transformer for Coherent ED Forecasting
21
Fig. 8: Sliding window strategy used to generate samples. Each sample consists of a 42-day encoder window (Tin ) followed by a 28-day prediction window (Tout ). Consecutive samples are shifted by one day.
Given a split with N consecutive days, the number of valid samples is computed as N − Tin − Tout + 1. Therefore, the 713 training days yield 644 training samples, while the 140-day validation and test periods each yield 71 samples.
B
Ablation study: Effect of the Coherence Weight α
B.1
Loss Imbalance and the Role of α
The hierarchical loss defined in (7) supervises fundamentally different numbers of time series across its two terms. The accuracy term, L(n) +L(r) +L(h) , operates over 87 individual series (81 hospital, 5 regional, and 1 national), whereas the coherence term Lcoh operates over 7 aggregation constrains: hospital-to-regional (5), hospital-to-national (1), and regional-to-national (1). This asymmetry has a direct consequence on the effective influence of α. The per-constraint contribution of the coherence term is α/7, while the per-series contribution of the accuracy term is (1 − α)/87. The effective ratio between them is therefore: α · 87 , (10) (1 − α) · 7 which departs from the nominal ratio α/(1 − α) implied by α alone. At α = 0.3, the coherence term already exerts approximately five times more influence per constraint than the accuracy term does per series, meaning that coherence strongly shapes the learned representations even at low α values. Conversely, even at α = 0.999, the effective ratio is approximately 12:1 rather than 999:1, confirming that the accuracy term is never truly suppressed and retains meaningful gradient signal throughout training. These considerations informed our choice of α and highlight that its interpretation must account for the structural imbalance between the two loss components.
22
F. Lino et al.
Table 2: Performance of HierSTT for different values of the coherence weight α. HAgE denotes Hierarchical Aggregation Error computed using WAPE. WAPE [%] Pred. HAgE [%] G.T. HAgE [%] Level MAE RMSE Per Level Avg. Aggregation Per Agg. Avg. Per Agg. Avg. Hospital 22.16 ± 0.42 33.5 ± 0.7 10.87 ± 0.20 Hosp.→Reg. 2.34 ± 0.74 5.56 ± 0.27 0 Regional 162.10 ± 4.92 228.53 ± 4.96 4.90 ± 0.15 7.13 ± 0.28 Hosp.→Nat. 3.28 ± 0.65 3.04 ± 0.61 4.95 ± 0.33 4.94 ± 0.25 National 929.89 ± 133.19 1110.10 ± 173.92 5.63 ± 0.79 Reg.→Nat. 3.49 ± 0.86 4.32 ± 0.25 Hospital 22.03 ± 0.18 33.22 ± 0.37 10.81 ± 0.09 Hosp.→Reg. 2.03 ± 0.42 5.56 ± 0.16 0.3 Regional 162.77 ± 4.55 229.70 ± 5.99 4.92 ± 0.14 7.06 ± 0.20 Hosp.→Nat. 3.11 ± 0.65 2.83 ± 0.53 4.93 ± 0.19 4.93 ± 0.14 National 897.78 ± 114.69 1064.59 ± 145.65 5.44 ± 0.67 Reg.→Nat. 3.34 ± 0.65 4.30 ± 0.17 Hospital 22.32 ± 0.31 33.75 ± 0.74 10.94 ± 0.16 Hosp.→Reg. 2.10 ± 0.62 5.45 ± 0.25 0.5 Regional 164.73 ± 5.05 230.73 ± 4.98 4.98 ± 0.16 7.22 ± 0.17 Hosp.→Nat. 3.20 ± 0.71 2.93 ± 0.40 4.84 ± 0.33 4.90 ± 0.23 National 942.28 ± 101.66 1126.85 ± 117.78 5.72 ± 0.61 Reg.→Nat. 3.49 ± 0.48 4.43 ± 0.21 Hospital 22.18 ± 0.26 33.42 ± 0.49 10.88 ± 0.12 Hosp.→Reg. 1.91 ± 0.50 5.53 ± 0.23 0.7 Regional 164.46 ± 6.70 231.54 ± 7.60 4.97 ± 0.21 7.19 ± 0.32 Hosp.→Nat. 3.38 ± 0.87 2.91 ± 0.54 4.96 ± 0.28 4.97 ± 0.22 National 941.08 ± 145.79 1131.79 ± 206.36 5.71 ± 0.90 Reg.→Nat. 3.43 ± 0.83 4.42 ± 0.27 Hospital 22.34 ± 0.47 34.23 ± 0.87 10.96 ± 0.23 Hosp.→Reg. 1.24 ± 0.44 5.72 ± 0.13 0.99 Regional 186.68 ± 10.05 260.83 ± 13.53 5.66 ± 0.31 7.46 ± 0.22 Hosp.→Nat. 1.79 ± 0.80 1.46 ± 0.64 5.00 ± 0.18 5.26 ± 0.16 National 948.26 ± 72.06 1134.55 ± 89.38 5.76 ± 0.43 Reg.→Nat. 1.35 ± 0.70 5.08 ± 0.31 Hospital 24.15 ± 0.61 37.07 ± 0.93 11.86 ± 0.31 Hosp.→Reg. 0.50 ± 0.16 6.23 ± 0.25 0.999 Regional 194.12 ± 7.05 276.42 ± 14.60 5.90 ± 0.22 7.74 ± 0.23 Hosp.→Nat. 0.73 ± 0.16 0.65 ± 0.14 5.41 ± 0.20 5.66 ± 0.20 National 899.54 ± 29.60 1076.37 ± 47.72 5.47 ± 0.19 Reg.→Nat. 0.72 ± 0.17 5.33 ± 0.16 α
B.2
Results Across α values
Table 2 reports the performance of HierSTT across different values of α, averaged over six fixed seeds (mean ± standard deviation). Each configuration is evaluated on forecasting accuracy (MAE, RMSE, WAPE per level and average WAPE) and hierarchical coherence (prediction-side HAgE and ground-truth HAgE across the three aggregation transitions). Several trends emerge. First, α = 0.3 achieves the lowest average WAPE (7.06±0.20%) and the best hospital and national level MAE, RMSE, and WAPE, making it the best overall operating point. Second, prediction-side HAgE decreases as α increases, from 3.04% at α = 0 to 0.65% at α = 0.999, confirming that the coherence objective progressively suppresses cross-level inconsistencies. Third, the accuracy-coherence trade-off remains mild up to α = 0.5 (average WAPE of 7.22%), but deteriorates for larger values: at α = 0.99, regional WAPE increases to 5.66% and at α = 0.999, hospital WAPE rises to 11.86%. This behavior is consistent with imbalance described in the previous subsection, where larger α values increase the relative influence of the 7 aggregation constraints with respect to the 87 directly supervised time series. Notably, even at α = 0.999, the model does not collapse to a degenerate solution. National-level WAPE remains at 5.47%, only slightly above the best value of 5.44% obtained at α = 0.3. This behavior is also consistent with the loss-imbalance analysis above: although α = 0.999 appears to almost entirely suppress the accuracy term in nominal weight, the effective per-series contribution of the accuracy loss remains meaningful once the imbalance between the 87 supervised time series and the 7 aggregation constraints is taken into account. Consequently, the model continues to receive direct forecasting supervision, while the coherence term strongly reduces cross-level inconsistencies.
Hierarchical Spatio-Temporal Transformer for Coherent ED Forecasting
23
Table 3: Forecasting performance across all models and hierarchy levels. Green = best, underline = second best per column. FA-Transformer stands for FutureAware Transformer, since the model uses future covariates in the decoder. Hospital
Hierarchical
Non-Hierarchical
Model Naïve ARIMA ETS LSTM Seq2seq-LSTM Seq2seq-GRU Transformer FA-Transformer TFT N-BEATS
MAE
RMSE
Regional WAPE [%]
MAE
RMSE
National WAPE [%]
MAE
RMSE
WAPE [%]
25.5 31.6 15.4 261.6 313.0 8.1 1251.1 1461.7 7.6 28.1 33.8 15.6 345.3 411.0 10.4 1721.7 2034.5 10.5 29.9 35.9 16.5 362.0 431.1 11.2 1879.4 2202.7 11.5 26.4 ± 2.7 32.7 ± 3.0 14.7 ± 1.0 359.9 ± 117.0 423.2 ± 127.9 12.4 ± 3.1 1360.4 ± 188.5 1708.2 ± 215.0 8.2 ± 1.1 24.5 ± 1.1 30.7 ± 1.2 14.1 ± 0.7 299.2 ± 132.9 357.9 ± 152.6 13.4 ± 8.4 1054.2 ± 235.2 1312.5 ± 281.3 6.4 ± 1.4 25.1 ± 2.8 31.2 ± 3.4 14.4 ± 1.0 345.9 ± 40.4 422.8 ± 50.3 15.0 ± 6.6 948.2 ± 102.9 1187.7 ± 128.2 5.7 ± 0.6 28.1 ± 1.4 34.8 ± 1.5 15.8 ± 0.4 498.6 ± 74.0 605.2 ± 82.4 15.7 ± 3.0 3301.9 ± 614.3 4358.4 ± 738.4 20.0 ± 3.7 26.9 ± 0.9 33.3 ± 1.0 14.7 ± 0.4 307.0 ± 80.8 354.4 ± 86.7 10.3 ± 2.7 1020.4 ± 435.1 1227.6 ± 467.4 6.2 ± 2.7 82.0 ± 10.8 89.8 ± 10.5 46.2 ± 10.0 381.3 ± 28.7 455.0 ± 29.6 12.1 ± 1.7 896.3 ± 119.3 1072.0 ± 134.0 5.4 ± 0.7 23.8 ± 0.9 29.7 ± 1.0 14.1 ± 0.6 330.9 ± 31.3 408.3 ± 35.4 10.6 ± 1.3 1149.8 ± 166.5 1402.1 ± 200.8 7.0 ± 1.0
Naïve 25.5 31.6 15.4 Top-Down ARIMA 30.0 35.9 16.4 ETS 32.7 39.1 17.7 Naïve 25.5 31.6 15.4 Middle-Out ARIMA 29.2 35.0 16.1 ETS 31.3 37.6 17.1 Naïve 25.5 31.6 15.4 Bottom-Up ARIMA 28.1 33.8 15.6 ETS 29.9 35.9 16.5 HierSTT (ours) 22.0 ± 0.2 33.2 ± 0.4 10.8 ± 0.1
261.6 363.1 393.7 261.6 345.3 362.0 261.6 321.0 335.9 162.8 ± 4.6
313.0 430.2 466.4 313.0 411.0 431.1 313.0 381.7 398.3 229.7 ± 6.0
8.1 10.9 12.1 8.1 10.4 11.2 8.1 9.8 10.2 4.9 ± 0.1
1251.1 1461.7 7.6 1721.7 2034.5 10.5 1879.4 2202.7 11.5 1251.1 1461.7 7.6 1612.6 1910.2 9.8 1706.1 2013.2 10.4 1251.1 1461.7 7.6 1549.9 1830.1 9.4 1629.7 1917.2 9.9 897.8 ± 114.7 1064.6 ± 145.7 5.4 ± 0.7
Based on these results, α = 0.3 is selected for all main experiments, as it provides the best trade-off between forecasting accuracy and hierarchical coherence.
C
Extended Experimental Results
C.1
Overall Forecasting Performance
Table 3 reports MAE, RMSE and WAPE across all hierarchy levels. Classical statistical baselines remain competitive, particularly Naïve forecasting, but are generally outperformed by deep learning models at the hospital level, where ED demand exhibits greater variability and complexity, and at national level. Among the hierarchical reconciliation approaches (top-down, middle-out, bottom-up), reconciliation rarely improves over the corresponding non-hierarchical baselines. The Nïve model yields identical results across all strategies because repeating the last observed 28-day preserves scale consistency by construction. The only exception is bottom-up reconciliation applied to ARIMA and ETS, which improves regional and national performance, suggesting that aggregating sufficiently structured lower-level forecasts can propagate useful signal upward. Overall, these results indicate that hierarchy alone is insufficient and effective cross-level integration during training is essential. Non-hierarchical deep learning models generally improve forecasting accuracy, although performance varies substantially across hierarchy levels. TFT achieves strong regional and national WAPE (12.1% and 5.4%, respectively) but performs poorly at the hospital level (WAPE 46.2%), highlighting the limitations of single-level optimization. N-BEATS achieves the best hospital-level
24
F. Lino et al.
Table 4: HAgE across aggregation transitions. Pred. measures inconsistency between predictions at different levels, while G.T. measures the discrepancy between aggregated lower-level predictions and observed higher-level values. ‘—’ indicate reconciliation baselines where prediction-side HAgE is zero by construction. Green = best, underline = second best. Hosp. → Reg.
Hierarchical
Non-Hierarchical
Model Naïve ARIMA ETS LSTM Seq2seq-LSTM Seq2seq-GRU Transformer FA-Transformer TFT N-BEATS
Hosp. → Nat.
Reg. → Nat.
Pred. [%] G.T. [%] Pred. [%] G.T. [%] Pred. [%] G.T. [%] — 8.0 — 7.6 — 7.6 3.7 9.7 3.9 9.4 2.1 9.8 4.1 10.2 4.8 9.9 3.0 10.4 8.7 ± 4.6 8.2 ± 1.6 10.1 ± 2.2 6.7 ± 1.2 11.6 ± 4.8 8.6 ± 3.1 7.0 ± 3.5 6.5 ± 0.7 6.0 ± 0.7 4.9 ± 0.3 6.1 ± 1.6 6.5 ± 2.3 7.8 ± 1.6 7.2 ± 1.7 6.0 ± 1.3 5.9 ± 1.2 7.5 ± 1.4 8.6 ± 0.4 14.6 ± 1.1 8.9 ± 1.0 20.8 ± 3.1 8.2 ± 1.4 18.4 ± 2.2 9.6 ± 1.2 6.9 ± 1.2 8.7 ± 0.3 9.9 ± 4.1 8.1 ± 0.5 9.2 ± 4.0 8.1 ± 2.4 20.8 ± 4.1 25.4 ± 2.9 22.7 ± 4.5 22.3 ± 3.7 7.8 ± 0.8 8.6 ± 0.5 8.1 ± 1.2 6.2 ± 0.4 8.1 ± 2.0 5.5 ± 0.4 10.9 ± 1.8 8.7 ± 0.6
Naïve — 8.0 — 7.6 — 7.6 Top-Down ARIMA — 11.0 — 10.5 — 10.5 ETS — 12.0 — 11.5 — 11.5 — 8.0 — 7.6 — 7.6 Naïve — 10.5 — 9.8 — 9.8 Middle-Out ARIMA ETS — 11.0 — 10.4 — 10.4 Naïve — 8.0 — 7.6 — 7.6 Bottom-Up ARIMA — 9.7 — 9.4 — 9.4 ETS — 10.2 — 9.2 — 9.2 HierSTT (ours) 2.0 ± 0.4 5.6 ± 0.2 3.1 ± 0.7 4.9 ± 0.2 3.3 ± 0.7 4.3 ± 0.2
RMSE (29.7) and TFT attains the best national level MAE (896.3). Seq2seq architectures perform consistently across levels, though without leading at any. HierSTT outperforms all baselines at the hospital and regional levels, reducing hospital WAPE by 23% relative to N-BEATS (10.8% vs. 14.1%) and regional MAE by 38% relative to the best non-hierarchical competitor (162.8 vs. 261.6). At the national level, it remains competitive with TFT. Notably, the TFT component of HierSTT preserve its standalone national level performance when trained jointly across all levels (MAE 897.8 vs. 896.3), showing that end-toend hierarchical training does not compromise higher-level forecasting. Overall, these results suggest that jointly modeling spatial dependencies and hierarchical context improves performance across all aggregation levels without sacrificing accuracy at any single level. C.2
Hierarchical Coherence Analysis
Table 4 reports HAgE across the three aggregation transitions. Prediction-side HAgE measures inconsistencies between forecasts at different hierarchy levels. Naïve forecasts are implicitly coherent because repeating historical observations preserves the aggregation structure, while reconciliation baselines achieve zero prediction-side HAgE through explicit post-hoc aggregation constraints. These values are therefore marked with ‘—’, since coherence is algebraically enforced rather than learned.
Hierarchical Spatio-Temporal Transformer for Coherent ED Forecasting
25
×104
1.9
Saturday
National - Sample 32 2
National ground-truth National prediction Sum of regional predictions Sum of hospital predictions
1.8
M1
1.7
1.6
1.5
1.4
1.3
0
5
10
15
20
25
30
Prediction day
Fig. 9: Example of the national-level forecast for the test sample with the lowest national WAPE. The figure compares the ground truth, the direct national prediction, the sum of regional predictions, and the sum of hospital predictions over the 28-day forecasting horizon. Vertical lines indicate Saturdays. Among statistical and deep learning models, the standard Transformer exhibits the worst coherence (14.6%, 20.8%, 18.4%), whereas ARIMA achieves the lowest prediction-side HAgE (3.7%, 3.9%, 2.1%), reflecting its smoother forecasts. In contrast, HierSTT learns hierarchical consistency directly through Lcoh , achieving prediction-side HAgE of 2.0%, 3.1%, and 3.3%, corresponding to the best results for hospital-to-regional and hospital-to-national aggregation, and second best for regional-to-national. Ground-truth HAgE measures the discrepancy between aggregated lowerlevel predictions and observed higher-level values. HierSTT also achieves the lowest errors across all transitions (5.6%, 4.9%, 4.3%), indicating that its lowerlevel forecasts remain both coherent and accurate when aggregated, unlike reconciliation baselines that enforce coherence but remain limited by the quality of their base forecasts. C.3
Qualitative Forecast Results
To complement the quantitative evaluation, we provide qualitative forecast examples generated by HierSTT in Fig. 9, 10, and 11. Examples were selected based on the lowest WAPE at each hierarchy level, allowing inspection of both forecasting accuracy and consistency between direct predictions and hierarchical aggregations. The added vertical Saturday markers further highlight the weekly structure present across hierarchy levels, supporting the inclusion of calendarrelated covariates. At the national and regional levels (Fig. 9 and 10), the model captures the weekly demand structure and the main rises and drops across the 28-day horizon. The aggregated lower-level predictions closely follow the direct forecasts,
F. Lino et al. Region 0 - Sample 44
1150
950
3600 3400
900 M1
1050 1000
3200
850
M1
1100 M1
Saturday
Saturday
1200
Region 2 - Sample 40
Region 1 - Sample 9 1000
Saturday
26
3000
800
950
2800
900
750
850
700
2600
800 10 15 20 Prediction day Region 3 - Sample 44
25
30
0
6500
5
10 15 20 Prediction day Region 4 - Sample 45
25
30
Saturday
5
Saturday
0
7500
2400 0
5
10
15 20 Prediction day
25
30
Regional ground-truth Regional prediction Sum of hospital predictions
7000 6000 M1
M1
6500 5500
6000
5000
5500
4500 0
5
10
15 20 Prediction day
25
30
5000 0
5
10
15 20 Prediction day
25
30
Fig. 10: Regional-level forecasts for the test sample with the lowest WAPE independently selected for each region. Each subplot shows the regional ground truth, the direct regional prediction, and the sum of the corresponding hospitallevel predictions. Vertical lines indicate Saturdays.
indicating strong hierarchical coherence throughout the forecasting window. In particular, the national example (Fig. 11) achieves a WAPE of 3.48%, while remaining highly consistent with both aggregated regional and hospital predictions. Across the five regional examples (Fig. 10), the aggregated hospital predictions closely follow the regional forecasts. Regional differences may partly reflect the hierarchical composition of each region. Regions with more hospitals aggregate a larger number of local demand signals, smoothing hospital-specific variability and facilitating coherent regional prediction. In contrast, regions with fewer hospitals, such as Region 0 and 1, are more sensitive to individual hospital fluctuations. At the hospital level (Fig. 11), the model captures local temporal dynamics despite the higher variability of individual hospital demand. Some sharp peaks are underestimated, particularly in lower-volume hospitals (Fig. 11a and Fig. 11b), but the main temporal patterns are preserved. We additionally include Hospital H44 (Fig. 11f) as an operational example containing hospital closure days. The model predicts near-zero demand during these periods, illustrating the benefit of including the hospital open/closed status as a future-known decoder variable.
Hierarchical Spatio-Temporal Transformer for Coherent ED Forecasting
Hospital H58 - Sample 43 | WAPE = 5.28%
240
Hospital ground-truth Hospital prediction Saturday
340
230
Hospital ground-truth Hospital prediction Saturday
320
220
300
210
M1
M1
Saturday
Saturday
Hospital H35 - Sample 57 | WAPE = 7.74% 250
27
200
280
190
260
180 240
170
220
160 0
5
10
15 Prediction day
20
25
30
0
5
10
(a)
25
30
Hospital H67 - Sample 39 | WAPE = 4.59% Saturday
Saturday
440 420
380
400 380 M1
360 M1
20
(b)
Hospital H76 - Sample 21 | WAPE = 4.10% 400
15 Prediction day
340
360 340
Hospital ground-truth Hospital prediction Saturday
320
320 300
300
Hospital ground-truth Hospital prediction Saturday
280 280 0
5
10
15 Prediction day
20
25
30
0
5
10
(c)
20
25
30
(d) Hospital H44 - Sample 19 | WAPE = 16.28%
45
Hospital ground-truth Hospital prediction Saturday
Saturday
Saturday
Hospital H69 - Sample 30 | WAPE = 4.25% 850
15 Prediction day
40 35
800
30 25
M1
M1
750
700
Hospital ground-truth Hospital prediction Saturday
20 15 10
650
5 600
0 0
5
10
15 Prediction day
(e)
20
25
30
0
5
10
15 Prediction day
20
25
30
(f)
Fig. 11: Hospital-level forecasts for representative hospitals selected from each region. For each region, the hospital example with the lowest WAPE was selected ((a) to (e)), and an additional example is shown for Hospital H44 (f), which includes days in which the hospital is closed. Vertical lines indicate Saturdays.