對於時間序列任意時刻的序列值Xt都是一個隨機變量，每一個隨機變量都會有均值和方差，記Xt的均值為ut，方差為任取t,s，定義序列的自協方差函數=E[(xt-ut)(Xs-us)]和自相關系數，之所以稱它們為自協方差函數和自相關系數，是因為它們衡量的是同一個事件在兩個不同時期（時刻t和時刻s）之間的相關程度，形象地講就是度量子集過去的行為對自己現在的影響。
如果事件序列在某一常數附近波動且波動范圍有限，即有常數均值和常數方差，並且延遲k期的序列變量的自協方差和自相關系數是相等的或者說延遲k期的序列變量之間的影響程度是一樣的，則稱為平穩序列。
平穩性的檢驗：
對序列平穩性的檢驗有兩種檢驗方法：一種是根據時序圖和自相關圖的特征做出判斷的圖檢驗，該方法操作簡單、應用廣泛；另一種是構造檢驗統計量進行的方法，目前最常用的方法是單位根檢驗。
l 時序圖檢驗：根據平穩時間序列的均值和方差都為常數的性質，平穩序列的時序圖顯示該序列值始終在一個常數附近隨機波動，而且波動的范圍有界；如果有明顯的趨勢性或周期性那它通常不是平穩序列。
l 自相關圖檢驗：平穩序列具有短期相關性，這個性質表明對平穩序列而言通常只有近期的序列值對現時值的影響比較明顯，間隔越遠的過去值對現時值的影響越小。隨着延遲期k的增加，平穩序列的自相關系數會比較快的衰減趨向於零，並在零附近隨機波動，而非平穩序列的自相關系數衰減的速度比較慢，這就是利用自相關圖進行平穩性檢驗的標准。
l 單位根檢驗：單位根檢驗是指檢驗序列中是否存在單位根，因為存在單位根就是非平穩性序列了。
from statsmodels.tsa.stattools import adfuller def check_stationary(ts,thresh=0.05): results = adfuller(ts) return results[1]<=thresh check_stationary(train_ts['co2'])
純隨機性檢驗：
如果一個序列是純隨機序列，那么它的序列值之間應該沒有任何關系，即滿足 ，這是一種理論上才會出現的理想狀態，實際上純隨機序列的樣本自相關系數不會絕對為零，但是很接近為零，並在零附近隨機波動。
純隨機性檢驗也稱為白噪聲檢驗，一般是構造檢驗統計量來檢驗序列的純隨機性，常用的檢驗統計量有Q統計量、LB統計量，由樣本各延遲期數的自相關系數可以計算得到檢驗統計量，然后計算出對應的p值，如果p值顯著大於顯著性水平，則表示該序列不能拒絕純隨機性的原假設，可以停止對該序列的分析。
from statsmodels.stats.diagnostic import acorr_ljungbox def check_stochastic(ts): p_value = acorr_ljungbox(ts,lags=1)[1] return p_value check_stochastic(train_ts['co2'])
差分運算：
差分運算具有強大的確定性信息提取能力，許多非平穩序列差分后會顯示出平穩序列的性質，這時稱這個非平穩序列為差分平穩序列。對差分平穩序列可以使用ARMA模型進行擬合。下圖顯示了差分前后的情況。
train_ts['diff1'] = train_ts['co2'].diff(1) train_ts.plot(figsize=(15,6))
對差分后的序列同樣進行平穩性與隨機性驗證。
check_stationary(train_ts['diff1'].dropna()) check_stochastic(train_ts['diff1'].dropna())
驗證通過后需要進行模型的定階，確定p、q：
第一種方法是人為識別的方法：將pacf、acf圖畫出，依據截尾還是拖尾進行查找
import statsmodels.api as sm fig = plt.figure(figsize=(12,8)) ax1 = fig.add_subplot(211) fig = sm.graphics.tsa.plot_acf(train_ts['diff1'].dropna(),lags=20,ax=ax1) ax2 = fig.add_subplot(212) fig = sm.graphics.tsa.plot_pacf(train_ts['diff1'].dropna(),lags=20,ax=ax2)
很多時候，並不能看出p、q的階數
第二種方法是相對最優模型識別：將p、d、q的可能值進行遍歷，一次帶入模型求最小的aic、bic、hqic
from statsmodels.tsa.arima_model import ARIMA ts = train_ts['co2'] import itertools p = d =q =range(0,3) pdq = list(itertools.product(p,d,q)) param_list = [] aic_list = [] bic_list = [] results_df = pd.DataFrame() for param in pdq: try: model = ARIMA(ts,order=param) results = model.fit() print 'ARIMA{} - AIC:{} - BIC:{} '.format(param,results.aic,results.bic) param_list.append(param) aic_list.append(results.aic) bic_list.append(results.bic) except: continue results_df['pdq'] = param_list results_df['aic'] = aic_list results_df['bic'] = bic_list results_df.loc[results_df['aic'].idxmin()]
取得p、q、d之后，帶入模型
best_result = ARIMA(ts,order=(2,1,2)).fit()
如果預測指定長期序列
pred_ts = best_result.predict(start='2000-01-01',end='2001-12-01',typ='levels')
但這樣的話，后期序列值會震盪開。這時可以采用滾動預測：
先創建要預測時間段的時間索引
test_range = pd.date_range(start='2000-01-01', end='2001-12-01', freq='MS')
利用已創建的時間索引，創建序列
pred_ts = pd.Series(index=test_range)
遍歷時間索引，對不斷迭代更新的序列，進行訓練、建模、預測
for cur_month in test_range: prev_month = cur_month - pd.Timedelta('1M') train_ts = proc_data_df.loc[:prev_month]['co2'] model = ARIMA(train_ts, order=(2, 1, 2)).fit() #pred = model.predict(start=cur_month, typ='levels') pred = model.forecast() pred_ts.loc[cur_month] = pred[0]
然后對測試序列、預測序列進行圖形上的比較
test_ts.plot(figsize=(15,8))
pred_ts.plot()
本站转载的文章为个人学习借鉴使用，本站对版权不负任何法律责任。如果侵犯了您的隐私权益，请联系我们删除。