Can the ARIMA function in statsmodel include covariates? I have already implemented a univariate time series forecasting using ARIMA in statsmodel. I know that in R we can put external reggressors using auto ARIMA. Can this be done with statsmodel's ARIMA implementation? 
ex: a store has different features like sales, dayOfWeek, promotion and schoolHoliday. I want to use all these features of day d-1 to predict sales of d. Can this be done using the ARIMA implementation of the stats model?
 A: I think you are looking for the 'exog' parameter. It allows you to include the predictors you mention into the model. Your line of code could look something like this:
exog = train.loc['dayofweek','promotion',...]
model = sm.tsa.statespace.SARIMAX(train.sales, order = (p,d,q),
                                  seasonal_order = (P,D,Q,m), 
                                  trend = 'n', exog = exog)
Ps: If you want your to include the lagged value of your exog variable, don't forget to shift your values by 1 day. This can be done using the dataframe.shift() method.
A: It is prudent to test for contemporaneous and lag effects for some user specified variables. Often there can be a lead effect where sales are affected the day before a price change at the next period. Care should be taken to detect and incorporate latent deterministic structure ala http://docplayer.net/12080848-Outliers-level-shifts-and-variance-changes-in-time-series.html as they have the potential to unmask effects that might otherwise be discarded as non-significant.
Identifying the form of the response between a dependent series and stochastic inputs is best handled by following http://www.math.cts.nthu.edu.tw/download.php?filename=569_fe0ff1a2.pdf&dir=publish&title=Ruey+S.+Tsay-Lec1 and https://web.archive.org/web/20160216193539/https://onlinecourses.science.psu.edu/stat510/node/75/
I don't believe these suggested identification approaches are part of the fitting software that you are considering and you might be better served by alternatives.
