SARIMAX Statsmodels fits a regression with ARMA errors. So the model is:
yt=βxt+ut
ut=μ+ηt−1+ζt
Since p and q orders only affect the second equation, why do I get different results (coefs) of the exogenous variables? What am i missing?
Note that while there are two separate equations, the fitting is not done in two separate steps. (Your question sounds like you are assuming this, which is quite understandable, but statsmodels.sarimax
does not work this way.) Rather, there is a single estimation step in which all the parameters are estimated - those for the regression and those for the SARIMA process on the residuals. And of course, if the SARIMA model order changes, then all parameters can (and usually will) change: those of the SARIMA error process and those of the regression on predictors.
Here is an earlier thread that covers regression with ARIMA errors (not SARIMA, but that doesn't make a difference here), in the R forecast
package (which the various ARIMA functionalities in Python are patterned on, because it's the gold standard), by the author of the forecast
package for R himself: ARIMAX vs. Regression With ARIMA Errors.