Yet another 機械学習で株価を予測する (1)

# 1. はじめに
　これまでもQiitaには株価を機械学習/ディープラーニングで予想する、という記事が投稿されてきました (例えば参考文献の1)。果たしてそれに付け加えることがあるのか、という気もしますが考え方の整理も兼ねて投稿したいと思います。主な方針としては

1. 日経平均(N225)がその日上がるか、下がるかを予測する
2. N225を含む指標の、始値 (Open)、高値 (High)、安値 (Low)、終値 (Close)のみを使う
3. 機械学習のライブラリとしてはscikit-learnを使う

です。1については今日の終値が昨日の終値より安いか高いか正確に予測できたとしても、それで利益を出すのは難しそうだからです。ただし予測の難易度は上がります(70%の精度は夢のまた夢です)。2はデータの入手性によります。3はディープラーニング系のライブラリ(tensorflowなど)は使わないということです。これはハードウェアの制約ではなく、実際にやってみた結果、ディープラーニングを使っても顕著に成績が良くならなかったからです (スキルが向上すればディープラーニングで結果を改善できるかもしれませんが)。

　動作の確認に利用したPythonのバージョンは3.6.3、numpyは1.13.3、scikit-learnは0.19.1、およびpandasの0.20.3を使っています。また例示するコードは個人的に使用される場合はご自由にお使いください。ただし使用に伴って生じたいかなる結果も自己責任でお願いします。

# 2. 今回の目標
　まずは前営業日の日経平均の指標から本日の日経平均が上がるか下がるかを予測します。
ただししきい値(例えば0.1 %)を設けて、それ以上上がるか下がるかを予測します。実際には上昇を予測する分類器と下落を予測する分類器の2台の分類器を1セットとして動作させます。またアルゴリズムとしては一般的に広く用いられているランダムフォレストを使ってみます。

# 3. データの準備
　日経平均の時系列データの入手法はいろいろありますが、現時点はアメリカの[Yahoo! Finance](https://finance.yahoo.com/)から無料でダウンロードできます([Nikkei 225](https://finance.yahoo.com/quote/%5EN225/history?p=%5EN225))。「^N225.csv」という名前でダウンロードできますが、私は「N225.csv」のように改名して使っています。

　ダウンロードしてきたデータを見ると

```php:N225.csv
Date,Open,High,Low,Close,Adj Close,Volume
2001-01-04,13898.089844,13990.570313,13667.679688,13691.490234,13691.490234,0
2001-01-05,13763.219727,13947.059570,13725.459961,13867.610352,13867.610352,0
2001-01-08,null,null,null,null,null,null
2001-01-09,13732.849609,13732.849609,13460.820313,13610.509766,13610.509766,0
2001-01-10,13593.160156,13593.160156,13349.150391,13432.650391,13432.650391,0
2001-01-11,13433.089844,13436.610352,13123.809570,13201.070313,13201.070313,0
```
のようになっています。今回はこのうち"Data"をインデックスとし、"Open"、"High"、"Low"、"Close"列のデータを予測に用います。また2001年1月8日は月曜日ですが、成人の日で日本の市場は休みなので"null"として表示されています。この行は不要なので取り除きます。以上の操作をpandasで行うコードは

```php:dummy01.py
import pandas as pd

df = pd.read_csv("N225.csv", na_values=["null"])
df["Date"] = pd.to_datetime(df["Date"])
df = df.set_index("Date")
df = df[["Open", "High", "Low", "Close"]]
df = df.dropna()
```
のようになります。以上の処理を行うとdfは

```php:dummy02.csv
                    Open          High           Low         Close
Date
2001-01-04  13898.089844  13990.570313  13667.679688  13691.490234
2001-01-05  13763.219727  13947.059570  13725.459961  13867.610352
2001-01-09  13732.849609  13732.849609  13460.820313  13610.509766
2001-01-10  13593.160156  13593.160156  13349.150391  13432.650391
2001-01-11  13433.089844  13436.610352  13123.809570  13201.070313
```
のようになっています。

# 4. データの加工
　次に後の学習で使いやすいようにデータを加工します。予測に使うデータは一日分なので株価の絶対値(例えば2万3456円)は必要ありません。むしろOpenが1000円でCloseが1010円の場合と、Openが10000円でCloseが10100円の場合は同じと考えて良いでしょう。そこでデータの規格化を行います。ここでは簡単に、その日のCloseで割ることにします。この方式ではCloseの値は常に1になります。

```php:dummy03.py
df["Open"] /= df["Close"]
df["High"] /= df["Close"]
df["Low"] /= df["Close"]
df["Close"] = 1.0
```
　また最も興味のある値は一日の間の変化率、つまりClose/Openですので、Result = Close/Open という列を新たに作ります。つまり、Open～Closeのデータを入力してResultが(1＋しきい値)より大きくなるか、(1－しきい値)より小さくなるかを予測できればよいわけです。ただし、その日のResultを予測しても意味がありませんから、1営業日後のResultを予測するようにします。そこでResultの列だけ1営業日前にずらします。

```php:dummy04.py
df["Result"] = 1/df["Open"].shift(-1)
df = df.dropna()
```
　これでdfは

```php:dummy05.csv
                Open      High       Low  Close    Result
Date
2001-01-04  1.015090  1.021844  0.998261    1.0  1.007585
2001-01-05  0.992472  1.005729  0.989749    1.0  0.991091
2001-01-09  1.008989  1.008989  0.989002    1.0  0.988192
2001-01-10  1.011949  1.011949  0.993784    1.0  0.982728
2001-01-11  1.017576  1.017842  0.994147    1.0  1.007666
```
のようになります。株価とResultの対応が1行ずれているのがわかると思います。dfのOpen～Closeまでの部分を入力にして、Resultの値と(1＋しきい値)あるいは(1－しきい値)の大小を学習させれば良いわけです。

# 5. データの分割
　データを学習データと評価データに分割します。ここでは最も単純にsklearn.model_selection.train_test_splitを使います。

```php:dummy06.py
from sklearn.model_selection import train_test_split

columns_input = ["Open", "High", "Low", "Close"]
X_train, X_test, r_train, r_test = train_test_split(df[columns_input], 
                                                    df["Result"], 
                                                    random_state=0,
                                                    test_size=0.3)
```
でdfの70%を学習データに、30%を評価データにします。ただしr_train、r_testはResultの値そのままなので、

```php:dummy07.py
threshold = 0.001
for polarity in ["positive", "negative"]:
    if polarity == "positive":
        y_train = [1 if r - 1.0 >= threshold else 0 for r in r_train]
        y_test = [1 if r - 1.0 >= threshold else 0 for r in r_test]
    else:
        y_train = [1 if 1.0 - r >= threshold else 0 for r in r_train]
        y_test = [1 if 1.0 - r >= threshold else 0 for r in r_test]
```
として上昇(positive)を学習する場合には(1＋しきい値)以上を1に、それ以外を0にします。下落(negative)は逆に(1ーしきい値)以下の時に1、それ以外の場合を0にします。

# 6. 学習
　学習部分のコードは

```php:dummy08.py
from sklearn.ensemble import RandomForestClassifier

clf.fit(X_train, y_train)
```
です。これでclfに学習結果が格納されます。


# 7. 結果出力
 学習データのスコアはclf.score(X_train, y_train)で、評価データのスコアはclf.score(X_test, y_test)で求めることができます。しかしこれではどれだけ利益が出るのかは評価できないので

```php:dummy09.py
if polarity == "positive":
    gain_pred = clf.predict(X_test) * (r_test - 1)
else:
    gain_pred = clf.predict(X_test) * (1 - r_test)
```
として上昇の場合と下落の場合の予想に対する損益を計算します。gain_predは1営業日ごとの損益を格納した行列ですが、最終的には平均を取って(値が小さいので)%表示します。

```php:dummy10.py
print('  {} training accuracy: {:.3f}'.format(polarity, clf.score(X_train, y_train)))
print('  {} test accuracy: {:.3f}'.format(polarity, clf.score(X_test, y_test)))
print('  {} mean gain: {:.3f} %'.format(polarity, gain_pred.mean()*100))
```
ソースコード全体は最後に載せます。

# 8. 実行結果
　今回作成したコードを実行した結果は

```
  positive training accuracy: 0.945
  positive test accuracy: 0.515
  positive mean gain: -0.008 %
  negative training accuracy: 0.942
  negative test accuracy: 0.476
  negative mean gain: -0.012 %
```
です。テストデータに対する予測精度は50%近辺となっています。また残念ながら上昇/下落ともに利益は出ないようです。しかしソース中のrandom_stateを1にして実行してみると

```
  positive training accuracy: 0.947
  positive test accuracy: 0.515
  positive mean gain: -0.013 %
  negative training accuracy: 0.945
  negative test accuracy: 0.519
  negative mean gain: 0.009 %
```
のようになり、下落側は利益が出そうです。このように今回の予測方法は利益が出るか出ないか際どく、実際に利益を出せる方法かどうかはもう少し精密に評価する必要があります。次回はそのあたりについて考えてみたいと思います。

# 参考文献
1. [日経平均株価の予測をTensorFlowで機械学習させ精度の向上を試みた　62％から68％位になった](https://qiita.com/verizi/items/a3642b44f0ef5ad299f5)

# ソース
```php:qiita01.py
#!/usr/bin/env python
# -*- coding: utf-8 -*-
# Copyright by troilus (2018)

import pandas as pd
import numpy as np
from sklearn.model_selection import train_test_split
from sklearn.ensemble import RandomForestClassifier


if __name__ == '__main__':
    table_root = "./tables"

    df = pd.read_csv("N225.csv", na_values=["null"])
    df["Date"] = pd.to_datetime(df["Date"])
    df = df.set_index("Date")
    df = df[["Open", "High", "Low", "Close"]]

    df = df.dropna()


    df["Open"] /= df["Close"]
    df["High"] /= df["Close"]
    df["Low"] /= df["Close"]
    df["Close"] = 1.0
    df["Result"] = 1/df["Open"].shift(-1)
    df = df.dropna()

    columns_input = ["Open", "High", "Low", "Close"]
    X_train, X_test, r_train, r_test = train_test_split(df[columns_input], 
                                                        df["Result"], 
                                                        random_state=0,
                                                        test_size=0.3)
    X_train = X_train.values
    X_test = X_test.values

    clf = RandomForestClassifier(n_estimators=5)
    threshold = 0.001
    for polarity in ["positive", "negative"]:
        if polarity == "positive":
            y_train = [1 if r - 1.0 >= threshold else 0 for r in r_train]
            y_test = [1 if r - 1.0 >= threshold else 0 for r in r_test]
        else:
            y_train = [1 if 1.0 - r >= threshold else 0 for r in r_train]
            y_test = [1 if 1.0 - r >= threshold else 0 for r in r_test]

        clf.fit(X_train, y_train)

        if polarity == "positive":
            gain_pred = clf.predict(X_test) * (r_test - 1)
        else:
            gain_pred = clf.predict(X_test) * (1 - r_test)


        print('  {} training accuracy: {:.3f}'.format(polarity, clf.score(X_train, y_train)))
        print('  {} test accuracy: {:.3f}'.format(polarity, clf.score(X_test, y_test)))
        print('  {} mean gain: {:.3f} %'.format(polarity, gain_pred.mean()*100))
```