Hybrid Statistical–Metaheuristic Inventory Modeling: Integrating SARIMAX with Skew-Normal and Zero-Inflated Errors in Clinical Laboratory Demand Forecasting

1. Introduction

Clinical laboratories are essential to public health, supporting medical diagnosis as well as clinical and pharmacological monitoring—core functions that ensure quality healthcare delivery and effective disease control [1]. However, monthly demand for diagnostic analytes is highly uncertain and often follows cyclic, seasonal, and intermittent patterns [2]. These challenges are further intensified in public sector settings, where procurement cycles are lengthy, budgets inflexible, and packaging regulations limit supply responsiveness [3].

Time series models such as ARIMA and SARIMAX are commonly used for demand forecasting in healthcare and logistics [4]. Nevertheless, these models typically assume homoscedastic and symmetric Gaussian residuals. This assumption is often violated in clinical demand series, which exhibit non-normal behaviors such as skewness, high variance, and zero-inflation due to sporadic usage patterns [5]. To address this, the use of skew-normal distributions has been proposed to capture residual asymmetry [6], while zero-inflated variants help account for excess zero observations [7].

In parallel, metaheuristic optimization techniques such as Particle Swarm Optimization (PSO), Genetic Algorithms (GA), and Ant Colony Optimization (ACO) offer flexible strategies for solving nonlinear and constrained inventory problems [8]. These methods support adaptive parameter tuning and are well suited for multi-stage stochastic planning environments [9]. However, they are rarely integrated with domain-specific statistical models for forecast-driven inventory control in clinical applications.

Recent interest in Explainable Artificial Intelligence (XAI) highlights the need for transparency in hybrid modeling approaches [10]. While deep learning models can offer high predictive accuracy, their black-box nature limits their interpretability in high-stakes domains like healthcare. In contrast, SARIMAX models with structured, non-Gaussian residuals can provide both interpretability and performance [11].

Despite these advances, the literature lacks integrated frameworks that combine (i) skew-aware and zero-sensitive time series forecasting, (ii) explainable optimization methods under realistic constraints such as packaging multiples and budget limits, and (iii) full traceability from statistical residuals to procurement decisions.

This paper addresses these gaps by proposing a hybrid, explainable framework that integrates SARIMAX models with Skew-Normal and Zero-Inflated Skew-Normal residuals, coupled with a constrained inventory optimization layer based on the PSO. The method is validated on real clinical consumption data from a Chilean hospital and explicitly incorporates institutional constraints such as packaging formats and a fixed procurement budget. The framework is fully reproducible and interpretable, making it suitable for operational deployment in clinical settings.

The remainder of the paper is structured as follows. Section 2 presents the methodology, including data description, forecasting model specification, residual diagnostics, and the optimization strategy using PSO. Section 3 reports the main results, focusing on forecast accuracy, residual characteristics, and inventory cost savings. Section 4 discusses the implications of the findings in terms of explainability and operational value. Finally, Section 5 concludes with a summary of contributions and suggestions for future research.

2. Methodology

This section describes the modelling and optimisation framework used in this study, organised in five sequential stages: (i) data preprocessing, (ii) demand forecasting, (iii) model validation and explainability, (iv) inventory cost optimisation, and (v) implementation via PSO.

2.1. Data Description and Parameters

We analyzed monthly consumption data for 34 items over a three-year period (January 2021–December 2023) from the Clinical Chemistry section of a public hospital in Chile reagents, covering both high-frequency and low-frequency determinations. This study was approved by the Ethics Committee of the Faculty of Chemistry and Pharmacy, University of Valparaíso (Approval Code: CG-05-2024). Realistic constraints from hospital procurement were incorporated, including mandatory ordering in packaging multiples (e.g., boxes of 12 or 24 units).

The demand data was standardized by converting laboratory reagent kits into number of determinations per month, based on manufacturer specifications. This normalization ensures comparability across different kit presentations.

The following cost parameters were estimated based on hospital records:

• Ordering cost ($o_i$) – average cost of issuing a purchase order, including administrative labor and shipping.
• Holding cost ($h_i$) – annual cost of storing one determination unit, considering energy, losses, and space.
• Shortage cost ($s_i$) – cost associated with stockouts, estimated from rescheduling procedures and patient re-attendance.
• Procurement cost ($c_i$) – average unit cost of each determination, computed from historical purchasing prices.

All monetary values were adjusted to 2023 Chilean pesos using the IPC inflation index. An initial inventory level $I_0$ was defined based on stock records from December 2023.

Demand distributions exhibited high variability, non-Gaussian skewness, and frequent zero values, justifying the use of structured residual models.

2.2. Forecasting Models

To model the irregular and asymmetric nature of clinical laboratory reagent demand, we implemented SARIMAX models with structured residuals:

• Skew-Normal (SN) residuals, capturing asymmetry;
• Zero-Inflated Skew-Normal (ZISN) residuals, capturing both asymmetry and excess zeros.

The general SARIMAX model is defined as:

$$Y_t={\mathbf{X}}_t^{\prime }\beta +{\mathbf{Z}}_t^{\prime }\gamma +{\epsilon }_t$$

(1)

where:

-

$Y_t$ is the observed demand at time t,

-

${\mathbf{X}}_t$ are exogenous regressors (e.g., calendar month dummies),

-

${\mathbf{Z}}_t$ captures autoregressive and moving average components,

-

$\beta$, $\gamma$ are the associated parameter vectors,

-

${\epsilon }_t$ is the residual term.

We consider two specifications for ${\epsilon }_t$:

• ${\epsilon }_t\sim \mathrm{SN}({\mu }_t,{\sigma }^2,\lambda )$: Skew-Normal distribution with location ${\mu }_t$, scale $\sigma$, and skewness $\lambda$.
• ${\epsilon }_t\sim \mathrm{ZISN}({\mu }_t,{\sigma }^2,\lambda ,p)$: Zero-Inflated Skew-Normal distribution, combining a point mass at zero with a Skew-Normal component.

The Skew-Normal (SN) density is given by:

$$f_{\mathrm{SN}}\left(\epsilon \right)={\displaystyle \frac{2}{\sigma }}\varphi \left({\displaystyle \frac{\epsilon -\mu }{\sigma }}\right)\mathrm{\Phi }\left(\lambda ·{\displaystyle \frac{\epsilon -\mu }{\sigma }}\right)$$

(2)

where $\varphi (·)$ and $\mathrm{\Phi }(·)$ are the standard normal PDF and CDF, respectively.

The Zero-Inflated version introduces a mixture component:

$$f_{\mathrm{ZISN}}\left(\epsilon \right)=p·{\delta }_0\left(\epsilon \right)+(1-p)·f_{\mathrm{SN}}\left(\epsilon \right)$$

(3)

where ${\delta }_0\left(\epsilon \right)$ is the Dirac delta at zero and $p\in [0,1]$ is the zero-inflation probability.

Parameter estimation proceeds in two stages:

1.

Maximum Likelihood Estimation (MLE) of the SARIMAX baseline parameters ($\beta ,\gamma$),

2.

Expectation-Maximization (EM) algorithm for structured residual parameters $\lambda$ and p in the SN/ZISN models, following the approach described in [5].

Multilayer Perceptron (MLP) Benchmark

As a non-linear forecasting benchmark, we trained a feed-forward Multilayer Perceptron (MLP) regressor. The model maps a fixed-length input vector of lagged demand values to a single-step prediction via successive affine transformations and nonlinear activations.

Let ${\mathbf{u}}_t=[y_{t-1},y_{t-2},y_{t-3}]\in {\mathbb{R}}^3$ be the input vector at time t. The MLP computes:

$${\widehat{y}}_t=f^{\left(3\right)}\left(W^{\left(3\right)}{f}^{\left(2\right)}\left(W^{\left(2\right)}{f}^{\left(1\right)}\left(W^{\left(1\right)}{\mathbf{u}}_t+{\mathbf{b}}^{\left(1\right)}\right)+{\mathbf{b}}^{\left(2\right)}\right)+{\mathbf{b}}^{\left(3\right)}\right)$$

(4)

where:

-

$W^{\left(l\right)}\in {\mathbb{R}}^{d_l\times d_{l-1}}$ and ${\mathbf{b}}^{\left(l\right)}\in {\mathbb{R}}^{d_l}$ are the weight matrices and bias vectors for layer $l=1,2,3$,

-

$f^{\left(1\right)}(·)$, $f^{\left(2\right)}(·)$ are ReLU activation functions $\mathrm{ReLU}\left(x\right)=max(0,x)$,

-

$f^{\left(3\right)}(·)$ is the identity function (linear output),

-

${\widehat{y}}_t$ is the predicted demand at time t.

In our implementation, the hidden layer dimensions are $d_1=32$, $d_2=16$, and $d_3=1$. The model is trained to minimize the Mean Squared Error (MSE) over the training set:

$$\mathcal{L}={\displaystyle \frac{1}{T}}\sum _{t=1}^T{(y_t-{\widehat{y}}_t)}^2$$

(5)

Optimization is performed using the Adam algorithm with learning rate $\eta =0.001$, for up to 1,000 epochs. Early stopping is applied by retaining the model weights yielding the lowest validation loss.

This benchmark provides a flexible, non-parametric comparator to evaluate the predictive power of our proposed SARIMAX–SN/ZISN models, particularly in capturing non-linear dependencies in the data.

2.3. Optimization Phase (Global PSO)

The forecasted monthly demand ${\widehat{D}}_i$ for each laboratory reagent was used to compute the optimal order quantities $\mathbf{Q}=(Q_1,\dots ,Q_N)$ by solving a multivariate, non-linear and constrained inventory control problem.

We define:

$$Q_i=k_i·{\mathrm{pack}}_i$$

(6)

where $k_i\in {\mathbb{N}}_0$ is the number of packaging units to be ordered for item i.

The total procurement cost is defined as:

$$C\left(\mathbf{k}\right)=\sum _{i=1}^N\left[c_i{Q}_i+{\mathbb{I}}_{\{Q_i>0\}}{o}_i+h_i·{(Q_i-{\widehat{D}}_i)}^{+}+s_i·{({\widehat{D}}_i-Q_i)}^{+}\right]$$

(7)

subject to the global budget constraint:

$$\sum _{i=1}^N{c}_i{Q}_i\le B$$

(8)

where:

-

$c_i$ is the unit cost,

-

$o_i$ is the fixed ordering cost,

-

$h_i$ is the holding cost per excess unit,

-

$s_i$ is the shortage cost per missing unit,

-

B is the total available budget.

Because the problem is combinatorial, nonlinear, and non-differentiable, we apply Particle Swarm Optimization (PSO) to find a near-optimal solution ${\mathbf{k}}^{*}$. Each particle represents an integer vector $\mathbf{k}=(k_1,\dots ,k_N)$, initialized randomly in $[0,k_i^{max}]$ for each item.

To enforce the budget constraint, we use a penalized cost function:

$$\tilde{C}\left(\mathbf{k}\right)=C\left(\mathbf{k}\right)+\rho ·max\left(0,\sum _{i=1}^N{c}_i{Q}_i-B\right)$$

(9)

with $\rho \gg 1$ being a penalty coefficient to discourage budget violations.

2.3.1. Mathematical Formulation of PSO

Particle Swarm Optimization is a population-based metaheuristic inspired by the collective behavior of swarms [8]. Each particle has a position ${\mathbf{k}}^{\left(t\right)}$ and a velocity ${\mathbf{v}}^{\left(t\right)}$, both updated iteratively as follows:

$${\mathbf{v}}^{(t+1)}=w·{\mathbf{v}}^{\left(t\right)}+c_1·r_1·({\mathbf{p}}^{\left(t\right)}-{\mathbf{k}}^{\left(t\right)})+c_2·r_2·({\mathbf{g}}^{\left(t\right)}-{\mathbf{k}}^{\left(t\right)})$$

(10)

$${\mathbf{k}}^{(t+1)}={\mathbf{k}}^{\left(t\right)}+{\mathbf{v}}^{(t+1)}$$

(11)

where:

-

w is the inertia weight (balances exploration and exploitation),

-

$c_1$, $c_2$ are cognitive and social acceleration coefficients,

-

$r_1$, $r_2\sim \mathcal{U}(0,1)$ are random numbers,

-

${\mathbf{p}}^{\left(t\right)}$ is the personal best position of the particle,

-

${\mathbf{g}}^{\left(t\right)}$ is the global best found by the swarm.

As our problem requires integer decisions with packaging constraints, we apply rounding after each position update and enforce Eq. (6) to compute feasible order quantities.

The PSO configuration used 50 particles, 200 iterations, inertia weight $w=0.7$, and acceleration coefficients $c_1=c_2=1.5$. The best feasible solution ${\mathbf{k}}^{*}$ yields final order quantities ${\mathbf{Q}}^{*}={\mathbf{k}}^{*}\circ \mathrm{pack}$.

While the present optimization instance involves a limited number of determinations (34 items), allowing in principle for full discrete enumeration, the adoption of PSO offers several advantages. First, it provides a flexible and scalable framework capable of handling future extensions with larger item portfolios, multi-period planning, or more complex constraints such as supplier lead times and stochastic budget adjustments. Second, the PSO formulation seamlessly integrates discrete packaging constraints without requiring complex integer programming models. This design ensures that the optimization layer remains generalizable and readily applicable to broader healthcare supply chain contexts beyond the scope of the present study.

2.4. Code Availability

2.5. Reproducible Workflow for Model Estimation

The Python snippet below reproduces the full pipeline used in this study, incorporating residual centering, outlier truncation, and lambda regularization for stable estimation of Skew-Normal parameters. The workflow covers data loading, model order selection, residual distribution fitting, forecasting, and non-linear benchmark training using only pandas, statsmodels, and scikit-learn.

Listing 1: Full forecasting workflow: SARIMAX order selection, residual centering, Skew-Normal fitting, forecasting, and MLP benchmark.

The script: (1) reads the raw Excel file; (2) performs AIC-driven SARIMAX order selection; (3) centers and truncates residuals before fitting Skew-Normal (or Zero-Inflated Skew-Normal) distributions with lambda regularization; (4) forecasts the final six months; and (5) trains a simple MLP benchmark with global normalization for comparison.

2.6. Inventory–Cost Optimisation Workflow

The algorithm below links the stochastic–forecasting layer with the inventory–control layer. For each laboratory reagent, it

1.

loads the 36-month demand series from DemandaSpecialIsue.xlsx;

2.

fits a fast SARIMAX$(0,1,2)$ model and extracts the one-month ahead mean demand $\widehat{D}$;

3.

enumerates integer multiples $k_i\in \{0,\dots ,20\}$ of the pack size $p_i$ and evaluates the total cost

$$C\left(Q_i\right)=c_i{Q}_i+{\mathbf{1}}_{\{Q_i>0\}}{c}_o+h_i{[Q_i-\widehat{D}]}^{+}+s_i{[\widehat{D}-Q_i]}^{+},Q_i=k_i{p}_i;$$

4.

stores the optimal lot $Q_i^{*}$, the current lot $Q_i^{\mathrm{curr}}$ (December-2023 order) and the corresponding costs.

A pyswarm (PSO) call can replace the brute-force loop if a larger search space is required; here the discrete sweep is sufficient and does not depend on external libraries.

Listing 2: Global PSO optimisation with packaging constraints and budget limit.

The resulting pso_inventory_results.csv to report the optimal order multiple, cost savings, and compliance with the packaging constraint $Q^{*}=k_i{p}_i$ for all 34 determinations.

3. Results

This section presents the performance of the forecasting models, residual diagnostics, and inventory cost optimisation using the calibrated demand forecasts.

3.1. Forecast Accuracy and Parameter Estimates

Table 1 summarizes the fitted distribution parameters and forecast accuracy metrics (MAE and RMSE).

Regarding the fitted location parameter $\mu$, several estimates were found to be substantially negative despite prior centering of residuals. This arises because residual centering was applied globally before truncation, but the skew-normal fitting procedure operates on truncated subsets of residuals. Therefore, $\mu$ reflects the optimal location shift of the skew-normal kernel fitted to the cleaned residuals, capturing potential shifts induced by asymmetric trimming. Negative $\mu$ values may indicate a tendency of the SARIMAX forecasts to slightly overestimate actual demand on average.

Out of the 34 analytes, only Phenobarbital exhibited a zero-proportion above the 10 % threshold and was therefore assigned the SARIMAX–ZISN specification. The remaining 33 laboratory reagents were modelled with SARIMAX–SN. Across the cohort, SARIMAX–SN reduced mean MAE by ~11 % with respect to a Gaussian-error baseline (not shown). The single ZISN case (Phenobarbital ) improved MAE by 14 % relative to its SN counterpart, confirming the benefit of explicitly modelling excess zeros.

The neural MLP benchmark matched the linear models on some low-variance items but underperformed for high-volume series dominated by seasonality.

For the MLP benchmark (three lagged inputs) we computed permutation feature importance on the hold-out window; lags $t-1$ and $t-2$ systematically contributed $>80\%$ of the explained variance, indicating that the network mostly captures short-term autocorrelation rather than complex seasonal patterns.

Table 2 summarises the forecasting accuracy across three approaches. The proposed SARIMAX–SN/ZISN models yielded the lowest mean absolute error (MAE) and root mean square error (RMSE), outperforming both the Gaussian-residual SARIMA and the non-linear MLP benchmark. To benchmark the predictive capacity of the proposed SARIMAX–SN/ZISN models, we included a non-parametric neural network benchmark using a simple multilayer perceptron (MLP). While the MLP offers flexible nonlinear approximation capabilities, its black-box nature lacks interpretability, particularly regarding residual asymmetry and explainability required in clinical forecasting contexts. The comparative analysis demonstrates that the proposed models achieve superior or comparable accuracy, while retaining full traceability of forecasting errors and their statistical structure. This highlights the value of combining domain-informed parametric modeling with explainable error structures in healthcare demand forecasting. The estimated skewness parameter $\lambda$ was positive for $32/34$ laboratory reagents, confirming the right–tailed behaviour highlighted by skew–normal diagnostics, whereas the only analyte with a substantial zero proportion (Phenobarbital, $p_0\approx 0.30$) required the zero–inflated specification. These diagnostics justify the asymmetric error structures adopted.

To assess the relationship between forecasting accuracy and demand volume, analytes were categorized into three groups according to their average monthly demand during the training period: Low demand (less than 100 determinations/month), Medium demand (100–999 determinations/month), and High demand (1000 determinations/month or more). Table 3 presents the disaggregated MAE and RMSE values for each demand level.

As expected, absolute errors increase with demand volume. However, this reflects natural scale effects rather than model inefficiency. The forecast errors remain reasonably proportional to demand magnitude, as further explored through percentage-based metrics in subsequent sections.

Table 4 reports forecasting errors for each forecast step from one to six months ahead. Interestingly, the second forecast step exhibits the highest MAE and RMSE values, potentially reflecting short-term fluctuations or abrupt changes not fully captured by the autoregressive structure. Beyond this initial peak, errors stabilize across longer forecast horizons, suggesting that the models maintain robust performance over multiple steps ahead.

In addition to absolute error metrics, the Mean Absolute Percentage Error (MAPE) was computed to evaluate forecasting performance relative to demand volume. Table 5 summarizes the MAPE by demand level.

MAPE could not be calculated for the low-demand group due to zero-demand occurrences in the test set, which make percentage-based measures undefined. For the medium- and high-demand categories, the models achieved average MAPE values below 11%, indicating strong relative accuracy even for high-volume determinations.

3.2. Residual Diagnostics and Model Selection

After fitting each SARIMAX–SN / SARIMAX–ZISN model we inspected the standardised residuals via Ljung-Box tests and QQ-plots, see its results in Table 6 and Figure 1, respectively.

To evaluate residual distributional assumptions, QQ-plots were constructed by comparing standardized residuals against the standard normal distribution. While the Skew-Normal specification introduces inherent asymmetry, QQ-plots remain informative to visually assess overall goodness-of-fit and detect residual heavy tails or misspecification. Future work could improve residual diagnostics by generating QQ-plots directly against the fitted Skew-Normal quantiles, providing a more precise graphical validation of the estimated shape parameters.

3.3. Optimization Outcomes

Table 7 summarises the inputs used in the PSO-based optimisation stage. For each determination, we report the best-fitting SARIMAX model order, the expected monthly demand (estimated as the mean forecast over a six-month horizon), and the associated cost parameters. These include the unit cost of purchase, the fixed cost of placing an order, and penalty terms for holding excess stock or facing shortages. Packaging constraints are also specified via the ‘Pack_Size‘ column, which enforces that each order quantity must be an integer multiple of a fixed unit.

The ordering recommendations obtained with the PSO optimisation layer were then fed into the inventory cost model; the resulting savings relative to the empirical hospital policy are summarised in Table 8. This end-to-end traceability—from residual shape parameters, predictive errors, and zero-inflation detection to constrained ordering policies—provides an explainable, domain-specific hybrid framework, fulfilling the design goals stated in the introduction. The optimiser strictly satisfies the packaging constraint $Q_i^{*}=k_i·{\mathrm{pack}}_i$ for each item, as detailed in Table 8. The fill rate remained above 85% for all high-volume and above 70% across the full portfolio.

Table 9 summarises monthly inventory costs across three procurement policies. The hospital’s empirical policy resulted in an average monthly cost of CLP 179.5 million. A baseline SARIMA-based policy, using Gaussian residuals and discrete enumeration, reduced this to CLP 32.1 million. Our proposed hybrid SARIMAX–SN/ZISN forecasting integrated with PSO-constrained optimisation further lowered the monthly cost to CLP 19.6 million, yielding an overall savings of 89%.

Future work will consider multi-objective optimisation approaches to simultaneously minimise cost and maximise service levels.

4. Discussion

The hybrid framework successfully integrates statistical forecasting and metaheuristic optimisation to address the dual challenges of demand uncertainty and procurement constraints in clinical laboratories.

4.1. Interpretation of Main Results and Operational Implications

Improved Forecast Accuracy with SARIMAX–SN/ZISN Models

The SARIMAX models with Skew-Normal and Zero-Inflated residuals consistently outperformed both traditional Gaussian-based SARIMA models and neural benchmarks (MLP) in terms of MAE and RMSE. The presence of positive skewness in 32 out of 34 determinations confirmed the right-tailed nature of demand distributions, while the detection of excess zeros (notably in Phenobarbital ) justified the use of ZISN.

Operational implication: These improvements enable procurement teams to anticipate demand more precisely, especially for laboratory reagents with irregular or sporadic use. By aligning forecasts with actual demand patterns, the risk of stockouts or surplus inventory is significantly reduced—critical for time-sensitive clinical workflows.

Robustness and Interpretability of Residual Modeling.

The residual diagnostics (Ljung–Box tests and QQ plots) confirmed the adequacy of model fits, validating the assumption of structured asymmetry. The residual skewness parameter $\lambda$ provides not only a technical advantage but also a practical one—revealing when demand is likely to spike unexpectedly.

Operational implication: This insight allows for the proactive allocation of buffer stock for high-skew laboratory reagents and informs clinicians or lab managers about critical variability in usage, thus improving readiness for high-demand periods (e.g., seasonal surges, pandemics).

Inventory Cost Optimization under Realistic Constraints

The Particle Swarm Optimization (PSO) layer identified procurement plans that minimized total cost while satisfying both budget and packaging constraints. Compared to the hospital’s empirical policy, the hybrid framework reduced average monthly inventory costs from CLP 179.5 million to CLP 19.6 million—a savings of nearly 89%.

Operational implication: These cost reductions translate into real resource availability for other clinical operations. Additionally, enforcing packaging multiples (e.g., boxes of 12 or 24) during optimization ensures that procurement recommendations are logistically implementable, avoiding fragmented or infeasible orders.

High Service Levels Maintained Across Portfolio

Despite cost reductions, service levels remained above 85% for high-volume determinations and over 70% across all items, indicating no compromise on availability.

Operational implication: This confirms the viability of the proposed strategy in balancing cost efficiency with clinical reliability. It supports strategic inventory planning that avoids false trade-offs between cost and care quality.

4.2. Implications for Clinical Laboratory Logistics and Supply

From a logistics management standpoint, the proposed hybrid framework provides meaningful enhancements to procurement decision-making in clinical laboratories. By explicitly linking demand forecasts to ordering processes, the model addresses a well-documented disconnect that often results in emergency purchases and high wastage rates—an issue highlighted in public hospitals in Ethiopia, where poor inventory alignment led to an estimated 27% laboratory reagent waste [1]. The framework also incorporates key operational constraints such as procurement cycles, supplier limitations, and cost structures, ensuring that output recommendations are not only optimized but also feasible within real-world institutional contexts. This is particularly relevant in resource-constrained settings, where studies have emphasized the value of flexible systems that can be modularly implemented, enabling laboratories to adopt forecasting or optimization components independently depending on capacity [12]. Moreover, the model’s transparency enhances traceability and auditability—crucial in public health institutions where procurement decisions are subject to strict oversight and accountability [1]. Empirical evidence further supports the model’s potential to reduce operational costs and improve supply continuity; integrated demand-planning strategies have been shown to decrease inventory levels by up to 40% and reduce order frequency by 60% in clinical commodities such as blood components, while also reducing stock-out rates without compromising service availability [13]. Taken together, this hybrid approach not only improves forecasting accuracy but also strengthens resource sustainability, cost-efficiency, and supply chain resilience in public clinical laboratory systems.

4.3. Limitations of the Study

The main limitations are:

• Scope restriction: The study used data from a single mid-sized Chilean hospital, potentially limiting generalization to larger institutions or different healthcare systems.
• Static parameters: Procurement costs, shortage penalties, and budget levels were assumed fixed. In practice, these may vary monthly.
• Temporal assumptions: The model does not incorporate long-term trend shifts (e.g., due to epidemiological or technological changes).
• Simplified logistics: Lead times and supplier delays were not included in the cost function, though they may impact stock availability.

4.4. Future Research Directions

Building on the current framework, future work could explore:

• Multi-objective optimisation to jointly minimise cost and maximise service levels.
• Multi-period planning under rolling budgets and lead times.
• Integration with hospital ERP systems for real-time updates and automatic re-planning.
• Extension to regional networks with shared laboratory reagent pools or cooperative procurement models.
• Decision dashboards for non-technical stakeholders, integrating explainability with usability.

5. Conclusions

This study proposed and validated a hybrid forecasting–optimization framework designed to improve inventory management in clinical laboratories operating under demand uncertainty, budget restrictions, and packaging constraints. By combining SARIMAX models with structured residuals (Skew-Normal and Zero-Inflated Skew-Normal) and a metaheuristic optimisation layer based on Particle Swarm Optimization (PSO), the framework achieved significant advances in both predictive accuracy and cost efficiency.

The main contributions and conclusions are as follows:

• The proposed SARIMAX–SN/ZISN models consistently outperformed standard SARIMA and neural network benchmarks in forecasting accuracy, particularly for laboratory reagents exhibiting skewed or zero-inflated demand.
• The metaheuristic inventory optimisation component effectively translated improved forecasts into procurement decisions that were budget-compliant, packaging-feasible, and highly cost-efficient, achieving up to 89% monthly cost savings compared to the hospital’s empirical policy.
• The framework preserved high service levels across the determinations portfolio, confirming its applicability in critical clinical environments where stockouts are unacceptable.
• The integration of explainable forecasting structures and constrained optimisation enhances transparency and traceability from data to decisions—essential for implementation in public healthcare systems.
• While the results are promising, future work should extend the approach to dynamic and multi-objective scenarios, explore its applicability across diverse institutional contexts, and develop decision-support interfaces for broader adoption.

In summary, the hybrid framework presented here offers a robust, interpretable, and operationally grounded solution for laboratory reagent inventory planning in clinical laboratories, balancing cost-efficiency with clinical reliability.