時間序列的理論
u 平穩時間序列
時間序列平穩性定義：
平穩時間序列分為：自回歸模型，滑動平均模型，自回歸滑動平均模型
自回歸模型：當前值由前p期值決定
滑動平均模型：
自回歸滑動平均模型：
根據模型的自相關圖，AR(p)模型的自相關系數隨着延遲階數的增加逐漸遞減，呈現拖尾狀態，而偏自相關系數隨着延遲階數的增加迅速減到0，呈現截尾狀態。MA(q)模型與AR(p)模型相反。ARMA模型自相關和偏自相關圖均是拖尾的。
模型的拖尾性和截尾性：
u 非平穩時間序列：
k階差分
一階差分：相距一期的兩個序列之間的減法運算稱為一階差分運算
二階差分：對一階差分后序列再進行一次一階差分運算稱為二階差分
k步差分：
相距k期的兩個序列值之間的減法運算稱為k步差分運算
差分方式的選擇：
實踐中，我們會根據序列的不同特點選擇合適的差分方式，常見情況有如下三種：
（1）具有顯著線性趨勢的序列，通常一階差分可以實現差分后平穩。
（2）具有曲線趨勢的序列，通常低階（二階或者三階）差分可以實現差分后平穩。
（3）具有固定周期的序列，通常進行步長為周期長度的差分運算，可以實現差分后平穩。
平穩時間序列的建模
u 平穩序列建模步驟
假如某個觀察值序列通過序列預處理可以判定為平穩非白噪聲序列，就可以利用ARMA模型對該序列進行建模。建模的基本步驟如下：
（1）求出該觀察值序列的樣本自相關系數（ACF）和樣本偏自相關系數（PACF）的值。
（2）根據樣本自相關系數和偏自相關系數的性質，選擇適當的ARMA（p，q）模型進行擬合。
（3）估計模型中位置參數的值。
（4）檢驗模型的有效性。如果模型不通過檢驗，轉向步驟（2），重新選擇模型再擬合。
（5）模型優化。如果擬合模型通過檢驗，仍然轉向不走（2），充分考慮各種情況，建立多個擬合模型，從所有通過檢驗的擬合模型中選擇最優模型。
（6）利用擬合模型，預測序列的將來走勢。
u 代碼實現：
#導入數據
discfile = 'd:/data/arima_data.xls'
forecastnum = 5
#讀取數據，指定日期列為指標，Pandas自動將“日期”列識別為Datetime格式
data = pd.read_excel(discfile, index_col = u'日期')
data = pd.DataFrame(data,dtype=np.float64)
data
#時序圖
plt.rcParams['font.sans-serif'] = ['SimHei'] #用來正常顯示中文標簽
plt.rcParams['axes.unicode_minus'] = False #用來正常顯示負號
data.plot()
plt.show()
#自相關圖
from statsmodels.graphics.tsaplots import plot_acf
plot_acf(data).show()
>>>觀察自相關圖可以看出是非平穩序列
#平穩性檢測
from statsmodels.tsa.stattools import adfuller as ADF
print( ADF(data[u'銷量']))
>>>
(1.813771015094526, 0.9983759421514264, 10L, 26L, {'5%': -2.981246804733728, '1%': -3.7112123008648155, '10%': -2.6300945562130176}, 299.46989866024177)返回值依次為adf、pvalue、usedlag、nobs、critical values、icbest、regresults、resstore
#由上面時間序列平穩性檢驗結果可知是非平穩時間序列
如何確定該序列能否平穩呢？主要看：
- 1%、%5、%10是不同程度拒絕原假設（原假設是存在單位根，即非平穩）的統計值，與ADF進行比較，ADF小於1%就可以極顯著的拒絕原假設，即說明數據是平穩的。
- P-value是否非常接近0.05，小於0.05，是拒絕原假設，即時間序列是平穩的
#差分后的結果
D_data = data.diff().dropna()
D_data.columns = [u'銷量差分']
D_data.plot() #時序圖
plt.show()
#自相關圖
plot_acf(D_data).show()
由自相關圖可以知道時間序列是1階截尾
#偏自相關圖
from statsmodels.graphics.tsaplots import plot_pacf
plot_pacf(D_data).show()
由偏自相關圖可以知道序列是拖尾的
ADF(D_data[u'銷量差分'])#平穩性檢測
#白噪聲檢驗
from statsmodels.stats.diagnostic import acorr_ljungbox
acorr_ljungbox(D_data, lags=1) #返回統計量和pvalue值
>>>
(array([11.30402222]), array([0.00077339])) #pvalue值遠小於0.05，序列是非白噪聲的
如何判斷是否是白噪聲序列？
原假設是序列是白噪聲序列，我們只看pvalue值，如果值小於0.05，就拒絕原假設，即時間序列是非白噪聲的
總結：根據上面的自相關圖和偏自相關圖，能夠大概估計出模型是ARIMA(0,1,1)，這種觀察圖的方法誤差較大
下面使用更科學的方法：根據信息准則判斷
#定階
from statsmodels.tsa.arima_model import ARIMA
pmax = int(len(D_data)/10) #一般階數不超過length/10
qmax = int(len(D_data)/10) #一般階數不超過length/10
bic_matrix = [] #bic矩陣
for p in range(pmax+1):
tmp = []
for q in range(qmax+1):
try: #存在部分報錯，所以用try來跳過報錯。
tmp.append(ARIMA(data, (p,1,q)).fit().bic)
except:
tmp.append(None)
bic_matrix.append(tmp) # bic_matrix是一個二維的矩陣
#從中可以找出最小值
bic_matrix = pd.DataFrame(bic_matrix)
p,q = bic_matrix.stack().idxmin() #先用stack展平，然后用idxmin找出最小值位置。
#打印最小值對應的p值和q值
print(u'BIC最小的p值和q值為：%s、%s' %(p,q))
#建立ARIMA(0, 1, 1)模型
model = ARIMA(data, (0,1,1)).fit()
#給出一份模型報告
model.summary()
#作為期5天的預測，返回預測結果、標准誤差、置信區間。
model.forecast(5)