This study proposes an extension of Seasonal Autoregressive Integrated Moving Average models with exogenous regressors (SARIMAX) by incorporating skew-normal and zero-inflated skew-normal error structures to better accommodate asymmetry and excess zeros in time series data. The proposed framework demonstrates improved flexibility and robustness compared to traditional Gaussian-based models. Simulation experiments reveal that the skewness parameter significantly affect forecasting accuracy, with reductions in mean absolute error (MAE) and root mean square error (RMSE) observed across both positively and negatively skewed scenarios. Notably, in negative-skew contexts, the model achieved an MAE of 0.40 and RMSE of 0.49, outperforming its symmetric-error counterparts. The inclusion of zero-inflation probabilities further enhances model performance in sparse datasets, yielding superior values in goodness-of-fit criteria such as the Akaike Information Criterion (AIC) and Bayesian Information Criterion (BIC). To illustrate the practical value of the methodology, a real-world case study is presented involving the modeling of optical density (OD600) data from Escherichia coli during stationary-phase growth. A SARIMAX(1,1,1) model with skew-normal errors was fitted to 200 time-stamped absorbance measurements, revealing significant positive skewness in the residuals. Bootstrap-derived confidence intervals confirmed the significance of the estimated skewness parameter (α=14.033 with 95% CI [12.07, 15.99]). The model outperformed the classical ARIMA benchmark in capturing the asymmetry of the stochastic structure, underscoring its relevance for biological, environmental, and industrial applications in which non-Gaussian features are prevalent.