聯絡方式

• Author: sunhailin-Leo
• E-mail: [email protected]
• 程式碼: 倉庫地址 -> github.com/sunhailin-L…

簡介

• 第一篇文章傳送門: Python爬取中國銀行外匯牌價(爬蟲 + PyFlux簡單預測分析)--(一)

• 如果資料還沒獲取到的話請移步看第一篇文章(本篇文章預設資料已儲存在資料庫中)

步入正題

• 上一篇文章末尾部分在講使用Pyflux這個庫對資料進行分析預測
• 總結兩個方面:
  • 優點:
    • Pyflux模型文件"一針見血"(建立在對時序分析有一定基礎的人, 能看懂部分核心公式)
  • 缺點:
    • 提供少量的資料分析API, 不像statsmodels提供了例如殘差分析等方法進行模型驗證調優的方法
• 本文將使用statsmodels對此前的資料進行分析。

資料前期準備

df['查詢時間'] = df['查詢時間'].apply(lambda x: x[:-9])
df['查詢時間'] = pd.to_datetime(df['查詢時間'], format="%Y-%m-%d")
df = df.groupby('查詢時間')['現匯賣出價'].mean()
df = df.to_frame()
print(df)
複製程式碼
 

資料平穩性校驗

• 直接上程式碼(差分畫圖看資料平穩性)
• 一共測試了1階和2階(不建議使用高階資料, 容易資料造成破壞)

# 差分圖
fig = plt.figure(figsize=(12, 8))
ax1 = fig.add_subplot(111)
# 裡面的1代表了差分階數
diff1 = df.diff(1)
diff1.plot(ax=ax1)
plt.show()
# fig.savefig('./picture/1.jpg')
複製程式碼

• 一階差分結果

• 二階差分結果

• 結論：由觀察可得出一階差分波動相對較小，較為平穩 (之後建模的資料使用1階差分的結果進行建模
  )

自相關圖和偏自相關圖

• 介紹一些自相關性和非自相關性:

  • 自相關性(ACF): 自相關性是指隨機誤差項的各期望值之間存在著相關關係，稱隨機誤差項之間存在自相關性（autocorrelation）或序列相關
  • 偏自相關性(PACF): 偏自相關是剔除干擾後時間序列觀察與先前時間步長時間序列觀察之間關係的總結（Partial autocorrelation）

• 接下來選擇40個滯後做自相關和偏自相關的分析

# statmodels 自相關圖 偏相關圖
diff1 = diff1.dropna()
fig = plt.figure(figsize=(12, 8))
ax1 = fig.add_subplot(211)
fig = sm.graphics.tsa.plot_acf(diff1, lags=40, ax=ax1)
ax2 = fig.add_subplot(212)
fig = sm.graphics.tsa.plot_pacf(diff1, lags=40, ax=ax2)
fig.show()
# fig.savefig('./picture/acf_pacf.jpg')
複製程式碼

• 接下來由看圖觀察定義ARMA模型的p(自相關係數), q(偏自相關係數)引數:

  • 自相關圖顯示滯後有三個階超出了置信邊界
  • 偏相關圖顯示在滯後2階時的偏自相關係數超出了置信邊界(其實波動還是沒有完全消除, 只是2階的時候比較突出)

• 可以嘗試如下p, q引數的模型, 並輸出AIC,BIC,HQIC的結果

  • AIC:赤池資訊量 akaike information criterion --> -2ln(L) + 2k
  • BIC:貝葉斯資訊量 bayesian information criterion --> -2ln(L) + ln(n) * k
  • HQIC:HQ準則 --> -2ln(L) + ln(ln(n)) * k
  • L是在該模型下的最大似然，n是資料數量，k是模型的變數個數。

• 嘗試結果如下:

  • p=2, q=0 --> -194.77200890459767 -179.81283725588074 -188.79262385129292
  • p=2, q=1 --> -197.01722373566554 -178.31825917476937 -189.54299241903462
  • p=2, q=2 --> -195.11076756537022 -172.6720100922948 -186.1416899854131
  • p=3, q=2 --> -201.37730533090803 -175.19875494565338 -190.91338148762475

• 模型選擇:

  • 在上述的結果看來選擇p=2, q=2為最優辦法, 雖然不一定是最適合, 但是在差分和根據自相關和偏自相關圖的結果看來, 可以嘗試使用(2, 2)引數進行建模

建模

• 建模的程式碼很短,就不做過多敘述,看下述程式碼

arma_mod22 = sm.tsa.ARMA(diff1, (3, 2)).fit()
# 輸出AIC,BIC和HQ準則結果
print(arma_mod22.aic, arma_mod22.bic, arma_mod22.hqic)
複製程式碼

模型驗證

1. 殘差驗證

# 殘差(輸出形式為DataFrame)
resid = arma_mod22.resid
# 殘差的ACF和PACF圖
fig = plt.figure(figsize=(12, 8))
ax1 = fig.add_subplot(211)
fig = sm.graphics.tsa.plot_acf(resid.values.squeeze(), lags=40, ax=ax1)
ax2 = fig.add_subplot(212)
fig = sm.graphics.tsa.plot_pacf(resid, lags=40, ax=ax2)
fig.show()
# fig.savefig('./picture/resid_acf_pacf.jpg')
複製程式碼

• 根據上圖可以發現雖然少數部分的階數下還有出現了超過置信區間的問題，但是總體看來序列殘差基本為白噪聲。

2. D-W驗證

• D-W驗證是檢驗自相關性的方法(只適用於檢測1階自相關性, 高階需要另闢蹊徑)

# 殘差D-W檢驗
resid_dw_result = sm.stats.durbin_watson(arma_mod22.resid.values)
# 1.9933441709003574 接近於 2 所以殘差序列不存在自相關性。
print(resid_dw_result)
複製程式碼

• 結果等於1.9933441709003574
• 結論：殘差序列不存在自相關性
• 理由：DW結果接近2的時候序列幾乎不存在自相關性(參考文章: www.doc88.com/p-543541848…)

3. 正態分佈

# 正態校驗 -> 基本符合正態分佈
fig = plt.figure(figsize=(12, 8))
ax = fig.add_subplot(111)
fig = sm.qqplot(resid, line='q', ax=ax, fit=True)
fig.show()
# fig.savefig('./picture/normal_distribution.jpg')
複製程式碼

• 除了右上角的部分異常點(這幾個點的資料很直觀的在接下來的Q檢驗中顯現出來)
• 由圖得出，絕大部分資料還是符合正態分佈

4. Q檢驗

# 殘差序列Ljung-Box檢驗，也叫Q檢驗
r, q, p = sm.tsa.acf(resid.values.squeeze(), qstat=True)
data = np.c_[range(1, 41), r[1:], q, p]
table = pd.DataFrame(data, columns=['lag', "AC", "Q", "Prob(>Q)"])
temp_df = table.set_index('lag')
print(temp_df)
# Prob(>Q)的最小值: 0.025734615668093132, 最大值: 0.9874705305611844, 均值: 0.2782013984159408
prob_q_min = temp_df['Prob(>Q)'].min()
prob_q_max = temp_df['Prob(>Q)'].max()
prob_q_mean = temp_df['Prob(>Q)'].mean()
print("Prob(>Q)的最小值: {}, 最大值: {}, 均值: {}".format(prob_q_min, prob_q_max, prob_q_mean))
複製程式碼

• 結果：Prob(>Q)的最小值: 0.025734615668093132, 最大值: 0.9874705305611844, 均值: 0.2782013984159408
• 由上圖可以看到最後10個數據的值幾乎在0.05左右，可以對應回正態分佈的那些離散點的值
• 根據Q校驗的方法, 在95%的置信區間內, 當Prob大於0.05當前殘差序列不存在自相關性, 因此滯後40個數據的模型屬於基本不存在自相關性。

5. 預測(預測從11月9日到11月14日的匯率變動)

• 由於此前將資料進行了1階差分，因此預測結果的資料也是1階差分的預測值

# 預測
predict_data = arma_mod22.predict('2018-11-09', '2018-11-14', dynamic=True)
# 畫預測圖
fig, ax = plt.subplots(figsize=(12, 8))
ax = diff1.ix['2018-01-01':].plot(ax=ax)
fig = arma_mod22.plot_predict('2018-11-09', '2018-11-14', dynamic=True, ax=ax, plot_insample=False)
fig.show()
# 結果預測
last_data = df['現匯賣出價'].values[-1]
# 還原結果
predict_data_list = predict_data.values.tolist()
restore_list = []
for d in predict_data_list:
    last_data = last_data + d
    restore_list.append(last_data)
predict_data = pd.DataFrame(restore_list, index=predict_data.index, columns=['現匯賣出價'])
print(predict_data)
複製程式碼

• 注:以上預測資料均為當天的匯率的均值

總結

• 以上的分析、建模、校驗、預測的流程只是簡單的一個流程.預測模型準確性還需要不斷驗證, 以上的關鍵程式碼僅提供了一套可以參考的流程.
• 時序資料預測與分析還有待學習，針對不同型別的資料和多變數的資料有著不同的模型, 以上步驟僅供參考.
• 目前使用過時序預測的庫有：statsmodels和pyflux, 之後考慮使用機器學習或者深度學習的方法進行預測.