Time Series prediction with multiple features in the input data

SARIMAX, as a generalisation of the SARIMA model, is designed to handle exactly this. From the docs,

Parameters:

• endog (array_like) – The observed time-series process y;
• exog (array_like, optional) – Array of exogenous regressors, shaped (nobs, k).

You could pass the holiday and promotion_day as an array of size (nobs, 2) to exog, which will inform the model of the exogenous nature of some of these observations.