How to extract seasonality of a time series when controlling for an additional regressor? My goal is to extract the seasonality of a product's grocery sales time series.
I know standard time series decomposition when a time series is expressed as $_=_+_+_$, the sum (for an additive model) of a trend component, a seasonal component and the remainder. Currently in Python, I use the seasonal_decompose function from the statsmodels library.
I know that wit a certain regularity, the seller gives a discount on the product in form of a promotion. Since I know the exact periods of the promotion, I would like to control for this information before decomposing the time series in order to make sure that the promotional effect is not falsely attributed to the seasonality. What methods/libraries could I use to control for this additional regressor?

A simple approach would be to first fit a regression of sales on explanatory variables capturing the promotion. This could be anything from a single Boolean dummy regressor to a larger model capturing price points, discount percentages, or any other attributes of the promotions. Promotions can get very complicated indeed. You can start with a simple model and analyze it for as long as your time box allows, just be sure not to overparameterize it too much.
Then take the residuals from that regression and decompose those, using STL or whatever other method you prefer, like Exponential Smoothing.
This approach is very similar to a "regression with ARIMA models", which is what ARIMA models typically fit in the presence of explanatory variables. Rob Hyndman's blog post on "The ARIMAX model muddle" gives a very nice overview about the differences between this and an ARIMAX model, much of which is applicable to your decomposition question.
