等价于 python 的 auto.arima()
Posted
技术标签:
【中文标题】等价于 python 的 auto.arima()【英文标题】：auto.arima() equivalent for python 【发布时间】：2014-05-11 07:26:43 【问题描述】：我正在尝试使用 ARMA ARIMA 模型预测每周销售额。我在statsmodels 中找不到用于调整 order(p,d,q) 的函数。目前 R 有一个函数forecast::auto.arima()，它将调整 (p,d,q) 参数。
如何为我的模型选择正确的订单？ python中是否有为此目的可用的库？
【问题讨论】：
【参考方案1】：您可以实施多种方法：
ARIMAResults 包括aic 和bic。根据他们的定义（参见here 和here），这些标准会影响模型中的参数数量。因此，您可以使用这些数字来比较模型。 scipy 还有optimize.brute，它在指定的参数空间上进行网格搜索。所以像这样的工作流程应该可以工作：
def objfunc(order, exog, endog):
from statsmodels.tsa.arima_model import ARIMA
fit = ARIMA(endog, order, exog).fit()
return fit.aic()
from scipy.optimize import brute
grid = (slice(1, 3, 1), slice(1, 3, 1), slice(1, 3, 1))
brute(objfunc, grid, args=(exog, endog), finish=None)
确保您使用finish=None 致电brute。
您可以从ARIMAResults 获得pvalues。因此，一种步进算法很容易实现，其中模型的度数在维度上增加，从而获得添加参数的最低 p 值。
使用ARIMAResults.predict 交叉验证替代模型。最好的方法是将时间序列的尾部（例如最近 5% 的数据）保留在样本之外，并使用这些点来获得拟合模型的测试误差。
【讨论】：
master 中已经有一个包装函数可以为您执行此操作。 statsmodels.sourceforge.net/devel/generated/… 这是我们目前最接近 auto.arima 的地方。 @jseabold 你知道源代码，但名字暗示它的 arma 不是 arima 是的，仅用于订单选择。目前还没有任何自动检查集成。 没关系，对于任何对此也感到困惑的人，exog/endog 是正在拟合的数据。 Brute 自动使用它的第二个参数作为函数的第一个参数，然后按照列出的顺序使用其他 args。 ARIMA（带集成）的订单选择似乎是available now。【参考方案2】：其实
def objfunc(order,*params ):
from statsmodels.tsa.arima_model import ARIMA
p,d,q = order
fit = ARIMA(endog, order, exog).fit()
return fit.aic()
from scipy.optimize import brute
grid = (slice(1, 3, 1), slice(1, 3, 1), slice(1, 3, 1))
brute(objfunc, grid, args=params, finish=None)
【讨论】：
虽然这段代码 sn-p 可以解决问题，但including an explanation 确实有助于提高帖子的质量。请记住，您正在为将来的读者回答问题，而这些人可能不知道您的代码建议的原因。也请尽量不要用解释性的 cmets 挤满你的代码，这会降低代码和解释的可读性！【参考方案3】：我编写了这些实用函数来直接计算 pdq 值 get_PDQ_parallel 需要三个以时间戳（日期时间）为索引的系列输入数据。 n_jobs 将提供并行处理器的数量。输出将是具有 aic 和 bic 值的数据帧，索引中的 order=(P,D,Q) p 和 q 范围是 [0,12] 而 d 是 [0,1]
import statsmodels
from statsmodels import api as sm
from sklearn.metrics import r2_score,mean_squared_error
from sklearn.utils import check_array
from functools import partial
from multiprocessing import Pool
def get_aic_bic(order,series):
aic=np.nan
bic=np.nan
#print(series.shape,order)
try:
arima_mod=statsmodels.tsa.arima_model.ARIMA(series,order=order,freq='H').fit(transparams=True,method='css')
aic=arima_mod.aic
bic=arima_mod.bic
print(order,aic,bic)
except:
pass
return aic,bic
def get_PDQ_parallel(data,n_jobs=7):
p_val=13
q_val=13
d_vals=2
pdq_vals=[ (p,d,q) for p in range(p_val) for d in range(d_vals) for q in range(q_val)]
get_aic_bic_partial=partial(get_aic_bic,series=data)
p = Pool(n_jobs)
res=p.map(get_aic_bic_partial, pdq_vals)
p.close()
return pd.DataFrame(res,index=pdq_vals,columns=['aic','bic'])
【讨论】：
【参考方案4】：可能的解决方案
df=pd.read_csv("http://vincentarelbundock.github.io/Rdatasets/csv/datasets/AirPassengers.csv")
# Define the p, d and q parameters to take any value between 0 and 2
p = d = q = range(0, 2)
print(p)
import itertools
import warnings
# Generate all different combinations of p, q and q triplets
pdq = list(itertools.product(p, d, q))
print(pdq)
# Generate all different combinations of seasonal p, q and q triplets
seasonal_pdq = [(x[0], x[1], x[2], 12) for x in list(itertools.product(p, d, q))]
print('Examples of parameter combinations for Seasonal ARIMA...')
print('SARIMAX: x '.format(pdq[1], seasonal_pdq[1]))
print('SARIMAX: x '.format(pdq[1], seasonal_pdq[2]))
print('SARIMAX: x '.format(pdq[2], seasonal_pdq[3]))
print('SARIMAX: x '.format(pdq[2], seasonal_pdq[4]))
Examples of parameter combinations for Seasonal ARIMA...
SARIMAX: (0, 0, 1) x (0, 0, 1, 12)
SARIMAX: (0, 0, 1) x (0, 1, 0, 12)
SARIMAX: (0, 1, 0) x (0, 1, 1, 12)
SARIMAX: (0, 1, 0) x (1, 0, 0, 12)
y=df
#warnings.filterwarnings("ignore") # specify to ignore warning messages
for param in pdq:
for param_seasonal in seasonal_pdq:
try:
mod = sm.tsa.statespace.SARIMAX(y,
order=param,
seasonal_order=param_seasonal,
enforce_stationarity=False,
enforce_invertibility=False)
results = mod.fit()
print('ARIMAx12 - AIC:'.format(param, param_seasonal, results.aic))
except:
continue
ARIMA(0, 0, 0)x(0, 0, 1, 12)12 - AIC:3618.0303991426763
ARIMA(0, 0, 0)x(0, 1, 1, 12)12 - AIC:2824.7439963684233
ARIMA(0, 0, 0)x(1, 0, 0, 12)12 - AIC:2942.2733127230185
ARIMA(0, 0, 0)x(1, 0, 1, 12)12 - AIC:2922.178151133141
ARIMA(0, 0, 0)x(1, 1, 0, 12)12 - AIC:2767.105066400224
ARIMA(0, 0, 0)x(1, 1, 1, 12)12 - AIC:2691.233398643673
ARIMA(0, 0, 1)x(0, 0, 0, 12)12 - AIC:3890.816777796087
ARIMA(0, 0, 1)x(0, 0, 1, 12)12 - AIC:3541.1171286722
ARIMA(0, 0, 1)x(0, 1, 0, 12)12 - AIC:3028.8377323188824
ARIMA(0, 0, 1)x(0, 1, 1, 12)12 - AIC:2746.77973129136
ARIMA(0, 0, 1)x(1, 0, 0, 12)12 - AIC:3583.523640623017
ARIMA(0, 0, 1)x(1, 0, 1, 12)12 - AIC:3531.2937768990187
ARIMA(0, 0, 1)x(1, 1, 0, 12)12 - AIC:2781.198675746594
ARIMA(0, 0, 1)x(1, 1, 1, 12)12 - AIC:2720.7023088205974
ARIMA(0, 1, 0)x(0, 0, 1, 12)12 - AIC:3029.089945668332
ARIMA(0, 1, 0)x(0, 1, 1, 12)12 - AIC:2568.2832251221016
ARIMA(0, 1, 0)x(1, 0, 0, 12)12 - AIC:2841.315781459511
ARIMA(0, 1, 0)x(1, 0, 1, 12)12 - AIC:2815.4011044132576
ARIMA(0, 1, 0)x(1, 1, 0, 12)12 - AIC:2588.533386513587
ARIMA(0, 1, 0)x(1, 1, 1, 12)12 - AIC:2569.9453272483315
ARIMA(0, 1, 1)x(0, 0, 0, 12)12 - AIC:3327.5177587522303
ARIMA(0, 1, 1)x(0, 0, 1, 12)12 - AIC:2984.716706112334
ARIMA(0, 1, 1)x(0, 1, 0, 12)12 - AIC:2789.128542154043
ARIMA(0, 1, 1)x(0, 1, 1, 12)12 - AIC:2537.0293659293943
ARIMA(0, 1, 1)x(1, 0, 0, 12)12 - AIC:2984.4555708516436
ARIMA(0, 1, 1)x(1, 0, 1, 12)12 - AIC:2939.460958374472
ARIMA(0, 1, 1)x(1, 1, 0, 12)12 - AIC:2578.7862352774437
ARIMA(0, 1, 1)x(1, 1, 1, 12)12 - AIC:2537.771484229265
ARIMA(1, 0, 0)x(0, 0, 0, 12)12 - AIC:3391.5248913820797
ARIMA(1, 0, 0)x(0, 0, 1, 12)12 - AIC:3038.142074281268
C:\Users\Dell\Anaconda3\lib\site-packages\statsmodels\base\model.py:496: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals
"Check mle_retvals", ConvergenceWarning)
ARIMA(1, 0, 0)x(0, 1, 0, 12)12 - AIC:2839.809192263449
ARIMA(1, 0, 0)x(0, 1, 1, 12)12 - AIC:2588.50367175184
ARIMA(1, 0, 0)x(1, 0, 0, 12)12 - AIC:2993.4630440139595
ARIMA(1, 0, 0)x(1, 0, 1, 12)12 - AIC:2995.049216326931
ARIMA(1, 0, 0)x(1, 1, 0, 12)12 - AIC:2588.2463284315304
ARIMA(1, 0, 0)x(1, 1, 1, 12)12 - AIC:2592.80110502723
ARIMA(1, 0, 1)x(0, 0, 0, 12)12 - AIC:3352.0350133621478
ARIMA(1, 0, 1)x(0, 0, 1, 12)12 - AIC:3006.5493366627807
ARIMA(1, 0, 1)x(0, 1, 0, 12)12 - AIC:2810.6423724894516
ARIMA(1, 0, 1)x(0, 1, 1, 12)12 - AIC:2559.584031948852
ARIMA(1, 0, 1)x(1, 0, 0, 12)12 - AIC:2981.2250436794675
ARIMA(1, 0, 1)x(1, 0, 1, 12)12 - AIC:2959.3142304724834
ARIMA(1, 0, 1)x(1, 1, 0, 12)12 - AIC:2579.8245645892207
ARIMA(1, 0, 1)x(1, 1, 1, 12)12 - AIC:2563.13922589258
ARIMA(1, 1, 0)x(0, 0, 0, 12)12 - AIC:3354.7462930846423
ARIMA(1, 1, 0)x(0, 0, 1, 12)12 - AIC:3006.702997636003
ARIMA(1, 1, 0)x(0, 1, 0, 12)12 - AIC:2809.3844175191666
ARIMA(1, 1, 0)x(0, 1, 1, 12)12 - AIC:2558.484602766447
ARIMA(1, 1, 0)x(1, 0, 0, 12)12 - AIC:2959.885810636943
ARIMA(1, 1, 0)x(1, 0, 1, 12)12 - AIC:2960.712709764296
ARIMA(1, 1, 0)x(1, 1, 0, 12)12 - AIC:2557.945907092698
ARIMA(1, 1, 0)x(1, 1, 1, 12)12 - AIC:2559.274166458508
ARIMA(1, 1, 1)x(0, 0, 0, 12)12 - AIC:3326.3285511700374
ARIMA(1, 1, 1)x(0, 0, 1, 12)12 - AIC:2985.868532151721
ARIMA(1, 1, 1)x(0, 1, 0, 12)12 - AIC:2790.7677149967103
ARIMA(1, 1, 1)x(0, 1, 1, 12)12 - AIC:2538.820635541546
ARIMA(1, 1, 1)x(1, 0, 0, 12)12 - AIC:2963.2789505804294
ARIMA(1, 1, 1)x(1, 0, 1, 12)12 - AIC:2941.2436984747465
ARIMA(1, 1, 1)x(1, 1, 0, 12)12 - AIC:2559.8258191422606
ARIMA(1, 1, 1)x(1, 1, 1, 12)12 - AIC:2539.712354465328
来自https://www.digitalocean.com/community/tutorials/a-guide-to-time-series-forecasting-with-arima-in-python-3
另见https://github.com/decisionstats/pythonfordatascience/blob/master/time%2Bseries%20(1).ipynb
【讨论】：
【参考方案5】：现在有一个合适的 python 包来做 auto-arima。 https://github.com/tgsmith61591/pmdarima
文档： http://alkaline-ml.com/pmdarima
用法示例：https://github.com/tgsmith61591/pmdarima/blob/master/examples/quick_start_example.ipynb
【讨论】：
图书馆文档可以在这里查阅：pyramid-arima.readthedocs.io/en/latest/index.html 仅供参考，python auto arima 已移至 pmdarima github.com/tgsmith61591/pmdarima，可能是因为它不会与另一个更受欢迎的同名库发生冲突【参考方案6】：def evaluate_arima_model(X, arima_order):
# prepare training dataset
train_size = int(len(X) * 0.90)
train, test = X[0:train_size], X[train_size:]
history = [x for x in train]
# make predictions
predictions = list()
for t in range(len(test)):
model = ARIMA(history, order=arima_order)
model_fit = model.fit(disp=0)
yhat = model_fit.forecast()[0]
predictions.append(yhat)
history.append(test[t])
# calculate out of sample error
error = mean_squared_error(test, predictions)
return error
# evaluate combinations of p, d and q values for an ARIMA model
def evaluate_models(dataset, p_values, d_values, q_values):
dataset = dataset.astype('float32')
best_score, best_cfg = float("inf"), None
for p in p_values:
for d in d_values:
for q in q_values:
order = (p,d,q)
try:
mse = evaluate_arima_model(dataset, order)
if mse < best_score:
best_score, best_cfg = mse, order
print('ARIMA%s MSE=%.3f' % (order,mse))
except:
continue
print('Best ARIMA%s MSE=%.3f' % (best_cfg, best_score))
# load dataset
def parser(x):
return datetime.strptime('190'+x, '%Y-%m')
import datetime
p_values = [4,5,6,7,8]
d_values = [0,1,2]
q_values = [2,3,4,5,6]
warnings.filterwarnings("ignore")
evaluate_models(train, p_values, d_values, q_values)
这将为您提供 p、d、q 值，然后将这些值用于您的 ARIMA 模型
【讨论】：
旧线程，但这可能不是选择 d 的最佳方式。传统上，这是通过使用单位根测试（例如 KPSS 测试）来完成的。 请注意，此代码只是从here 复制而来，侵犯了该网站的版权。【参考方案7】：目前，我们可以直接使用 PyPI 中的pyramid-arima 包。
【讨论】：
【参考方案8】：在 conda 中，使用conda install -c saravji pmdarima 进行安装。
用户saravji已将其放入anaconda云中。
然后使用，
from pmdarima.arima import auto_arima
（注意名称pyramid-arima改为pmdarima）。
【讨论】：
以上是关于等价于 python 的 auto.arima()的主要内容，如果未能解决你的问题，请参考以下文章
利用Auto ARIMA构建高性能时间序列模型（附Python和R代码）