【Python】俺の株価予測 【HFT】

# 1. モチベーション
* Qiitaの記事を見ていてもいろいろな株価予測についての記事がありますが、本稿では僕も自分なりに本気で株価の予測をしてみようと思いました。
* ここでの株価予測は数分程度先の**短期の予測**です。
* **取引所の板を見ていて、買いが強いなーとか、売りが多いなーっていう印象を受けることがありますが、それを元に株価を予測しよう**というのが今回のアプローチです。
* 参考にしたのは、["Deep Convolutional Neural Networks for Limit Order Books"](https://arxiv.org/pdf/1808.03668.pdf)という論文です。Limit Order Book (LOB)は板情報という意味です。
* また、データセットは[FI-2010](https://etsin.avointiede.fi/dataset/urn-nbn-fi-csc-kata20170601153214969115)というフィンランドのヘルシンキ証券取引所のデータを使用しています。

# 2. 板情報について
* 板情報については[ザイ・オンライン](https://diamond.jp/articles/-/121153)の記事が分かりやすく解説してくれています。簡単に説明すると、以下のような売りと買いのニーズを集約したものが板になります。取引所には様々な人の売買注文が集まってきます。●●円で○○株買いたい（売りたい）というような注文をその値段ごとに集約し、その株数合計を表示してくれています。

|売数量(ASK)|株価|買数量(BID)|
|--:|---|:--|
|500|670||
|400|669||
|600|668||
||667|300|
||666|1,200|
||665|400|

* 株価がなぜどうやって変動するのかについてはいろいろな考え方がありますが、売りたい人たちの売り数量よりも買いたい人たちの売り数量が多い場合、株価が上がっていくというふうにここでは考えます。
* 上の板の例で、ここに100株買いたい人が現れて新たに注文を出す場合を考えてみます。注文の方法はだいたい２パターンに分かれます。
* ひとつめは、668円に600株の売り注文が出ているのに自分の買い注文をぶつける方法です。この時、板は下のように変化します。

|ASK|株価|BID|
|--:|---|:--|
|500|670||
|400|669||
|**500**|668||
||667|300|
||666|1,200|
||665|400|

* もうひとつは例えば667円に買い注文を出す方法です。必ず約定するとは限りませんが、667円に売りの注文が出てきた時に買いが成立します。このとき、板は下のように変化します。

|ASK|株価|BID|
|--:|---|:--|
|500|670||
|400|669||
|600|668||
||667|**400**|
||666|1,200|
||665|400|

* ここで言いたいことは以下の通りです。
 * 買いたい人が増えた時板情報は、**ASKの数量が減るかBIDの数量が増える**という反応になる。
 * 売りたい人が増えた時板情報は、**ASKの数量が増えるかBIDの数量が減る**という反応になる。
* **板情報を丁寧に追っていけば、その後の株価推移について何かしらの示唆が得られる**ような気がしてきませんか？？

# 3. FI-2010データセット
## FI-2010データセットとは
* [Benchmark Dataset for Mid-Price Forecasting of Limit Order Book Data with Machine Learning Methods](https://arxiv.org/pdf/1705.03233.pdf)という論文に詳しく説明があり、以下そこからの抜粋です。
* フィンランドのヘルシンキ取引所から取った板情報のデータです。
* 2010年6月の1日から14日までの期間のデータです。
* 対象の銘柄は5銘柄です。ケスコ(KESBV)、オウトクンプ(OUT1V)、サンポコンツェルン(SAMPO)、ロータールーキー(RTRKS)、バルチラ(WET1V)になります。僕は見事に全銘柄知らない銘柄でした。
* 板に変化があったごとにデータをサンプルしているようです。（それにしてはデータ数が少ない気がするのでここは僕の理解が違うかも。ただ、2010年のことなので分かりません。）
* ダウンロードできるのは正規化された後のデータです。以下の3種類の正規化が用意されています。
 *  Z-score<br>$\quad x_i^{Zscore} = \frac{x_i - x_{mean}}{x_{std}}$<br>$\quad \rm where \quad x_{mean} = \frac{1}{N} \sum_{j=1}^{N} x_j, \quad x_{std} = \sqrt{\frac{1}{N} \sum_{j=1}^{N} (x_j - x_{mean})^2}$
 * Min-Max Scaling<br> $\quad x_i^{(MM)} = \frac{x_i - x_{min}}{x_{max} - x_{min}}$
 * Decimal Precision<br> $\quad x_i^{DP} = \frac{x_i}{10^k}$<br>$\quad$where k is the integer that will give the maximum value for $|x_i^{(DP)}|<1$
* [ここのリンク](https://etsin.avointiede.fi/dataset/urn-nbn-fi-csc-kata20170601153214969115)のData AvailabilityのAccess this dataset freely.をクリックするとダウンロードできます。

## データの概観
* 以下ではDecimal Precisionで正規化されたデータを見て行きます。これが一番データの数字を直感的に理解しやすいと思うためです。まずは全データ10日分のうち、初日のデータを読み込んできます。

```python
data = pd.read_csv('Train_Dst_Auction_DecPre_CF_1.txt', 
                   header=None, delim_whitespace=True)
print(data.shape)
#=> (149, 47342)
```
* 全部で47,342件のデータがあり、それぞれ149の要素で構成されていることが分かります。
* 全体像が掴みにくいのでヒートマップで見てみます。

```python
plt.figure(figsize=(20,10))
plt.imshow(data, interpolation='nearest', vmin=0, vmax=0.75, 
           cmap='jet', aspect=data.shape[1]/data.shape[0])
plt.colorbar()
plt.grid(False)
plt.show()
```
![save.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/191326/244dd423-7966-76c3-c706-84fc0fb0cfab.png)

* こう見ると、確かに５銘柄のデータが横につながって入っている様子が分かります。上から144個の行は特徴量、最後の5行はラベルを表しています。また、具体的な特徴量については論文に以下のような記述があります。
<img width="737" alt="Screen Shot 2020-01-11 at 21.37.06.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/191326/6b447d07-fe1b-26bf-c25a-1a5eaf1531c7.png">
* 基本的な板情報は$u_1$の40行に入っています。それ以降は板情報を加工して作成した特徴量となっています。板情報についてどのようにデータが入っているか詳細に見てみます。試しに最初の1列目を使います。

```python
lob = data.iloc[:40,0].values
lob_df = pd.DataFrame(lob.reshape(10,4), 
                      columns=['ask','ask_vol','bid','bid_vol'])
print(lob_df)
```
|  | ask | ask_vol | bid | bid_vol |
|--:|--:|--:|--:|--:|
| 0 | 0.2631 | 0.00392 | 0.2616 | 0.00663 |
| 1 | 0.2643 | 0.00028 | 0.2615 | 0.00500 |
| 2 | 0.2663 | 0.00165 | 0.2614 | 0.00500 |
| 3 | 0.2664 | 0.00500 | 0.2613 | 0.00043 |
| 4 | 0.2667 | 0.00039 | 0.2612 | 0.00646 |
| 5 | 0.2710 | 0.00700 | 0.2611 | 0.00200 |
| 6 | 0.2745 | 0.00200 | 0.2609 | 0.00199 |
| 7 | 0.2749 | 0.00487 | 0.2602 | 0.00081 |
| 8 | 0.2750 | 0.00300 | 0.2600 | 0.00197 |
| 9 | 0.2769 | 0.01000 | 0.2581 | 0.01321 |

* ここまでくると理解しやすいです。1番上に最良ASK、BIDがあって、次第に最良気配から遠いところの板情報になっていってます。<br>$\quad u_1 = \\{P_i^{ask},V_i^{ask},P_i^{bid},V_i^{bid}\\}_{i=1}^{10}$<br>という記述の通り、40行のデータは『1番目（最良）ASK価格、ASK数量、BID価格、BID数量、2番目のASK価格、ASK数量、BID価格、BID数量、、、、』という形で格納されています。

# 4. モデル
## 学習データとラベル
* ここからは実際に使う機械学習モデルの説明をしていきます。学習するにあたって学習データ$\mathbb X$と対応するラベル$\mathbb y$が必要になりますが、これを構成する$(\mathbb x_t, y_t)$についてです。
* **学習に用いるデータ**ですが、ある時点$t$の板データを<br>$\quad v_t = \\{P_{t,i}^{ask},V_{t,i}^{ask},P_{t,i}^{bid},V_{t,i}^{bid}\\}_{i=1}^{10}$<br>とし、これを直近の$p$個集めたものを１つの学習データ（$\mathbb x\_t$）とします。具体的には、<br>$\quad \mathbb x\_t = \begin{pmatrix} v\_{t-p+1} \\\v\_{t-p+2} \\\ \vdots \\\ v\_t\end{pmatrix} = \begin{pmatrix} P\_{t-p+1,1}^{ask} & V\_{t-p+1,1}^{ask} & P\_{t-p+1,1}^{bid} & V\_{t-p+1,1}^{bid} & P\_{t-p+1,2}^{ask} & \cdots & P\_{t-p+1,10}^{bid} & V\_{t-p+1,10}^{bid} \\\ P\_{t-p+2,1}^{ask} & V\_{t-p+2,1}^{ask} & P\_{t-p+2,1}^{bid} & V\_{t-p+2,1}^{bid} & P\_{t-p+2,2}^{ask} & \cdots & P\_{t-p+2,10}^{bid} & V\_{t-p+2,10}^{bid} \\\ \vdots & \vdots & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\\ P\_{t,1}^{ask} & V\_{t,1}^{ask} & P\_{t,1}^{bid} & V\_{t,1}^{bid} & P\_{t,2}^{ask} & \cdots & P\_{t,10}^{bid} & V\_{t,10}^{bid}\end{pmatrix}$<br>という$p×40$の行列になります。これをCNNで畳み込みしたのち、LSTMに通すので、最も古いデータが１行目、$t$時点のデータは1番下の行に入っています。
* **ラベル**（$y_t$）は$t$時点以降の$k$期間の仲値の平均が閾値$\alpha$を基準に上昇しているか、下落しているか、横ばいかをもとに割り振ります。まず、仲値（$p_t$）とはそれぞれの時点での最良ASKとBIDの平均を指すので、<br>$\quad p_t = \frac{P_{t,1}^{ask} + P_{t,1}^{bid}}{2}$<br> となります。さらに、期間$k$における仲値の平均値（$m_{+}(t)$）とその騰落率（$l_t$）は、<br>$\quad m_{+}(t) = \frac{1}{k} \sum_{i=1}^{k} p_{t+i}, \quad l_t = \frac{m_{+}(t) - p_t}{p_t}$<br>と表すことができます。最後に閾値（$\alpha$）を基準にして、<br>$\quad y_t =\left\\{ \begin{array}{}1, & l_t>\alpha \\\ -1, & l_t<-\alpha \\\ 0, & \rm otherwise\end{array}\right.$<br>という形にラベリングします。

## モデルアーキテクチャ
* 先にモデルの例をあげておくと以下の様になります。まず各時点での板情報をCNNで畳み込み、最後にLSTMで時系列の関係を処理するという流れになります。
<img width="558" alt="Screen Shot 2020-01-12 at 15.03.41.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/191326/4265d5cb-8fd6-3793-ac8e-5e8340d39375.png">
 * **input_1**レイヤから順に見ていきます。Output Shapeは（バッチサイズ、参照する板情報数（$p$）、ひとつの板情報に含まれるデータ数（40）、1）という意味です。ここでは$p=50$としています。
 * **conv2d_1**レイヤは最初の畳み込み層です。カーネルのサイズは$1×2$、ストライドも$1×2$としています。$\mathbb x_t$には板情報の価格と数量のペアが羅列されていますが、ここではこのひとつひとつのペアに対して畳み込みを行なっています。カーネル数は8を指定しています。ここでは行方向の畳み込みをしているのみなので、Output Shapeの参照している板情報数は$p(=50)$のままで変わらず、ひとつの板情報に含まれるデータ数が半減しています。
 * **conv2d_2**レイヤもカーネルのサイズが$1×2$、ストライドが$1×2$の畳み込み層です。カーネル数は8を指定しています。ここの畳み込みは板情報の「ASKサイドの価格・数量をconv2d_1で畳み込んだもの」と「BIDサイドの価格・数量をconv2d_1で畳み込んだもの」に対して畳み込みをしており、最良気配なら最良気配についてASK、BIDの価格・数量という４つの数字をひとまとめにする役割を果たしています。conv2d_1とconv2d_2のウエイト次第では対応するASK、BIDの加重平均仲値を計算するような振る舞いをすることを考えると少しイメージが掴めます。なお、この加重平均仲値はマイクロプライス（$p_t^{(micro)}$）とも呼ばれ、以下の様に定義されます。<br> $\quad p_{t,i}^{(micro)} = \frac{P_\{t,i}^{ask} V_\{t,i}^{ask} + P_\{t,i}^{bid} V_\{t,i}^{bid}}{V_\{t,i}^{ask} + V_\{t,i}^{bid}}$
 * **conv2d_3**レイヤはカーネルのサイズが$1×10$、ストライドが$1$の畳み込み層です。カーネル数は8を指定しています。そもそもASK、BID両サイドとも10本の板を参照しているため、conv2d_2のアウトプットはひとつの板情報ごとに10個の数字からなっています。この10の数字をひとまとめに畳み込むのがこのレイヤです。これである時点での板情報がひとつの数字に畳み込まれました。
 * **reshape_1**レイヤはconv2d_3のアウトプットを次のlstm_1にフィードしてやるためにあります。
 * **ltsm_1**レイヤではここまで畳み込まれてきた板情報についてその時系列での関係をキャプチャしようとしています。ユニット数は8を指定しています。
 * **dense_1**レイヤはLSTMのアウトプットを受けるシンプルな隠れ層です。
 * **dense_2**レイヤはこのネットワークの出力層となります。ラベルの種類数に合わせて出力は3つ、activationにはsoftmaxを使っています。

# 5. 実装
## データの前処理
* ここでは先ほど読み込んだデータのうち、最もサンプル数が多い5つ目の銘柄のデータを用いてモデルを動かしていきます。

```python
# 板情報は最初の40行に入っている。5つめの銘柄のデータとして29738~47294を指定。
lob = data.iloc[:40, 29738:47294].T.values
# ここで価格・数量ごとに標準化する。
lob = lob.reshape(-1,2)
lob = (lob - lob.mean(axis=0)) / lob.std(axis=0)
lob = lob.reshape(-1,40)
lob_df = pd.DataFrame(lob)
# 標準化されていない仲値を計算する。
lob_df['mid'] = (data.iloc[0,29738:47294].T.values + data.iloc[2,29738:47294].T.values) / 2
```
* 仲値をプロットしてみると以下のようになりました。まあよくある株価チャートです。現値に近い全ての呼値に指値が出ていないこともあり細かい振動が多くなっている印象があります。株価が一方向に動き続けているデータではない点はモデルを作成する上でよさそうです。
![save.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/191326/59ea0893-2d44-a9ac-ecd1-031774cb2219.png)
* 次にラベルを作成していきます。

```python
# パラメータを指定する。
p = 50
k = 50
alpha = 0.0003
# パラメータをもとに仲値からラベルを作成する。
lob_df['lt'] = (lob_df['mid'].rolling(window=k).mean().shift(-k)-lob_df['mid'])/lob_df['mid']
lob_df = lob_df.dropna()
lob_df['label'] = 0
lob_df.loc[lob_df['lt']>alpha, 'label'] = 1
lob_df.loc[lob_df['lt']<-alpha, 'label'] = -1
```

* 今回のパラメータによるラベリングは以下のようなイメージになります。それなりに株価が上昇する際には上昇シグナルを、下落するときには下落シグナルを発しています。
![save.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/191326/405fc8ac-bd32-bce0-798f-cf61c473f1cb.png)


* 必要なライブラリを読み込みます。

```python
from sklearn.model_selection import train_test_split
from keras.utils import to_categorical
from keras.layers import Conv2D, Dense, Reshape, Input, LSTM
from keras import Model, backend
import tensorflow as tf
```

* 学習データを作成していきます。

```python
# 学習データを作成します。
X = np.zeros((len(lob_df)-p+1, p, 40, 1))
lob = lob_df.iloc[:,:40].values
for i in range(len(lob_df)-p+1):
    X[i] = lob[i:i+p,:].reshape(p,-1,1)
y = to_categorical(lob_df['label'].iloc[p-1:], 3)
print(X.shape, y.shape)
#=> (17457, 50, 40, 1) (17457, 3)
```

* 最後に学習データとテストデータに分割します。

```python
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)
```

* これで前処理は完了です！

## モデル構築
* ここからはkerasを使ってニューラルネットワークのモデルを作成していきます。SequentialよりもFunctional APIに慣れているのでこっちで記述します。

```python
tf.reset_default_graph()
backend.clear_session()

inputs = Input(shape=(p,40,1))
x = Conv2D(8, kernel_size=(1,2), strides=(1,2), activation='relu')(inputs)
x = Conv2D(8, kernel_size=(1,2), strides=(1,2), activation='relu')(x)
x = Conv2D(8, kernel_size=(1,10), strides=1, activation='relu')(x)
x = Reshape((p, 8))(x)
x = LSTM(8, activation='relu')(x)
x = Dense(16, activation='relu')(x)
outputs = Dense(3, activation='softmax')(x)

model = Model(inputs=inputs, outputs=outputs)
model.compile(optimizer='adam', loss='categorical_crossentropy', metrics=['accuracy'])
```

## いざ学習！
* 後は学習させるのみです。

```python
epochs = 50
batch_size = 256
history = model.fit(X_train, y_train,
                    epochs=epochs,
                    batch_size=batch_size,
                    verbose=1,
                    validation_data=(X_test, y_test))
```
> Epoch 100/100
13965/13965 [==============================] - 5s 326us/step - loss: 0.6526 - acc: 0.6808 - val_loss: 0.6984 - val_acc: 0.6595

* エポックごとのlossとaccuracyは以下のようになりました。そこそこうまく学習が進んでいる様子が分かります。

![save.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/191326/e20badf1-058a-396a-39e2-b306479e00ad.png)
![save.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/191326/86f6c1da-dd80-2e11-46be-01506eceb76b.png)

# 6. 考察
* 学習データに対してはAccuracyが0.6808、**テストデータに対しては0.6595**という思いのほか良好な結果が得られました。
* テストデータの結果をヒートマップで可視化してみます。正しく分類できていることもいいことですが、実用化を考えると「株価が上昇するときに下落シグナルを出してしまうこと」、「下落するときに上昇シグナルを出してしまうこと」が損失につながるので最も問題となります。この点に関して、まったく逆のシグナルを出してしまっているケースはかなり少なく、ポジティブな結果となりました。
![save.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/191326/e5beab9a-306c-1c30-7368-c293ce56fec4.png)

* 今後の課題
 * ハイパーパラメータチューニング。
 * 他の銘柄のデータも合わせて学習をする。 
 * 他の取引日のデータも利用して学習をする。（一般に、ある１日のデータをもとにモデルを作ってもそれを別日に使うとパフォーマンスが落ちることが分かっている。）
 * 最近の日本の株式市場でも同様のパフォーマンスを出せるか検証する。