摘要：作者天瓊，數(shù)據(jù)游戲優(yōu)勝隊(duì)伍成員介紹本文整理記錄了參與的一次小型數(shù)據(jù)分析競賽數(shù)據(jù)游戲，競賽目標(biāo)是預(yù)測年月日股閉市時(shí)招商銀行的股價(jià)。日發(fā)現(xiàn)的數(shù)據(jù)有錯(cuò)誤，需要手工矯正日該數(shù)據(jù)恢復(fù)正常。而函數(shù)，是對樣本外的數(shù)據(jù)進(jìn)行預(yù)測。

作者：天瓊，「數(shù)據(jù)游戲」優(yōu)勝隊(duì)伍成員

本文整理記錄了參與的一次小型數(shù)據(jù)分析競賽「數(shù)據(jù)游戲」，競賽目標(biāo)是預(yù)測2019年5月15日A股閉市時(shí)招商銀行600036的股價(jià)。

主要思路是利用ARIMA算法做時(shí)間序列預(yù)測。

使用的數(shù)據(jù)是公開的數(shù)據(jù)集 tushare。

拿到題目和數(shù)據(jù)之后，首先結(jié)合既往經(jīng)歷，覺得想要預(yù)測準(zhǔn)股價(jià)，本身是一個(gè)不可能的事情，尤其是A股。

影響股價(jià)的因素非常復(fù)雜而且不透明，以及金融投資領(lǐng)域具有的反身性理論，使得這次預(yù)測更多偏向于實(shí)驗(yàn)性質(zhì)，同時(shí)對競爭結(jié)果不要有過高的期望。

預(yù)測得準(zhǔn)，是你的運(yùn)氣；預(yù)測的爛，也不會(huì)影響你從中學(xué)到什么。學(xué)習(xí)第一，比賽第二吧。

鑒于以上，本次預(yù)測只使用了close的時(shí)間序列。更多的數(shù)據(jù)其實(shí)并沒有什么用。

獲取時(shí)間序列數(shù)據(jù) ；

觀察數(shù)據(jù)是否為平穩(wěn)時(shí)間序列；

對于非平穩(wěn)時(shí)間序列要先進(jìn)行d階差分運(yùn)算，轉(zhuǎn)化為平穩(wěn)時(shí)間序列；

對處理得到的平穩(wěn)時(shí)間序列，求它的階數(shù)p,q;

根據(jù)ARIMA算法建模，ARIMA(data, order=(p, d, q))

模型檢驗(yàn)和調(diào)優(yōu)

預(yù)測

初學(xué)的小伙伴們可能對這個(gè)過程并不熟悉，沒關(guān)系，先背下來。

對Python不熟悉的小伙伴們， 我給大家總結(jié)了幾句車轱轆話，大家先強(qiáng)行了解下。

獲取時(shí)間序列數(shù)據(jù):

data = pd.read_excel(‘600036.xlsx’, index=None)
train = data[‘close’]

觀察數(shù)據(jù)是否為平穩(wěn)時(shí)間序列；對于非平穩(wěn)時(shí)間序列要先進(jìn)行d階差分運(yùn)算，轉(zhuǎn)化為平穩(wěn)時(shí)間序列；

adf_data = sts.adfuller(train)
diff = train.diff(1)

對處理得到的平穩(wěn)時(shí)間序列，求它的階數(shù)p, q

sm.tsa.arma_order_select_ic(train, max_ar=8, max_ma=8, ic=[‘a(chǎn)ic’, ‘bic’, ‘hqic’])

根據(jù)ARIMA算法建模

ARMAModel = sm.tsa.ARIMA(train, order=(4,1,2)).fit()

模型檢驗(yàn)和調(diào)優(yōu)

train_shift = train.shift(1)
pred_recover = predicts.add(train_shift)
np.sqrt( sum( (pred_recover -train) ** 2)/train.size )

預(yù)測

f = ARMAModel.forecast(3)

以上，是本文的核心代碼，大家如果一時(shí)看不懂，可以跳過。
可以看看下面更詳細(xì)的步驟。

# 導(dǎo)入必須的模塊
import tushare as ts  #使用的公開的數(shù)據(jù)
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import statsmodels.api as sm
import statsmodels.tsa.stattools as sts

import wconfig as wc  #自定義打印輸出控制模塊，與本次競賽無關(guān)
wc.dispy() 

# 導(dǎo)入數(shù)據(jù)，從去年1月開始即可
data = ts.get_hist_data("600036", start="2018-01-08", end="2019-05-10").sort_index()
#data = ts.get_hist_data("600036", start="2018-04-01").sort_index().reset_index()

# 數(shù)據(jù)安全和源數(shù)據(jù)備份
#data.to_excel("600036.xlsx")

# 源數(shù)據(jù)可靠性檢查。12日發(fā)現(xiàn)tushare的數(shù)據(jù)有錯(cuò)誤，需要手工矯正（13日該數(shù)據(jù)恢復(fù)正常）。
print("tushare中5月10日close股價(jià)為: %.2f元，與實(shí)際不符！" % (data["close"]["2019-05-10"]))
data["close"]["2019-05-10"]=33.61
print("當(dāng)日實(shí)際收盤價(jià)應(yīng)為: %.2f元" % data["close"]["2019-05-10"])

# 只取 close 字段作為訓(xùn)練數(shù)據(jù)
train = data["close"]
train.index = pd.to_datetime(train.index)  # 將字符串索引轉(zhuǎn)換成時(shí)間索引
train.tail()
train.tail().index

tushare中5月10日close股價(jià)為: 33.48元，與實(shí)際不符！
當(dāng)日實(shí)際收盤價(jià)應(yīng)為: 33.61元

ARIMA算法要求時(shí)間序列穩(wěn)定，所以在建模之前，要先檢驗(yàn)時(shí)間序列的穩(wěn)定性。
adfuller就是用來干這個(gè)的。
adfuller 全稱 Augmented Dickey–Fuller test， 即擴(kuò)展迪基-福勒檢驗(yàn)，用來測試平穩(wěn)性。

def tagADF(t):
    result = pd.DataFrame(index=["Test Statistic Value", "p-value", "Lags Used", 
                                 "Number of Observations Used", 
                                 "Critical Value(1%)", "Critical Value(5%)", "Critical Value(10%)"],
                          columns=["value"]
    )
    result["value"]["Test Statistic Value"]=t[0]
    result["value"]["p-value"]=t[1]
    result["value"]["Lags Used"]=t[2]
    result["value"]["Number of Observations Used"] = t[3]
    result["value"]["Critical Value(1%)"]=t[4]["1%"]
    result["value"]["Critical Value(5%)"]=t[4]["5%"]
    result["value"]["Critical Value(10%)"]=t[4]["10%"]
    print("t is:", t)
    return result

adfuller檢驗(yàn)是檢查時(shí)間序列平穩(wěn)性的統(tǒng)計(jì)測試之一。 這里的零假設(shè)是：序列 train 是非平穩(wěn)的。
測試結(jié)果包括測試統(tǒng)計(jì)和差異置信水平的一些關(guān)鍵值。 如果測試結(jié)果中的 P-value 小于臨界值，我們可以拒絕原假設(shè)并說該序列是平穩(wěn)的。

adf_data = sts.adfuller(train)
tagADF(adf_data)

檢驗(yàn)結(jié)果顯示，p-value=0.414, 遠(yuǎn)遠(yuǎn)大于5%的臨界值，說明零假設(shè)是成立的，即序列 train 是非平穩(wěn)的。

t is: (-1.7325346908056185, 0.4144323576685054, 0, 322, {‘1%’: -3.4508226600665037, ‘5%’: -2.870558121868621, ‘10%’: -2.571574731684734}, 523.9067372199033)

為了讓時(shí)間序列平穩(wěn)，需要對 train 序列做差分運(yùn)算：

# df.diff 差分運(yùn)算，默認(rèn)是后一行減前一行
# http://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.DataFrame.diff.html?highlight=diff#pandas.DataFrame.diff
diff = train.diff(1).dropna()  # df.dropna 刪除有缺失值的行
plt.figure(figsize=(11,6))  # 指定顯示大小
plt.plot(diff, label="Diff")  # 繪制趨勢圖
plt.legend(loc=0)  # 顯示圖例，loc指定圖例位置，0為最佳位置。

關(guān)于時(shí)間序列穩(wěn)定性的判斷標(biāo)準(zhǔn)，可參考這篇博客：https://blog.csdn.net/u012735...

檢驗(yàn)差分后數(shù)據(jù)的平穩(wěn)性，和第一次驗(yàn)證方法相同

adf_Data1 = sts.adfuller(diff)
tagADF(adf_Data1)  # p-value很小，零假設(shè)不成立，因此，diff數(shù)據(jù)序列符合平穩(wěn)性要求。

在這里，先簡單解釋下自己對 ARIMA 算法的理解。不正確的地方請大家指點(diǎn)。

ARIMA算法認(rèn)為時(shí)間序列上 t 時(shí)刻的值由2部分構(gòu)成，第一部分是由之前p項(xiàng)歷史值決定的，比如15日的收盤價(jià)是12日，13日，14日的收盤價(jià)的線性回歸，用AR§表示。

但是這個(gè)線性回歸的輸出值肯定和15日的實(shí)際收盤價(jià)有一個(gè)誤差，假設(shè)是e。

所以第二部分就是如果表示這個(gè) e 。

這個(gè)e 可以認(rèn)為是之前q項(xiàng)誤差的線性回歸，用MA(q)表示。

兩部分合起來就是ARMA。

然后為了讓時(shí)間序列平穩(wěn)，再加個(gè)d 階差分操作功能，就變成了ARIMA算法。

所以在使用 ARIMA 算法之前，需要先確定（p, d, q）的值。

確定ARIMA的階數(shù)p,q

# ARMA(p,q)是AR(p)和MA(q)模型的組合，關(guān)于p和q的選擇，一種方法是觀察自相關(guān)圖ACF和偏相關(guān)圖PACF, 
# 另一種方法是通過借助AIC、BIC等統(tǒng)計(jì)量自動(dòng)確定。 
ic = sm.tsa.arma_order_select_ic(
    train, 
    max_ar=8, 
    max_ma=8, 
    ic=["aic", "bic", "hqic"]
)
ic

注意ARIMA的參數(shù)中，輸入數(shù)據(jù) 應(yīng)該是原始數(shù)據(jù) train，ARIMA 會(huì)根據(jù) d 的值，對原始數(shù)據(jù)做 d 階差分運(yùn)算。
d的含義是，輸入序列需要先經(jīng)過一個(gè)d階的差分，變成一個(gè)平穩(wěn)序列后才能進(jìn)行數(shù)據(jù)擬合。

ARMAModel = sm.tsa.ARIMA(train, order=(4,1,2)).fit()  # order=(p,d,q)
# fittedvalues和diff對比
plt.figure(figsize=(11, 6))
plt.plot(diff, "r", label="Orig")
plt.plot(ARMAModel.fittedvalues, "g",label="ARMA Model")
plt.legend()

# 樣本內(nèi)預(yù)測
predicts = ARMAModel.predict()

# 因?yàn)轭A(yù)測數(shù)據(jù)是根據(jù)差分值算的，所以要對它一階差分還原
train_shift = train.shift(1)  # shift是指series往后平移1個(gè)時(shí)刻
pred_recover = predicts.add(train_shift).dropna()  #這里add是指兩列相加，按index對齊

# 模型評(píng)價(jià)指標(biāo) 1：計(jì)算 score
delta = ARMAModel.fittedvalues - diff
score = 1 - delta.var()/train.var()
print("score:
", score)

# 模型評(píng)價(jià)指標(biāo) 2：使用均方根誤差（RMSE）來評(píng)估模型樣本內(nèi)擬合的好壞。
#利用該準(zhǔn)則進(jìn)行判別時(shí)，需要剔除“非預(yù)測”數(shù)據(jù)的影響。
train_vs = train[pred_recover.index]  # 過濾沒有預(yù)測的記錄
plt.figure(figsize=(11, 6))
train_vs.plot(label="Original")
pred_recover.plot(label="Predict")
plt.legend(loc="best")
plt.title("RMSE: %.4f"% np.sqrt(sum((pred_recover-train_vs)**2)/train_vs.size))
plt.show()

# 局部數(shù)據(jù)觀察
train_t = train_vs.tail(15)
pred_t = pred_recover.tail(15)
plt.figure(figsize=(11, 6))
train_t.plot(label="Original")
pred_t.plot(label="Predict")
plt.legend(loc="best")
plt.title("RMSE: %.4f"% np.sqrt(sum((pred_t-train_t)**2)/train_t.size))
plt.show()

使用 forecast 對樣本外的時(shí)間序列進(jìn)行預(yù)測。

關(guān)于foreast和predict的區(qū)別：
predict 可以對樣本內(nèi)和樣本外的進(jìn)行預(yù)測，結(jié)果是一樣的。
舉例說明：forecast（10），表示對未來10個(gè)點(diǎn)進(jìn)行預(yù)測，但是可以用model.fittedvalues查看樣本內(nèi)點(diǎn)的擬合值；
而predict(start,end)里面的參數(shù)0表示樣本內(nèi)的第一個(gè)數(shù)，以此類推。
如果想要預(yù)測樣本外的數(shù)，需要將start設(shè)置為len(data)+1,即數(shù)據(jù)長度+1，才表示預(yù)測樣本外的第一個(gè)數(shù)字。
而 forecast函數(shù)，是對樣本外的數(shù)據(jù)進(jìn)行預(yù)測。
但是這兩個(gè)函數(shù)的預(yù)測結(jié)果是一樣的。
另外，需要提到的是，ARIMA算法一般只能預(yù)測一點(diǎn)點(diǎn)，越長越不準(zhǔn)確，即便是簡單的正弦函數(shù)也不能準(zhǔn)確預(yù)測。

# 預(yù)測15日close股價(jià)，即10日之后的第三個(gè)交易日的收盤價(jià)
# 但是通過上面的局部數(shù)據(jù)觀察發(fā)現(xiàn)，預(yù)測的數(shù)據(jù)趨勢會(huì)延遲1個(gè)交易日，所以就取f[0][1]
f = ARMAModel.forecast(3)  # 樣本外預(yù)測
print("5月15日close時(shí)的股價(jià)為：%.2f 元" % f[0][1])


#----結(jié)束-----

5月15日close時(shí)的股價(jià)為：33.82 元。

順便安利下這次參與的小競賽「數(shù)據(jù)游戲」，是由個(gè)人發(fā)起、異步社區(qū)（郵電出版社）贊助的活動(dòng)。

總體來說，這個(gè)比賽氛圍適中，沒有太大的心理負(fù)擔(dān)。

所以感覺比較適合打算初入的小伙伴，參與之后，這樣的經(jīng)驗(yàn)對以后的學(xué)習(xí)進(jìn)階和職業(yè)發(fā)展都有好處。

學(xué)習(xí)更多數(shù)據(jù)科學(xué)知識(shí)請關(guān)注微信公眾號(hào)：read_csv

參與數(shù)據(jù)科學(xué)活動(dòng)請加 QQ 群：759677734