QuantXで機械学習を用いた株価予測（特徴量エンジニアリング編）

# はじめに
こんにちは、材料工学修士1年の学生です。鉄鋼材リサイクルへの画像解析応用を目指し、春から海外大のコンピュータサイエンス専攻へ交換留学しています。
この度Smart Tradeでインターンすることになり、株価予測モデルを作ってみました。

# 今日のお題
QuantX Factoryでアルゴリズムを作ってみる。
先人の方の機械学習モデル（ランダムフォレストモデル）を参考に、特徴量エンジニアリングを追加してみました。完成した[コード](https://factory.quantx.io/developer/b851b0dccf024e33a24c49462fe18990)はこちら。

# QuantX Factoryとは
smart trade社が無料で提供しているシステムトレードアルゴリズムをブラウザ上でpythonを使って開発できるプラットフォーム。

以下の理由から、pythonで株価予測をやってみたい人におすすめです。
- Web上のプラットフォームのため開発環境構築が不要であること
- 各銘柄の終値や高値などの各種データが用意されておりデータセットの用意不要
- 開発したトレードアルゴリズムは審査を通過すればQuantX Storeで販売可能
- 新エンジン移行に伴い、終値、始値、指値での注文が選択可能になるなど、より多様な戦略を組めるようになったこと

データの収集、更新だけに大量の時間を持っていかれる時代もこれにて終わりですね。
QuantXの詳しい使い方については[こちら](https://factory.quantx.io/handbook/ja/)。

# 本題
## １、目的と手法
- 既存のランダムフォレストモデルを使用した株価予測モデルへ、追加で特徴量エンジニアリングを実施し精度向上を試みます。
- エンジン移行に伴って、新エンジンにて使用可能なアルゴリズムへ修正を行います。
（データが命の機械学習の使用において、長期間バックテストの実行可能な新エンジンへの移行は重要です。）

### 参考にする株価予測モデル
機械学習を用いた株価予測モデルとして、**@sarubobo29**さんの[QuantXで機械学習使って株価予測してみたpart1](https://qiita.com/sarubobo29/items/e8219c2c30d500de2bf9)を使用させていただきました。「教師あり学習を使って、10日後の株価が現在の価格よりも「上がっているか(下がっているか)どうか」を予測するアルゴリズム。」になります。

ちなみに**@sarubobo29**さんの[QuantXで機械学習使って株価予測してみたpart2](https://qiita.com/sarubobo29/items/b94ccb6c1305c157fafb)へも適用を試みましたが、こちらは相性があまりよくありませんでした。


### 実施する特徴量エンジニアリング
機械学習モデル精度向上のため、以下の二点で特徴量エンジニアリングを行います。

#### i. 交互作用特徴量と多項式特徴量の生成
今回は交互作用特徴量と多項式特徴量を生成し、特徴量として追加していきます。交互作用特徴量とは特徴量の交互作用を表すための特徴量で、特徴量同士の積で表します。例えば、x1、x2という特徴量がある場合、x1 * x2という特徴量を新たに作成します。多項式特徴量とは特徴量の多項式表現を新たな特徴量とし、追加する方法です。例えば、x1、x2の特徴量の2次の表現を作ると、x1^2、x2^2という特徴量が新たに作成されます。

#### ii. 自動特徴量選択
訓練データへの過学習を防ぐために、自動特徴量選択を行います。自動特徴量選択とは、機械学習のモデルを使用において、全特徴量の中から有効な特徴量の組み合わせを自動で探索するプロセスのことを表します。今回はその中でも決定木ベースモデルなどの教師あり学習モデルを1つ用いて、特徴量の重要性を判断し重要なものだけを残す、モデルベース特徴量選択を行います。学習モデルとしてランダムフォレストモデルを用い、特徴量選択を行いました。

## ２、コード
### 修正点
#### 準備段階での修正点
scikit-learnから新たに、特徴量生成時に使用するPolynomialFeatures関数と、特徴量選択時に使用するSelectFromModel関数をインポートします。

```
from sklearn.feature_selection import SelectFromModel
from sklearn.preprocessing import PolynomialFeatures
```

#### 特徴量エンジニアリングを行う上での追加点(下で個別に説明します。)

```
# 特徴量の追加
poly = PolynomialFeatures(degree = 3)
poly.fit(df_x)
df_x = poly.transform(df_x)

# ここから手順3
# 学習データとテストデータに分割(前半5割を学習データ，後半5割をテストデータ)
X_train, X_test, y_train, y_test = train_test_split(df_x, y, train_size=0.5, shuffle = False)

#　特徴量の選択（モデルはランダムフォレスト）
select = SelectFromModel(RandomForestClassifier(n_estimators = 300, random_state = 1), threshold = "1.40*median")
select.fit(X_train, y_train)
X_train = select.transform(X_train)
X_test = select.transform(X_test)

# 学習データを使って学習(モデルはランダムフォレスト)
clf = RandomForestClassifier(random_state=1, n_estimators = 100, 
max_leaf_nodes = 10, max_depth=6, max_features=None)
clf = clf.fit(X_train, y_train)
```

まず、今回はPolynomialFeatures関数を用いて、3次以下の交互作用特徴量と、多項式特徴量を新たな特徴量として追加します。degree値を設定することで、その値以下の次元の交互作用特徴量と多項式特徴量を、生成することができます。


```
# 特徴量の追加
poly = PolynomialFeatures(degree = 3)
poly.fit(df_x)
df_x = poly.transform(df_x)
```

次に、特徴量選択を行います。SelectFromModelを用い、特徴量を選択します。今回は特徴量の重要性を判断するモデルとして、RandomForestClassifierを使用しました。特徴量選択では、ある基準を設定し、それ以上の重要性を持つ特徴量のみを選択します。今回は、この特徴量選択の基準として、median（中央値）を用い、threshold = "1.40*median"を採用しています。

```
#　特徴量の選択（モデルはランダムフォレスト）
select = SelectFromModel(RandomForestClassifier(n_estimators = 300, random_state = 1), threshold = "1.40*median")
select.fit(X_train, y_train)
X_train = select.transform(X_train)
X_test = select.transform(X_test)
```
#### 新エンジンでの使用のための修正点
シグナル名を"market:sig"へ修正し、新エンジンでも使用可能にします。

```
return {
  "market:sig": market_sig,
}
```

#### 訓練データ、テストデータのスコアをデバックするコードの追加
特徴量の追加や選択により、過剰適合や適合不足が起きることを防ぐために、訓練データ、テストデータに対する精度をデバック表示させます。

```
ctx.logger.debug("Training score: {:.3f}".format(clf.score(X_train, y_train)))
ctx.logger.debug("Test score: {:.3f}".format(clf.score(X_test, y_test)))
```

### 全体のコード
```
# Sample Algorithm
# ライブラリーのimport
# 必要ライブラリー
import maron
import maron.signalfunc as sf
import maron.execfunc as ef
# 追加ライブラリー
# 使用可能なライブラリに関しましては右画面のノートをご覧ください①
import pandas as pd
import numpy as np
import talib as ta
from sklearn.model_selection import train_test_split
from sklearn.ensemble import RandomForestClassifier
from sklearn.feature_selection import SelectFromModel
from sklearn.preprocessing import PolynomialFeatures

# オーダ方法(目的の注文方法に合わせて以下の2つの中から一つだけコメントアウトを外してください)
# オーダー方法に関しましては右画面のノートをご覧ください②
#ot = maron.OrderType.MARKET_CLOSE # シグナルがでた翌日の終値のタイミングでオーダー
ot = maron.OrderType.MARKET_OPEN   # シグナルがでた翌日の始値のタイミングでオーダー
#ot = maron.OrderType.LIMIT        # 指値によるオーダー

# 銘柄、columnsの取得
# 銘柄の指定に関しては右画面のノートをご覧ください③
# columnsの取得に関しては右画面のノートをご覧ください④
def initialize(ctx):
  # 設定
  ctx.logger.debug("initialize() called")
  ctx.configure(
    channels={               # 利用チャンネル
      "jp.stock": {
        "symbols": [
          "jp.stock.1321", #日経225連動型上場投資信託
        ],
        "columns": [
          #"close_price",     # 終値
          "close_price_adj", # 終値(株式分割調整後)
          #"volume_adj",     # 出来高
          #"txn_volume",     # 売買代金
        ]
      }
    }
  )
  
  # シグナル定義
  def _my_signal(data):
    # 欠損値補間
    cp = data["close_price_adj"].fillna(method='ffill')
    
    # ctx.logger.debug(cp)
       
    # ここから手順1(教師データの生成)
    # 教師データを入れる型を用意
    y = pd.DataFrame(0, index = cp.index, columns = cp.columns)
    y = cp.diff(-10)    # ここで何日後を予測するかを決めれる
    y = y['jp.stock.1321']<0    # True, Falseに変換
    
    # ここから手順2(株価データから特徴抽出)
    # 特徴量を入れるための準備
    rsi7 = pd.DataFrame(0, index=cp.index, columns=["rsi7"])
    sma14 = pd.DataFrame(0,index=cp.index, columns=["sma14"])
    sma7 = pd.DataFrame(0,index=cp.index, columns=["sma7"])
    sma7_diff = pd.DataFrame(0,index=cp.index, columns=["sma7_diff"])
    sma14_diff = pd.DataFrame(0,index=cp.index, columns=["sma14_diff"])
    sma_diff = pd.DataFrame(0,index=cp.index, columns=["sma_diff"])
    bb_mid = pd.DataFrame(0,index=cp.index, columns=["bb_mid"])
    bb_up = pd.DataFrame(0,index=cp.index, columns=["bb_up"])
    bb_low = pd.DataFrame(0,index=cp.index, columns=["bb_low"])
    bb_diff = pd.DataFrame(0,index=cp.index, columns=["bb_diff"])
    mom = pd.DataFrame(0,index=cp.index, columns=["mom"])
    
    for (sym,val) in cp.items():
      rsi7['rsi7'] = ta.RSI(cp[sym].values.astype(np.double), timeperiod=7)
      sma14['sma14'] = ta.SMA(cp[sym].values.astype(np.double), timeperiod=14)
      sma7['sma7'] = ta.SMA(cp[sym].values.astype(np.double), timeperiod=7)
      bb_up['bb_up'], bb_mid['bb_mid'], bb_low['bb_low'] = ta.BBANDS(cp[sym].values.astype(np.double), timeperiod=14)
      mom['mom'] = ta.MOM(cp[sym].values.astype(np.double), timeperiod=14)
      
    sma7_diff['sma7_diff'] = sma7.diff()
    sma14_diff['sma14_diff'] = sma14.diff()
    sma_diff['sma_diff'] = sma7['sma7'] - sma14['sma14']
    bb_diff['bb_diff'] = bb_mid['bb_mid'] - bb_up['bb_up']
    
    # 特徴量の結合
    df_x = pd.concat([rsi7, sma7_diff, sma14_diff, sma_diff, 
    bb_diff, mom], axis = 1)
    ctx.logger.debug(df_x)

    # はじめのNAN削除
    df_x = df_x[14:]
    y = y[14:]
    
    # 特徴量の追加
    poly = PolynomialFeatures(degree = 3)
    poly.fit(df_x)
    df_x = poly.transform(df_x)
    
    # ここから手順3
    # 学習データとテストデータに分割(前半5割を学習データ，後半5割をテストデータ)
    X_train, X_test, y_train, y_test = train_test_split(df_x, y, train_size=0.5, shuffle = False)
    
    #　特徴量の選択（モデルはランダムフォレスト）
    select = SelectFromModel(RandomForestClassifier(n_estimators = 300, random_state = 1), threshold = "1.40*median")
    select.fit(X_train, y_train)
    X_train = select.transform(X_train)
    X_test = select.transform(X_test)
    
    # 学習データを使って学習(モデルはランダムフォレスト)
    clf = RandomForestClassifier(random_state=1, n_estimators = 100, 
    max_leaf_nodes = 10, max_depth=6, max_features=None)
    clf = clf.fit(X_train, y_train)

    # ここから手順4
    # テストデータを使って予測
    pred = clf.predict(X_test)
    
    # market_sigのインデックスを合わせる
    test = np.ones(len(X_train)+14, dtype=np.bool) * False
    market_sig = np.hstack((test, pred))
    market_sig = pd.DataFrame(data = market_sig, columns=cp.columns, index=cp.index)
    
    # ctx.logger.debug(market_sig)feature = clf.feature_importances_
    # ctx.logger.debug(feature)
    
    feature = clf.feature_importances_
    ctx.logger.debug(market_sig)
    
    # normalize signals
    market_sig[market_sig == True] = 1.0
    market_sig[market_sig == False] = 0.0
    ctx.logger.debug(market_sig)
    
    # スコア表示
    ctx.logger.debug("Training score: {:.3f}".format(clf.score(X_train, y_train)))
    ctx.logger.debug("Test score: {:.3f}".format(clf.score(X_test, y_test)))


    return {
      "market:sig": market_sig,
    }
      
  # シグナル登録
  ctx.regist_signal("my_signal", _my_signal)

def handle_signals(ctx, date, current):
  '''
  current: pd.DataFrame
  '''
  # 1321を取引するため
  bull = ctx.getSecurity("jp.stock.1321")
  
  # シグナルを取得
  market_signal = current["market:sig"][0]

  if market_signal == 1.0:
    bull.order_target_percent(0.5, comment="BULL BUY")
  else:
    bull.order_target_percent(0, comment="BULL SELL")
```
## ３、テスト結果
<img width="819" alt="スクリーンショット 2019-07-04 11.01.27.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/383100/26f31be0-9a15-2314-db39-97de157be580.png">

損益率に着目すると、同期間でのバックテストで、14.60% → 22.60%へ上昇することができました。SharpeRatioなども改善されたようです。

<img width="301" alt="スクリーンショット 2019-07-03 19.46.06.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/383100/49c4e598-6af4-8585-88cb-5b19394edd90.png">

訓練データ、テストデータへのスコアをみてみると、訓練データへフィットしすぎており、過学習が起きていることがわかります。追加する特徴量の数や、そもそものランダムフォレストモデルの変数等を操作しましたが、多少の改善しか達成できず、今後の課題です。


## ４、改善案


- モデル面
過学習を防ぐための、各モデルのパラメータ調整
その他のモデルの使用や、複数モデルのアンサンブル

- 特徴量
ベーシックな指標を特徴量として追加

- その他
新エンジンで、より長期でのバックテストを行いtrainデータを増量

## ５、最後に
記事を読んでいただき、ありがとうございました。まだまだ知識が足りず、ミスなどあればぜひご指摘頂ければ幸いです。より良いアルゴリズム作成に向け、全力で取り組んでいきます！