auto.arima()相当于python
我试图预测使用ARMA ARIMA模型的每周销售量。 我无法find调整statsmodels
的顺序(p,d,q)的statsmodels
。 目前R有一个函数auto.arima()
,它将调整(p,d,q)参数。
我如何去为我的模型select正确的顺序? Python中有没有用于此目的的库?
您可以实施多种方法:
-
ARIMAResults
包括aic
和bic
。 根据他们的定义(见这里和这里 ),这些标准惩罚模型中参数的数量。 所以你可以使用这些数字来比较模型。 另外scipy有在指定的参数空间上进行网格search的optimize.brute
。 所以这样的工作stream程应该工作: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
。 -
您可以从
ARIMAResults
获得pvalues
。 因此,一种前向algorithm易于实现,其中模型的程度在增加的参数获得最小的p值的维度上增加。 -
使用
ARIMAResults.predict
来交叉validation替代模型。 最好的办法是保持时间序列的尾部(比如最近5%的数据)不在样本中,用这些点来获得拟合模型的testing误差 。
解决scheme
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('ARIMA{}x{}12 - 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://github.com/decisionstats/pythonfordatascience/blob/master/time%2Bseries%20(1).ipynb
我写了这些效用函数来直接计算pdq值get_PDQ_parallel需要三个input数据,它们是以timestamp(datetime)作为索引的序列。 n_jobs将提供多个并行处理器。 输出将是具有aic和bic值的dataframe,其中索引p中的order =(P,D,Q),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'])
其实
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)