trainSeting.head()
Out[36]: 
date
2017-10-01 126.4
2017-10-02  82.4
2017-10-03  78.1
2017-10-04  51.1
2017-10-05  90.9
Name: sales, dtype: float64

plt.plot(trainSeting)

'''進(jìn)行ADF檢驗(yàn)
adf_test的返回值
Test statistic：代表檢驗(yàn)統(tǒng)計(jì)量
p-value：代表p值檢驗(yàn)的概率
Lags used：使用的滯后k，autolag=AIC時(shí)會(huì)自動(dòng)選擇滯后
Number of Observations Used：樣本數(shù)量
Critical Value(5%) : 顯著性水平為5%的臨界值。
(1)假設(shè)是存在單位根，即不平穩(wěn)；
(2)顯著性水平，1%：嚴(yán)格拒絕原假設(shè)；5%：拒絕原假設(shè)，10%類(lèi)推。
(3)看P值和顯著性水平a的大小，p值越小，小于顯著性水平的話(huà)，就拒絕原假設(shè)，認(rèn)為序列是平穩(wěn)的；大于的話(huà)，不能拒絕，認(rèn)為是不平穩(wěn)的
(4)看檢驗(yàn)統(tǒng)計(jì)量和臨界值，檢驗(yàn)統(tǒng)計(jì)量小于臨界值的話(huà)，就拒絕原假設(shè)，認(rèn)為序列是平穩(wěn)的；大于的話(huà)，不能拒絕，認(rèn)為是不平穩(wěn)的
'''
#滾動(dòng)統(tǒng)計(jì)
def rolling_statistics(timeseries):
 #Determing rolling statistics
 rolmean = pd.rolling_mean(timeseries, window=12)
 rolstd = pd.rolling_std(timeseries, window=12)
 
 #Plot rolling statistics:
 orig = plt.plot(timeseries, color='blue',label='Original')
 mean = plt.plot(rolmean, color='red', label='Rolling Mean')
 std = plt.plot(rolstd, color='black', label = 'Rolling Std')
 plt.legend(loc='best')
 plt.title('Rolling Mean & Standard Deviation')
 plt.show(block=False)
 
##ADF檢驗(yàn)
from statsmodels.tsa.stattools import adfuller
def adf_test(timeseries):
 rolling_statistics(timeseries)#繪圖
 print ('Results of Augment Dickey-Fuller Test:')
 dftest = adfuller(timeseries, autolag='AIC')
 dfoutput = pd.Series(dftest[0:4], index=['Test Statistic','p-value','#Lags Used','Number of Observations Used'])
 for key,value in dftest[4].items():
  dfoutput['Critical Value (%s)'%key] = value #增加后面的顯著性水平的臨界值
 print (dfoutput)
 
adf_test(trainSeting) #從結(jié)果中可以看到p值為0.1097>0.1，不能拒絕H0,認(rèn)為該序列不是平穩(wěn)序列

Results of Augment Dickey-Fuller Test:
Test Statistic    -5.718539e+00
p-value      7.028398e-07
#Lags Used      0.000000e+00
Number of Observations Used 6.200000e+01
Critical Value (1%)   -3.540523e+00
Critical Value (5%)   -2.909427e+00
Critical Value (10%)   -2.592314e+00
dtype: float64

'''acorr_ljungbox(x, lags=None, boxpierce=False)函數(shù)檢驗(yàn)無(wú)自相關(guān)
lags為延遲期數(shù)，如果為整數(shù)，則是包含在內(nèi)的延遲期數(shù)，如果是一個(gè)列表或數(shù)組，那么所有時(shí)滯都包含在列表中最大的時(shí)滯中
boxpierce為T(mén)rue時(shí)表示除開(kāi)返回LB統(tǒng)計(jì)量還會(huì)返回Box和Pierce的Q統(tǒng)計(jì)量
返回值：
lbvalue:測(cè)試的統(tǒng)計(jì)量
pvalue:基于卡方分布的p統(tǒng)計(jì)量
bpvalue:((optionsal), float or array) – test statistic for Box-Pierce test
bppvalue:((optional), float or array) – p-value based for Box-Pierce test on chi-square distribution
'''
from statsmodels.stats.diagnostic import acorr_ljungbox
def test_stochastic(ts,lag):
 p_value = acorr_ljungbox(ts, lags=lag) #lags可自定義
 return p_value

####自相關(guān)圖ACF和偏相關(guān)圖PACF
import statsmodels.api as sm
def acf_pacf_plot(ts_log_diff):
 sm.graphics.tsa.plot_acf(ts_log_diff,lags=40) #ARIMA,q
 sm.graphics.tsa.plot_pacf(ts_log_diff,lags=40) #ARIMA,p
 
acf_pacf_plot(trainSeting) #查看數(shù)據(jù)的自相關(guān)圖和偏自相關(guān)圖

##借助AIC、BIC統(tǒng)計(jì)量自動(dòng)確定
from statsmodels.tsa.arima_model import ARMA
def proper_model(data_ts, maxLag): 
 init_bic = float("inf")
 init_p = 0
 init_q = 0
 init_properModel = None
 for p in np.arange(maxLag):
  for q in np.arange(maxLag):
   model = ARMA(data_ts, order=(p, q))
   try:
    results_ARMA = model.fit(disp=-1, method='css')
   except:
    continue
   bic = results_ARMA.bic
   if bic < init_bic:
    init_p = p
    init_q = q
    init_properModel = results_ARMA
    init_bic = bic
 return init_bic, init_p, init_q, init_properModel
 
proper_model(trainSeting,40)

#在statsmodels包里還有更直接的函數(shù)：
import statsmodels.tsa.stattools as st
order = st.arma_order_select_ic(ts_log_diff2,max_ar=5,max_ma=5,ic=['aic', 'bic', 'hqic'])
order.bic_min_order
'''
我們常用的是AIC準(zhǔn)則，AIC鼓勵(lì)數(shù)據(jù)擬合的優(yōu)良性但是盡量避免出現(xiàn)過(guò)度擬合(Overfitting)的情況。所以?xún)?yōu)先考慮的模型應(yīng)是AIC值最小的那一個(gè)模型。
為了控制計(jì)算量，我們限制AR最大階不超過(guò)5，MA最大階不超過(guò)5。 但是這樣帶來(lái)的壞處是可能為局部最優(yōu)。
timeseries是待輸入的時(shí)間序列，是pandas.Series類(lèi)型，max_ar、max_ma是p、q值的最大備選值。
order.bic_min_order返回以BIC準(zhǔn)則確定的階數(shù)，是一個(gè)tuple類(lèi)型

################################模型######################################
# AR模型，q=0
#RSS是殘差平方和
# disp為-1代表不輸出收斂過(guò)程的信息，True代表輸出
from statsmodels.tsa.arima_model import ARIMA
model = ARIMA(trainSeting,order=(1,0,0)) #第二個(gè)參數(shù)代表使用了二階差分
results_AR = model.fit(disp=-1)
plt.plot(trainSeting)
plt.plot(results_AR.fittedvalues, color='red') #紅色線(xiàn)代表預(yù)測(cè)值
plt.title('RSS:%.4f' % sum((results_AR.fittedvalues-trainSeting)**2))#殘差平方和

############################預(yù)測(cè)未來(lái)走勢(shì)##########################################
# forecast方法會(huì)自動(dòng)進(jìn)行差分還原，當(dāng)然僅限于支持的1階和2階差分
forecast_n = 12 #預(yù)測(cè)未來(lái)12個(gè)天走勢(shì)
forecast_AR = results_AR.forecast(forecast_n)
forecast_AR = forecast_AR[0]
print (forecast_AR)

##將預(yù)測(cè)的數(shù)據(jù)和原來(lái)的數(shù)據(jù)繪制在一起，為了實(shí)現(xiàn)這一目的，我們需要增加數(shù)據(jù)索引，使用開(kāi)源庫(kù)arrow:
import arrow
def get_date_range(start, limit, level='day',format='YYYY-MM-DD'):
 start = arrow.get(start, format) 
 result=(list(map(lambda dt: dt.format(format) , arrow.Arrow.range(level, start,limit=limit))))
 dateparse2 = lambda dates:pd.datetime.strptime(dates,'%Y-%m-%d')
 return map(dateparse2, result)
 
# 預(yù)測(cè)從2017-12-03開(kāi)始，也就是我們訓(xùn)練數(shù)據(jù)最后一個(gè)數(shù)據(jù)的后一個(gè)日期
new_index = get_date_range('2017-12-03', forecast_n)
forecast_ARIMA_log = pd.Series(forecast_AR, copy=True, index=new_index)
print (forecast_ARIMA_log.head())
##繪圖如下
plt.plot(trainSeting,label='Original',color='blue')
plt.plot(forecast_ARIMA_log, label='Forcast',color='red')
plt.legend(loc='best')
plt.title('forecast')