版權說明：本公号轉載文章旨在學習交流，不用于任何商業用途，版權歸原作者所有，如有異議，敬請背景聯絡我們，議定合作或删除，我們将第一時間按版權法規定妥善處理，非常感謝！
作為一種技術手段，預測在金融、證券領域的應用非常廣泛，尤其是對股票價格的預測。我們介紹一下獲得股票資料的方法，并基于此對資料進行預處理，接着使用資料分析方法，建立基礎特征，進一步建構預測模型，且基于新資料驗證模型效果。拟使用VAR及LSTM兩種算法建立預測模型。
擷取股票資料
股票資料通常可從新浪股票、雅虎股票等網頁上擷取，此外還有一些炒股軟體，如同花順、通達信等都提供了非常清楚的股票資料展示和圖表呈現。如果要獲得實時的股票資料，可以考慮使用新浪股票提供的接口擷取資料。以大秦鐵路(股票代碼：601006)為例，如果要擷取它的最新行情，隻需通路新浪的股票資料接口(具體可以百度)，該接口會傳回一串文本，例如：
1 var hq_str_sh601006="大秦鐵路,6.980,6.960,7.010,7.070,6.950,7.010,7.020,121033256,847861533.000,18900, 7.010,214867,7.000,66500,6.990,386166,6.980,336728,6.970,273750,7.020,836066,7.030,630800,7.040,936306,7.050,579400,7.060,2016-03-18,15:00:00,00";
這個字元串由許多資料拼接在一起，不同含義的資料用逗号隔開了，按照程式員的思路，順序号從0開始。
0：，股票名字1：<< span="">6.980>，今日開盤價2：<< span="">6.960>，昨日收盤價3：<< span="">7.010>，目前價格4：<< span="">7.070>，今日最高價5：<< span="">6.950>，今日最低價6：<< span="">7.010>，競買價，即“買一”報價7：<< span="">7.020>，競賣價，即“賣一”報價8：<< span="">121033256>，成交的股票數，由于股票交易以一百股為基本機關，是以在使用時，通常把該值除以一百9：<< span="">847861533.000>，成交金額，機關為“元”，為了一目了然，通常以“萬元”為成交金額的機關，是以通常把該值除以一萬10：<< span="">18900>，“買一”申請4695股，即47手11：<< span="">7.010>，“買一”報價12：<< span="">214867>，“買二”13：<< span="">7.000>，“買二”14：<< span="">66500>，“買三”15：<< span="">6.990>，“買三”16：<< span="">386166>，“買四”17：<< span="">6.980>，“買四”18：<< span="">336728>，“買五”19：<< span="">6.970>，“買五”20：<< span="">273750>，“賣一”申報3100股，即31手21：<< span="">7.020>，“賣一”報價(22,23),(24,25),(26,27),(28,29)分别為“賣二”至“賣四的情況”30：<< span="">2016-03-18>，日期31：<< span="">15:00:00>，時間
這個接口對于JavaScript程式非常友善，如果要檢視該股票的日K線圖，可通路新浪股票的K線圖接口(具體可百度)，便可得到日K線圖。
日K線圖 如果要檢視該股票的分時線，可通路連結新浪股票的分時線圖接口(具體可百度)，便可得到分時線圖。
分時線圖 對于周K線和月K線的查詢，可分别通路新浪股票的周K線圖和月K線圖的接口(具體可百度)。Python中我們可以使用pandas_datareader庫來擷取股票資料，預設是通路yahoofinance的資料，其中包括上證和深證的股票資料，還有港股資料，該庫隻能擷取股票的曆史交易記錄資訊：如最高價、最低價、開盤價、收盤價以及成交量，無法擷取個股的分筆交易明細曆史記錄。上證代碼是ss，深證代碼是sz，港股代碼是hk，比如茅台：6000519.ss，萬科000002.sz，長江實業0001.hk。這裡以貴州茅台股票為例，說明pandas_datareader庫中股票資料的擷取方法及簡單的可視化，代碼如下：
1import pandas as pd 2import pandas_datareader.data as web 3import datetime as dt 4data = web.DataReader('600519.ss','yahoo', dt.datetime(2019,8,1),dt.datetime(2019,8,31)) 5data.head() 6 High Low Open Close Volume Adj Close 7# Date 8# 2019-08-01 977.000000 953.020020 976.51001 959.299988 3508952 959.299988 9# 2019-08-02 957.979980 943.000000 944.00000 954.450012 3971940 954.45001210# 2019-08-05 954.000000 940.000000 945.00000 942.429993 3677431 942.42999311# 2019-08-06 948.000000 923.799988 931.00000 946.299988 4399116 946.29998812# 2019-08-07 955.530029 945.000000 949.50000 945.000000 2686998 945.0000001314kldata=data.values[:,[2,3,1,0]] # 分别對應開盤價、收盤價、最低價和最高價15from pyecharts import options as opts16from pyecharts.charts import Kline1718kobj = Kline().add_xaxis(data.index.strftime("%Y-%m-%d").tolist()).add_yaxis("貴州茅台-日K線圖",kldata.tolist()).set_global_opts(19 yaxis_opts=opts.AxisOpts(is_scale=True),20 xaxis_opts=opts.AxisOpts(is_scale=True),21 title_opts=opts.TitleOpts(title=""))22kobj.render()
貴州茅台股票日K線圖如圖。
為給定時間序列的财務圖表，代碼中對象data包含6個屬性，依次為Open(開盤價)、High(最高價)、Low(最低價)、Close(收盤價)、Volume(成交量)、Adjusted(複權收盤價)。基于收盤價的重要性，可從收盤價的曆史資料中分割訓練集、驗證集、測試集，使用适當的特征，建立預測模型，并實施預測。
基于VAR算法的預測
向量自回歸(VAR)模型就是非結構化的多方程模型，它的核心思想不考慮經濟理論，而直接考慮經濟變量時間時序之間的關系，避開了結構模組化方法中需要對系統中每個内生變量關于所有内生變量滞後值函數模組化的問題，通常用來預測相關時間序列系統和研究随機擾動項對變量系統的動态影響。VAR模型類似聯立方程，将多個變量包含在一個統一的模型中，共同利用多個變量資訊，比起僅使用單一時間序列的ARIMA等模型，其涵蓋的資訊更加豐富，能更好地模拟現實經濟體，因而用于預測時能夠提供更加貼近現實的預測值。此處拟基于貴州茅台股票資料，建立VAR的預測模型。使用後30天的資料作為驗證集，剩餘的資料用于建立預測模型。本節從VAR模型的平穩性檢驗出發，依次完成VAR模型的定階及模組化預測，最終通過分析驗證集上的準确率來評估預測效果。 ▊ 平穩性檢驗 隻有平穩的時間序列才能夠直接建立VAR模型，是以在建立VAR模型之前，首先要對變量進行平穩性檢驗。通常可利用序列的自相關分析圖來判斷時間序列的平穩性，如果序列的自相關系數随着滞後階數的增加很快趨于0，即落入随機區間，則序列是平穩的；反之，序列是不平穩的。另外，也可以對序列進行ADF檢驗來判斷平穩性。對于不平穩的序列，需要進行差分運算，直到差分後的序列平穩後，才能建立VAR模型。此處首先提取用于建立預測模型的基礎資料，并對其進行機關根檢驗，對應的Python代碼如下：
1import statsmodels.tsa.stattools as stat 2import pandas_datareader.data as web 3import datetime as dt 4import pandas as pd 5import numpy as np 6 7data = web.DataReader('600519.ss','yahoo', dt.datetime(2014,1,1),dt.datetime(2019,9,30)) 8subdata = data.iloc[:-30,:4] 9for i in range(4):10 pvalue = stat.adfuller(subdata.values[:,i], 1)[1]11 print("名額 ",data.columns[i]," 機關根檢驗的p值為：",pvalue)12# 名額 High 機關根檢驗的p值為：0.995520228085040113# 名額 Low 機關根檢驗的p值為：0.994250943975568914# 名額 Open 機關根檢驗的p值為：0.993854819399032315# 名額 Close 機關根檢驗的p值為：0.9950049124079876
可以看到，p值都大于0.01，是以都是不平穩序列。現對subdata進行1階差分運算，并再次進行機關根檢驗，對應的Python代碼如下：
1subdata_diff1 = subdata.iloc[1:,:].values - subdata.iloc[:-1,:].values2for i in range(4):3 pvalue = stat.adfuller(subdata_diff1[:,i], 1)[1]4 print("名額 ",data.columns[i]," 機關根檢驗的p值為：",pvalue)5# 名額 High 機關根檢驗的p值為：0.06# 名額 Low 機關根檢驗的p值為：0.07# 名額 Open 機關根檢驗的p值為：0.08# 名額 Close 機關根檢驗的p值為：0.0
如結果所示，對這4個名額的1階差分單獨進行機關根檢驗，其p值都不超過0.01，是以可以認為是平穩的。 ▊ VAR模型定階 接下來就是為VAR模型定階，可以讓階數從1逐漸增加，當AIC值盡量小時，可以确定最大滞後期。我們使用最小二乘法，求解每個方程的系數，并通過逐漸增加階數，為模型定階，Python代碼如下：
1# 模型階數從1開始逐一增加 2rows, cols = subdata_diff1.shape 3aicList = [] 4lmList = [] 5 6for p in range(1,11): 7 baseData = None 8 for i in range(p,rows): 9 tmp_list = list(subdata_diff1[i,:]) + list(subdata_diff1[i-p:i].flatten())10 if baseData is None:11 baseData = [tmp_list]12 else:13 baseData = np.r_[baseData, [tmp_list]]14 X = np.c_[[1]*baseData.shape[0],baseData[:,cols:]]15 Y = baseData[:,0:cols]16 coefMatrix = np.matmul(np.matmul(np.linalg.inv(np.matmul(X.T,X)),X.T),Y)17 aic = np.log(np.linalg.det(np.cov(Y - np.matmul(X,coefMatrix),rowvar=False))) + 2*(coefMatrix.shape[0]-1)**2*p/baseData.shape[0]18 aicList.append(aic)19 lmList.append(coefMatrix)2021#對比檢視階數和AIC22pd.DataFrame({"P":range(1,11),"AIC":aicList})23# P AIC24# 0 1 13.58015625# 1 2 13.31222526# 2 3 13.54363327# 3 4 14.26608728# 4 5 15.51243729# 5 6 17.53904730# 6 7 20.45733731# 7 8 24.38545932# 8 9 29.43809133# 9 10 35.785909
如上述代碼所示，當p=2時，AIC值最小為13.312225。是以VAR模型定階為2，并可從對象lmList[1]中擷取各名額對應的線性模型。 ▊ 預測及效果驗證 基于lmList[1]中擷取各名額對應的線性模型，對未來30期的資料進行預測，并與驗證資料集進行比較分析，Python代碼如下：
1p = np.argmin(aicList)+1 2n = rows 3preddf = None 4for i in range(30): 5 predData = list(subdata_diff1[n+i-p:n+i].flatten()) 6 predVals = np.matmul([1]+predData,lmList[p-1]) 7 # 使用逆差分運算，還原預測值 8 predVals=data.iloc[n+i,:].values[:4]+predVals 9 if preddf is None:10 preddf = [predVals]11 else:12 preddf = np.r_[preddf, [predVals]]13 # 為subdata_diff1增加一條新記錄14 subdata_diff1 = np.r_[subdata_diff1, [data.iloc[n+i+1,:].values[:4] - data.iloc[n+i,:].values[:4]]]1516#分析預測殘差情況17(np.abs(preddf - data.iloc[-30:data.shape[0],:4])/data.iloc[-30:data.shape[0],:4]).describe()18# High Low Open Close19# count 30.000000 30.000000 30.000000 30.00000020# mean 0.010060 0.009380 0.005661 0.01373921# std 0.008562 0.009968 0.006515 0.01367422# min 0.001458 0.000115 0.000114 0.00013023# 25% 0.004146 0.001950 0.001653 0.00278524# 50% 0.007166 0.007118 0.002913 0.01041425# 75% 0.014652 0.012999 0.006933 0.02230526# max 0.039191 0.045802 0.024576 0.052800
從上述代碼第17行可以看出這4個名額的最大百分誤差率分别為3.9191%、4.5802%、2.4576%、5.28%，最小百分誤差率分别為0.1458%、0.0115%、0.0114%、0.013%，進一步，繪制二維圖表觀察預測資料與真實資料的逼近情況，Python代碼如下：
1import matplotlib.pyplot as plt 2plt.figure(figsize=(10,7)) 3for i in range(4): 4 plt.subplot(2,2,i+1) 5 plt.plot(range(30),data.iloc[-30:data.shape[0],i].values,'o-',c='black') 6 plt.plot(range(30),preddf[:,i],'o--',c='gray') 7 plt.ylim(1000,1200) 8 plt.ylabel("$"+data.columns[i]+"$") 9plt.show()10v = 100*(1 - np.sum(np.abs(preddf - data.iloc[-30:data.shape[0],:4]).values)/np.sum(data.iloc[-30:data.shape[0],:4].values))11print("Evaluation on test data: accuracy = %0.2f%% \n" % v)12# Evaluation on test data: accuracy = 99.03%
該預測效果如下圖，其中黑色實線為真實資料，灰色虛線為預測資料，使用VAR模型進行預測的效果總體還是不錯的，平均準确率為99.03%。針對多元時間序列的情況，VAR模型不僅考慮了其他名額的滞後影響，計算效率還比較高，從以上代碼可以看到，對于模型的拟合，直接使用的最小二乘法，這增加了該模型的适應性。
預測效果
基于LSTM算法的預測
本節主要基于LSTM算法對貴州茅台股票資料進行預測，該算法非常擅長序列資料的模組化，由于引入了遺忘門等更為複雜的内部處理單元來處理上下文資訊的存儲與更新，這樣既可以消除梯度問題的困擾，也可以對存在短期或長期依賴的資料模組化，該算法在文本、語音等序列資料模型中廣泛使用。本節從LSTM模組化的資料要求及網絡結構設計講起，通過設定合理的參數，通過訓練得到模型，并基于該模型進行預測，最後将結果與真實資料進行比較，評估預測效果。 ▊ 資料要求 本節使用LSTM算法對貴州茅台股票資料進行預測，可基于前N條樣本對目前樣本進行預測，是以該模型不需要像DNN那樣，将曆史資料進行複雜轉換，将基礎資料稍加處理就能用于訓練模型。對基礎資料的處理即為對該資料進行重新封裝，将樣本前N期的集合與目前樣本對應上，分别得到訓練資料的輸入與輸出。
所示資料對應關系(具體資料為示意) ▊ 資料預處理 首先，需要将基礎資料重構為包含曆史3周特征資料的基礎資料，以預測日的High(最高價)、Low(最低價)、Open(開盤價)、Close(收盤價)4個名額作為輸出資料。這裡我們使用2014年1月1日至2019年8月31日的貴州茅台股票資料作為訓練資料，使用2019年整個9月的資料作為測試資料，來驗證模型效果。用Python将對全體資料進行标準化，并将基礎資料的特征進行重構，代碼如下：
1SEQLEN = 21 2dim_in = 4 3dim_out = 4 4pred_len = 30 5vmean = data.iloc[:,:4].apply(lambda x:np.mean(x)) 6vstd = data.iloc[:,:4].apply(lambda x:np.std(x)) 7t0 = data.iloc[:,:4].apply(lambda x:(x-np.mean(x))/np.std(x)).values 8X_train = np.zeros((t0.shape[0]-SEQLEN-pred_len, SEQLEN, dim_in)) 9Y_train = np.zeros((t0.shape[0]-SEQLEN-pred_len, dim_out),)10X_test = np.zeros((pred_len, SEQLEN, dim_in))11Y_test = np.zeros((pred_len, dim_out),)12for i in range(SEQLEN, t0.shape[0]-pred_len):13 Y_train[i-SEQLEN] = t0[i]14 X_train[i-SEQLEN] = t0[(i-SEQLEN):i]15for i in range(t0.shape[0]-pred_len,t0.shape[0]):16 Y_test[i-t0.shape[0]+pred_len] = t0[i]17 X_test[i-t0.shape[0]+pred_len] = t0[(i-SEQLEN):i]
如上述代碼所示，SEQLEN表示使用前期資料的長度，dim_in表示輸入資料的次元，dim_out表示輸出資料的次元，pred_len表示預測資料的長度。第5~7行代碼對資料進行zscore标準化，将資料映射到标準正态分布。第12~17行代碼對基礎資料進行重構，分别得到訓練資料X_train、Y_train以及測試資料X_test、Y_test。 ▊ 網絡結構設計 經嘗試，我們使用近3周的曆史資料來訓練LSTM模型，同時，設定隐含層神經元的數量為64。是以，我們可以将LSTM神經網絡按下面的結構進行設計(圖中N可取21，即3周對應的天數)。
LSTM神經網絡結構 ▊ 建立模型 現基于Keras搭建LSTM神經網絡，并基于訓練集對模型進行訓練，Python代碼如下：
1from keras.layers import LSTM, Dense 2from keras.models import Sequential 3model = Sequential() 4model.add(LSTM(64, input_shape=(SEQLEN, dim_in),activation='relu',recurrent_dropout=0.01)) 5model.add(Dense(dim_out,activation='linear')) 6model.compile(loss = 'mean_squared_error', optimizer = 'rmsprop') 7history = model.fit(X_train, Y_train, epochs=200, batch_size=10, validation_split=0) 8# Epoch 1/200 9# 1350/1350 [==============================] - 1s 1ms/step - loss: 0.044710# Epoch 2/20011# 1350/1350 [==============================] - 1s 737us/step - loss: 0.005912# Epoch 3/20013# 1350/1350 [==============================] - 1s 743us/step - loss: 0.004314# ......15# Epoch 200/20016# 1350/1350 [==============================] - 1s 821us/step - loss: 9.2794e-04
如上述代碼所示，我們使用rmsprop算法來優化模型。由于目前的模組化場景是數值預測，是以使用MSE(均方誤差)來定義損失函數。算法經過200次疊代，loss從0.0447降到了9.2794e-04。我們可以基于得到的模型進行進一步預測。 ▊ 預測實作 基于上文得到的模型，進一步編寫Python代碼，對X_test對應的輸出資料進行預測。需要注意的是，直接得到的預測結果是處于标準化的資料空間中的，需要将其還原成原始資料空間的值，結果才有意義。對應的Python代碼如下：
1preddf=model.predict(X_test)*vstd.values+vmean.values
如上述代碼所示，将模型的預測結果pred_y乘以vstd再加上vmean，即可對資料進行還原。preddf即是最終得到的預測資料，可列印其值，代碼如下：
1preddf 2# array([[1069.35781887, 1038.57915742, 1056.77147186, 1053.83827734], 3# [1070.65142282, 1039.58533719, 1057.34561875, 1054.85567074], 4# [1083.58529328, 1052.70457308, 1070.78824637, 1067.49741882], 5# 6# [1186.19297789, 1161.52758381, 1172.33666591, 1170.44623263], 7# [1181.42680223, 1155.14778501, 1166.5726204 , 1165.00336968], 8# [1186.75600881, 1160.84733425, 1172.37636963, 1170.09819923]]) 910preddf.shape11# (30, 4)
如上述代碼所示，preddf是一個的二維資料，包含了2019年9月整月的預測結果。 ▊ 效果評估 對貴州茅台股票資料預測的效果評估可以采用兩種方法。一種方法是對預測的結果與真實結果進行繪圖比較，通過直覺觀察可以知道預測效果，如果預測曲線與真實曲線完全重合或相當接近，則說明預測效果較好；反之，則說明預測模型還需要改進。另一種方法是基于貴州茅台股票資料預測的誤差累計值來計算一個誤差率，進而得到平均精度水準，該值越大說明整體預測效果也就越好，該值越小說明預測模型還存在優化空間。編寫Python代碼，同時實作預測結果與真實資料的對比圖，以及計算累計誤差，進而全面地評估預測效果，代碼如下：
1import matplotlib.pyplot as plt 2plt.figure(figsize=(10,7)) 3for i in range(4): 4 plt.subplot(2,2,i+1) 5 plt.plot(range(30),data.iloc[-30:data.shape[0],i].values,'o-',c='black') 6 plt.plot(range(30),preddf[:,i],'o--',c='gray') 7 plt.ylim(1000,1200) 8 plt.ylabel("$"+data.columns[i]+"$") 9plt.show()10v = 100*(1 - np.sum(np.abs(preddf - data.iloc[-30:data.shape[0],:4]).values)/np.sum (data.iloc[-30:data.shape[0],: 4].values))11print("Evaluation on test data: accuracy = %0.2f%% \n" % v)12# Evaluation on test data: accuracy = 99.01%
預測評估對比圖如下。
我們可以看到，黑色實線為真實資料，灰色虛線為預測資料，橫坐标為日期下标，縱坐标為對應的股票價格。使用LSTM模型進行預測的效果總體還是不錯的，平均準确率為99.01%。對于多元時間序列資料，可嘗試使用LSTM模型，該模型能夠記憶曆史較長的重要資訊，可有效識别曆史資料中存在的規律和模式，如今廣泛應用于包含大量序列資料的場景中。