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

# 1. 今回の目的
1. [Yet another 機械学習で株価を予測する (1)](https://qiita.com/troilus/items/42cf69cf84273f24bc37)
2. [Yet another 機械学習で株価を予測する (2)](https://qiita.com/troilus/items/6ec99c1f8c9c171733c8)

　1では日経平均の日足データから、翌営業日の日経平均が始値から終値にかけて上昇するか下落するか予想するプログラムを作りました。また2ではそのプログラムの検証を簡単な交差検定で行いました。いずれも予想される利益は0近辺で、利益を出すのは難しそうでした。
　そこで今回は入力データの種類を増やし、予測の精度の改善を目指します。追加する指標としてはニューヨーク証券取引所のNYSE総合 (NYA)を使いましょう。また複数のテーブルを読み込みますので、共通な操作は関数としてまとめることにします。

# 2. 複数の入力データを使う場合の注意点
　今回は複数の入力データを使います。それに伴い、入力データが1つだったときには気にする必要のなかった点に気を使う必要が出てきます。以下でそれについて簡単にまとめます。

## a. 列名
　N225のデータだけを使っていたときには学習データの列名をそのまま __Open__、__High__、__Low__、__Close__としていましたが、複数の指標のデータを一つのテーブルにまとめるときにそのままだとわけがわからなくなってしまいます。そこで以降は __N225_Open__、__N225_High__、__N225_Low__、__N225_Close__のような形式にします。

## b. 前後関係
　東京市場とニューヨーク市場の順序関係を考えると例えば

__8月1日の東京市場__ ->  __8月1日のニューヨーク市場__ -> __8月2日の東京市場__

のようになります。都合の良いことにオーバーラップはありません。従って __8月2日の東京市場__ の予測をする場合は __8月1日の東京市場__ と __8月1日のニューヨーク市場__ のデータを使えば良いことになります。一方、__8月2日のニューヨーク市場__ を予測する場合には

__8月1日のニューヨーク市場__ -> __8月2日の東京市場__ -> __8月2日のニューヨーク市場__

という順序関係ですので、入力データは __8月1日のニューヨーク市場__ と __8月2日の東京市場__ ということになり、日付をうまく合わせるような処理が必要になります。今回は東京市場の予測をしたいので、この処理は必要ありません。


## c. 休日の処理
　もう一つ問題なるのは、東京とニューヨークの(日本とアメリカの)休日の違いです。2018年を例に取ると一番劇的なのが12月で

```
                 日本            アメリカ
12/23 (日)    天皇誕生日  
------------------------------------------------
12/24 (月)     振替休日             平日
            (東京市場休場)
------------------------------------------------
12/25 (火)       平日            クリスマス
                            (ニューヨーク市場休場)
-------------------------------------------------
12/26 (水)       平日              平日

```
という順序関係になっています。このように面倒な場合は予測しない、というのも一つの手ではあるのですが、せっかくのデータなので使うことにしましょう。今回は直近の営業日のデータを使って予測することにします。予測したいデータ <= 入力データの形式で表すと

__12月25日の東京市場__ <= __12月21日の東京市場__ ＋ __12月24日のニューヨーク市場__
__12月26日の東京市場__ <= __12月25日の東京市場__ ＋ __12月24日のニューヨーク市場__

とします。このようなデータを作るため、前回までは結果の方を1営業日前にずらしていましたが、今回は結果を固定して入力データの方を後ろにずらします。具体的には

1. N225のデータに対して __Ratio__ ≡ __N225_Close__/__N225_Open__を計算する
2. 各指標のデータのfrequencyを毎日にする。このとき、nanは直近のデータで埋める (ffill)
3. 各指標のデータの__Open__、__High__、__Low__、__Close__を1日後ろにずらす
4. __Ratio__ と各指標のデータを結合する。このとき index は __Ratio__ の日付とする

以上の処理でそれぞれの指標の直前のデータと、評価したい結果が同じ日付に紐づけされるはずです。
　ちょっとややこしい処理が必要になりますが、今回これ以上面倒なことはありません。

### 訂正
最初に投稿したときには休日処理が
1. N225のデータに対して __Ratio__ ≡ __N225_Close__/__N225_Open__を計算する
2. 各指標のデータの__Open__、__High__、__Low__、__Close__を1営業日前にずらす (コード上は後ろでした)
3. ずらした各指標のデータのfrequencyを毎日にする。このとき、nanは直近のデータで埋める (ffill)
4. __Ratio__ と各指標のデータを結合する。このとき index は __Ratio__ の日付とする
となっていましたが、これだと休日を挟んだ場合必要以上に過去のデータを使ってしまいます。


# 3. 今回の変更点
　今回の変更点をまとめます。

## a. 関数化
　データテーブルが複数になったので、読み込みと処理部分を関数化します。読み込み部分は

```php:dummy01.py
def read_table(index='N225'):
    csv = index + '.csv'
    df = pd.read_csv(csv, na_values=["null"])
    df["Date"] = pd.to_datetime(df["Date"])
    df = df.set_index("Date")
    df = df[["Open", "High", "Low", "Close"]]
    return df.dropna()
```
　またデータを正規化し、1営業日後ろにずらし、周期を1日にして列名を変更する関数は

```php:dummy02.py
def normalize_table(table, index='N225'):
    suffix_list = ["Open", "High", "Low", "Close"]
    table["Open"] /= table["Close"]
    table["High"] /= table["Close"]
    table["Low"] /= table["Close"]
    table["Close"] = 1.0
    table = table.resample("D").ffill()
    table = table[suffix_list].shift(1)
    table.columns = [index+"_"+suffix for suffix in suffix_list]
    return table
```
です。

## 訂正
　ここは最初

```php:dummy02a.py
    table = table[suffix_list].shift(1)
    table = table.resample("D").ffill()
```
のようになっていましたが、修正しました。
(ソースも修正してあります)。


## b. データの処理
　以上の関数を使ってデータ作成を変更します。

```php:dummy03.py

indices = ['N225', 'NYA']

df_dict = {}
for index in indices:
    df_dict[index] = read_table(index)
df_ratio = df_dict[indices[0]]['Close']/df_dict[indices[0]]['Open']
df_ratio = df_ratio.to_frame("Ratio")
    
for index in indices:
    df_dict[index] = normalize_table(df_dict[index], index=index)
df = pd.concat([df_dict[index] for index in indices] + [df_ratio], axis=1, join='inner')
df = df[start:end]

columns_input = [index + suffix for index in indices\
                  for suffix in ["_Open", "_High", "_Low", "_Close"]]
```
　まずは読み込む指標を __indices__というリストに記述します。このリストの先頭にある指標について、予測を行います。この例では"N225" (日経平均)と"NYA" (NYSE総合)を入力として用い、N225の予測を行います。
　次にCSVファイルの読み込みを行います。__[indicesの指数].csv__ という名前のファイルを読み込み、__df_dict__ という辞書に格納します。細かいことを言うと、__indices__ の中身は __df_dict__ の"key"と同じなので、以降は __indices__ は必要なさそうですが、Pythonの辞書は順序を保存しませんのでこの後でも __indices__ を使っています ("OrderedDict"は順序を保存する辞書です)。
　また __indices[0]__ に相当するテーブルの __Close__/__Open__ を計算して、__df_ratio__ というテーブルに格納しています。以後、このテーブルのindex (日付)が日時の基準となります。
　次に各テーブルに __normalize_table__ を適用し、__start__ から __end__ までの分を取り出します。
　最後に __X_train__、__X_test__ の入力に使う列のリストを __columns_input__ として作成します。


# 4. 実行結果
　__(修正後の実行結果に差し替えてあります)。__
　実行前にデータをダウンロードしておきます。今回のプログラムを実行するためには日経平均のデータを __N225.csv__、NYSE総合のデータを __NYA.csv__ としてプログラムと同じフォルダに保存しておきます。
　前回までのプログラムとの整合性を見るために、まず

```php:dummy04.py
indices = ["N225"]
```
としてN225のデータだけで性能評価しましょう。すると

```
Indices: ['N225']
Start date: 2008-01-01
End date: 2017-12-31
 threshold: 0.10 %
 n_estimators: 5
 n_splits: 10
   positive training accuracy: 0.943
   positive test accuracy: 0.512±0.030
   positive mean gain: -0.020±0.029 %
   negative training accuracy: 0.941
   negative test accuracy: 0.523±0.043
   negative mean gain: 0.011±0.029 %
```
という結果になりました。前回の結果とは大きく違っていないので、同じように動作していそうです。次に

```php:dummy05.py
indices = ["N225", "NYA"]
```
として実行します。すると

```
Indices: ['N225', 'NYA']
Start date: 2008-01-01
End date: 2017-12-31
 threshold: 0.10 %
 n_estimators: 5
 n_splits: 10
   positive training accuracy: 0.950
   positive test accuracy: 0.520±0.041
   positive mean gain: 0.015±0.050 %
   negative training accuracy: 0.952
   negative test accuracy: 0.533±0.028
   negative mean gain: 0.031±0.029 %
```
という結果になりました。上昇側の平均利益が 0.015±0.050 %、下落側の平均利益が 0.031±0.029 %となり、下側は標準偏差と同じくらいの平均利益となっています。標準偏差と同程度なので「確実に利益が出る」とは言えませんが、「利益が期待できる」くらいの感じでしょうか。

　これでようやく利益を期待できそうなモデルができたので、次回はパラメーターを振ってさらに利益を大きくできないか検討してみることにします。

# ソース

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

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

def read_table(index='N225'):
    csv = index + '.csv'
    df = pd.read_csv(csv, na_values=["null"])
    df["Date"] = pd.to_datetime(df["Date"])
    df = df.set_index("Date")
    df = df[["Open", "High", "Low", "Close"]]
    return df.dropna()

def normalize_table(table, index='N225'):
    suffix_list = ["Open", "High", "Low", "Close"]
    table["Open"] /= table["Close"]
    table["High"] /= table["Close"]
    table["Low"] /= table["Close"]
    table["Close"] = 1.0
    table = table.resample("D").ffill()
    table = table[suffix_list].shift(1)
    table.columns = [index+"_"+suffix for suffix in suffix_list]
    return table


if __name__ == '__main__':
    ###########################################################
    # 読み込む指標を指定する。
    # 先頭の指標の予測をする。
    indices = ["N225", "NYA"]
    
    # 以下のブロックでパラメーターを指定する
    # この日以降のデータを使う
    start = '2008-01-01'

    # この日以前のデータを使う
    end = '2017-12-31'

    # 上昇または下落を判定するしきい値
    threshold = 0.001

    # RandomForestのn_estimator
    n_estimators = 5

    # KFoldの分割数
    n_splits = 10
    ##########################################################

    print("Indices: {}".format(indices))
    print("Start date: {}".format(start))
    print("End date: {}".format(end))
    print(" threshold: {:.2f} %".format(threshold * 100.0))
    print(" n_estimators: {}".format(n_estimators))
    print(" n_splits: {}".format(n_splits))
    
    df_dict = {}
    for index in indices:
        df_dict[index] = read_table(index)
    df_ratio = df_dict[indices[0]]['Close']/df_dict[indices[0]]['Open']
    df_ratio = df_ratio.to_frame("Ratio")
    
    for index in indices:
        df_dict[index] = normalize_table(df_dict[index], index=index)
    df = pd.concat([df_dict[index] for index in indices] + [df_ratio], axis=1, join='inner')
    df = df[start:end]

    columns_input = [index + suffix for index in indices\
                        for suffix in ["_Open", "_High", "_Low", "_Close"]]

    acc_train = {'positive': [], 'negative': []}
    acc_test  = {'positive': [], 'negative': []}
    gain      = {'positive': [], 'negative': [], 'total': {}}

    kf = KFold(n_splits=n_splits)

    for train, test in kf.split(df):
        X_train, r_train = df[columns_input].iloc[train], df["Ratio"].iloc[train]
        X_test,  r_test  = df[columns_input].iloc[test], df["Ratio"].iloc[test]
        X_train, X_test = X_train.values, X_test.values

        clf = RandomForestClassifier(n_estimators=n_estimators)
        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)
            acc_train[polarity].append(clf.score(X_train, y_train))
            acc_test[polarity].append(clf.score(X_test, y_test))
            if polarity == "positive":
                gain[polarity].append((clf.predict(X_test) * (r_test - 1)).mean())
            else:
                gain[polarity].append((clf.predict(X_test) * (1 - r_test)).mean())

    for polarity in ["positive", "negative"]:
        print('   {} training accuracy: {:.3f}'.format(polarity, np.array(acc_train[polarity]).mean()))
        print('   {} test accuracy: {:.3f}±{:.3f}'.format(polarity, 
                                                           np.array(acc_test[polarity]).mean(),
                                                           np.array(acc_test[polarity]).std()))
        print('   {} mean gain: {:.3f}±{:.3f} %'.format(polarity, 
                                                 np.array(gain[polarity]).mean()*100,
                                                 np.array(gain[polarity]).std()*100))



```