This is somewhat an abstract question. In Mathematica there is a function TimeSeriesModelFit which automatically finds the best autoregressive model fitting the data. Is there a way to do the same in Python? Might it be there is a package for this?
This is the data from Mathematica itself as an example:
data = {5., 9., 8., 10., 6.1, 10.4, 9.1, 11.6, 7.5, 12.1, 10.4, 13.5,
9., 14.1, 11.9, 15.7, 10.8, 16.4, 13.7, 18.3, 12.9, 19., 15.8, 21.2,
15.3, 22.1, 18.3, 24.6};
statsmodels.tsa.ar_model.ar_select_order can do that:
>>> from statsmodels.tsa import ar_model
>>> search = ar_model.ar_select_order(endog=data, maxlag=10, ic="aic", glob=True)
>>> search.ar_lags
array([ 1, 4, 9, 10])
It looks up until lag 10; but since glob is enabled, it considers every 2**10 = 1024 possible lag combinations. Note that this means it doesn't try only 10 models (i.e., AR(1) up to AR(10)), but it turns on and off lags (which is also seen from the above output's discontinous lag suggestion).
statsmodels lets you pass these not necessarily continuous lag terms for forming a new model when refitting, e.g.,
>>> mod = ar_model.AutoReg(endog=data, lags=[1, 4, 9, 10])
>>> res = mod.fit()
>>> res.aic
-6.380532624300298
For the models that are not necessarily purely autoregressive, there is pmdarima package that implements automatic selection of S-ARIMA-X models. It is based on auto.arima of R. I'm not familiar with ARCH-like models but it seems the arch package's implementation can be wrapped in a loop to try out different models and select with some criterion.