joda!!

jodaによるプログラミング芸

仮想通貨MLBot つくってみた い ② テンプレートコード

はじめに

5倍が
f:id:jodawithforce:20220107210545p:plain

20倍に!?!?!
f:id:jodawithforce:20220107210630p:plain

こんにちは。 jodaと申します。

機械学習を使用した仮想通貨botをつくってみた(い) ということで、
まずは機械学習を利用して、勝てるロジックを見つけていこうと思います。

簡単な自己紹介をさせていただきますと、
自分はここ2年ほど、C++をベースとした為替の自動売買システム(EA)を作って稼働させておりました。 (結果は微妙です)
心機一転、世を賑わせる仮想通貨にも手を出そう!
ということで、仮想通貨botをつくっていこうと思ったのですが、
肝心のトレードロジックの検証方法がわからない状況でした。
今回は、まずpnadasを使って検証を行います。
その後、冒頭の画像のように、 機械学習を使ってトレードロジックを改善したいと思います。
(実際の執行が冒頭の画像の様になるとは限りません)


概要

こちらの記事は、
jodawithforce.hatenablog.com
の改善版となっております。

本記事では、richmanbtcさんのチュートリアルを参考に、
OHLCVデータを分析する際に、テンプレートとなるコードを記述いたします。

トレードにおけるエッジ(収益機会)は、広まると消滅してしまいます。
エッジ保護のため、実際にモデルに行った変更に関して、詳細な記述はせず、
モデル改善による結果を記載いたします。


目次


目的

仮想通貨自動売買で、機械学習を利用し、リターン予測を向上する。


実行環境

Google Colaboratory
Python 3.7.12
crypto-data-fetcher 0.0.17
tblib 1.7.0


ディレクトリ構成

最終的に以下のようなディレクトリ構成になります。

Mydrive/
 ┗Colab Notebooks/
    ┗deliverables/
         ┣df_y.pkl
         ┣df_ohlcv.pkl
         ┣df_fit.pkl
         ┣df_features.pkl
         ┣df_for_learn.pkl
         ┗deliverables.ipynb


設計

こちらの記事をベースにさせていただきまして、試行錯誤していきます。
機械学習で仮想通貨取引をするチュートリアル。python - Qiita

  1. 現在レートから一定距離はなれたところに指値注文を入れる。
  2. 特徴量としてテクニカル指標、目的変数として損益を用い、モデルを学習させる。
  3. 注文結果の予測が正の値である時、注文を入れるようにする。


流れ

  1. Google Driveのマウント
  2. ライブラリのインストール、インポート
  3. crypto-data-fetcher を利用して OHLCVを取得
  4. 手数料カラムを追加
  5. talibを用いてテクニカル分析を行い、特徴量を作成していく
  6. 目的変数の計算
  7. モデルの学習とOOS予測値計算
  8. テストデータで検証

となってございます。


コード作成

それではコードを作成していきましょう。
以下のコードに特徴量を足したり、サンプリングしたり様々な変更を加えてモデル改善いたします。
変更が行われること前提として、特徴量の増加などに、対応しやすいような記述を心がけます。
コードは、

github.com

で公開しております。 (コードがおいてあるだけで恐縮です、、、、)






Google Driveのマウント

以下のコードを実行しないとGoogle Drive上のファイルが利用できません。
遷移したりリンクが出力されたりするので、指示に従いましょう。

from google.colab import drive
drive.mount('/content/drive')





ライブラリのインストール

まずは必要ライブラリをインストールしていきます。

以下のライブラリをインストールしていきます。

  • crypto_data_fetcher

  • talib

  • lightgbm

  • optuna

crypto_data_fetcherはrichmanbtcさんのgitからインストールさせていただきました。
talibに関してはそのままpipインストールできませんでしたので、
以下の記事を参考にさせていただきました。

Python:TA-Libでテクニカル分析、Plotlyでローソク足の描画 - Investment Tech Hack
lightgbmは、特徴量重要度を表示しようとするとエラーになってしまったので、
アップデートします。
コードは以下のようになります。

!pip install ccxt 
!pip install "git+https://github.com/richmanbtc/crypto_data_fetcher.git@v0.0.17#egg=crypto_data_fetcher"
!curl -L http://prdownloads.sourceforge.net/ta-lib/ta-lib-0.4.0-src.tar.gz -O && tar xzvf ta-lib-0.4.0-src.tar.gz
!cd ta-lib && ./configure --prefix=/usr && make && make install && cd - && pip install ta-lib
!pip install -U lightgbm
!pip install optuna  

その他のライブラリに関してはすでにGoogleColaboratoryにインストールされている模様。
ライブラリのインストールが終了しましたので、 インポートしてみます。

from datetime import datetime
import calendar
from pprint import pprint
import pytz

import ccxt
from crypto_data_fetcher.gmo import GmoFetcher
import joblib
import lightgbm as lgb
import matplotlib.pyplot as plt
import numba
import numpy as np
import pandas as pd
from scipy.stats import ttest_1samp
import scipy.stats
import seaborn as sns
import talib
import optuna
import lightgbm as lgb_o

from sklearn.metrics import accuracy_score, confusion_matrix, classification_report
from sklearn.ensemble import BaggingRegressor
from sklearn.metrics import confusion_matrix
from sklearn.linear_model import RidgeCV, LinearRegression
from sklearn.model_selection import cross_val_score, KFold, TimeSeriesSplit
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_squared_error 
from sklearn.metrics import r2_score 

sns.set()

エラーが出なければOKです。 エラーが出た場合、足りないライブラリを入れていきます。






OHLCVの取得

まずは解析するデータを取得していきます。 OHLCVとは、

一定期間(1分、5分、1時間など)の,
始値(open)、高値(high)、安値(low)、終値(close)、出来高(volume)
を表すデータです。

crypto_data_fetcherは、binance bybit ftx gmo などの
価格データ等を取得することができます。

詳細はこちら

(当初はccxtを用いて、bitgetのOHLCVを取得しておりました。
しかしながら、 ccxtは 500本までという制限があり
データ分析を行うには足りなかったのでこちらを利用させていただいています)

memory = joblib.Memory('/tmp/gmo_fetcher_cache', verbose=0)
fetcher = GmoFetcher(memory=memory)

# GMOコインのBTC/JPYレバレッジ取引 ( https://api.coin.z.com/data/trades/BTC_JPY/ )を取得
# 初回ダウンロードは時間がかかる

#ダウンロードは1度で良いので、ランタイム再起動時に再度データを取得しないようにコメントアウト

df = fetcher.fetch_ohlcv(
    market='BTC_JPY', # 市場のシンボルを指定
    interval_sec=15 * 60, # 足の間隔を秒単位で指定。この場合は15分足
)

df.to_pickle('/content/drive/My Drive/Colab Notebooks/deliverables/df_ohlcv.pkl')





maker手数料カラムを追加

feeの単位には注意が必要。
%ではなく、割合表記にしないと、目的変数の計算でバグります。
詳細は、私の別記事にて。

jodawithforce.hatenablog.com

GMOさんのメイカー手数料は、現在は0となっております。

過去にはマイナス手数料であったこともありますが、
モデルの学習自体には関係ないため、一旦0で統一いたします。

df["fee"] = 0 #割合表記!!!!

作成されたDFは以下のようになりました。
(当時の時刻からの過去データを取得するので、
セルを実行するたびにDFは更新されていきます。)
f:id:jodawithforce:20220106000544p:plain






特徴量作成

テクニカル指標計算ライブラリのTA-Libを利用して特徴量を作成していきます。
talibは非常に便利で、関数一つでテクニカル指標を計算することができます。

例えば、以下のコードの
df['ADX'] = talib.ADX(high, low, close, timeperiod=14)
に着目してみます。
こちらはtalib.ADX()一つで、
dfのtimestampに対応するADXが算出できます。
引数には、高値(high)、安値(low)、終値(close)、そして算出期間を記述します。

どのテクニカル指標がモデルの予測精度を向上させるか不明ですので、
ひとまずメジャーなテクニカル指標をまとめてdfに追加することにします。

def calc_features(df):
    open = df['op']
    high = df['hi']
    low = df['lo']
    close = df['cl']
    volume = df['volume']
    
    orig_columns = df.columns

    hilo = (df['hi'] + df['lo']) / 2
    df['BBANDS_upperband'], df['BBANDS_middleband'], df['BBANDS_lowerband'] = talib.BBANDS(close, timeperiod=5, nbdevup=2, nbdevdn=2, matype=0)
    df['BBANDS_upperband'] -= hilo
    df['BBANDS_middleband'] -= hilo
    df['BBANDS_lowerband'] -= hilo
    df['DEMA'] = talib.DEMA(close, timeperiod=30) - hilo
    df['EMA'] = talib.EMA(close, timeperiod=30) - hilo
    df['HT_TRENDLINE'] = talib.HT_TRENDLINE(close) - hilo
    df['KAMA'] = talib.KAMA(close, timeperiod=30) - hilo
    df['MA'] = talib.MA(close, timeperiod=30, matype=0) - hilo
    df['MIDPOINT'] = talib.MIDPOINT(close, timeperiod=14) - hilo
    df['SMA'] = talib.SMA(close, timeperiod=30) - hilo
    df['T3'] = talib.T3(close, timeperiod=5, vfactor=0) - hilo
    df['TEMA'] = talib.TEMA(close, timeperiod=30) - hilo
    df['TRIMA'] = talib.TRIMA(close, timeperiod=30) - hilo
    df['WMA'] = talib.WMA(close, timeperiod=30) - hilo

    df['ADX'] = talib.ADX(high, low, close, timeperiod=14)
    df['ADXR'] = talib.ADXR(high, low, close, timeperiod=14)
    df['APO'] = talib.APO(close, fastperiod=12, slowperiod=26, matype=0)
    df['AROON_aroondown'], df['AROON_aroonup'] = talib.AROON(high, low, timeperiod=14)
    df['AROONOSC'] = talib.AROONOSC(high, low, timeperiod=14)
    df['BOP'] = talib.BOP(open, high, low, close)
    df['CCI'] = talib.CCI(high, low, close, timeperiod=14)
    df['DX'] = talib.DX(high, low, close, timeperiod=14)
    df['MACD_macd'], df['MACD_macdsignal'], df['MACD_macdhist'] = talib.MACD(close, fastperiod=12, slowperiod=26, signalperiod=9)
    # skip MACDEXT MACDFIX たぶん同じなので
    df['MFI'] = talib.MFI(high, low, close, volume, timeperiod=14)
    df['MINUS_DI'] = talib.MINUS_DI(high, low, close, timeperiod=14)
    df['MINUS_DM'] = talib.MINUS_DM(high, low, timeperiod=14)
    df['MOM'] = talib.MOM(close, timeperiod=10)
    df['PLUS_DI'] = talib.PLUS_DI(high, low, close, timeperiod=14)
    df['PLUS_DM'] = talib.PLUS_DM(high, low, timeperiod=14)
    df['RSI'] = talib.RSI(close, timeperiod=14)
    df['STOCH_slowk'], df['STOCH_slowd'] = talib.STOCH(high, low, close, fastk_period=5, slowk_period=3, slowk_matype=0, slowd_period=3, slowd_matype=0)
    df['STOCHF_fastk'], df['STOCHF_fastd'] = talib.STOCHF(high, low, close, fastk_period=5, fastd_period=3, fastd_matype=0)
    df['STOCHRSI_fastk'], df['STOCHRSI_fastd'] = talib.STOCHRSI(close, timeperiod=14, fastk_period=5, fastd_period=3, fastd_matype=0)
    df['TRIX'] = talib.TRIX(close, timeperiod=30)
    df['ULTOSC'] = talib.ULTOSC(high, low, close, timeperiod1=7, timeperiod2=14, timeperiod3=28)
    df['WILLR'] = talib.WILLR(high, low, close, timeperiod=14)

    df['AD'] = talib.AD(high, low, close, volume)
    df['ADOSC'] = talib.ADOSC(high, low, close, volume, fastperiod=3, slowperiod=10)
    df['OBV'] = talib.OBV(close, volume)

    df['ATR'] = talib.ATR(high, low, close, timeperiod=14)
    df['NATR'] = talib.NATR(high, low, close, timeperiod=14)
    df['TRANGE'] = talib.TRANGE(high, low, close)

    df['HT_DCPERIOD'] = talib.HT_DCPERIOD(close)
    df['HT_DCPHASE'] = talib.HT_DCPHASE(close)
    df['HT_PHASOR_inphase'], df['HT_PHASOR_quadrature'] = talib.HT_PHASOR(close)
    df['HT_SINE_sine'], df['HT_SINE_leadsine'] = talib.HT_SINE(close)
    df['HT_TRENDMODE'] = talib.HT_TRENDMODE(close)

    df['BETA'] = talib.BETA(high, low, timeperiod=5)
    df['CORREL'] = talib.CORREL(high, low, timeperiod=30)
    df['LINEARREG'] = talib.LINEARREG(close, timeperiod=14) - close
    df['LINEARREG_ANGLE'] = talib.LINEARREG_ANGLE(close, timeperiod=14)
    df['LINEARREG_INTERCEPT'] = talib.LINEARREG_INTERCEPT(close, timeperiod=14) - close
    df['LINEARREG_SLOPE'] = talib.LINEARREG_SLOPE(close, timeperiod=14)
    df['STDDEV'] = talib.STDDEV(close, timeperiod=5, nbdev=1)

    return df

df = pd.read_pickle('/content/drive/My Drive/Colab Notebooks/deliverables/df_ohlcv_with_fee.pkl')
df = df.dropna()
df = calc_features(df)
display(df)
df.to_pickle('/content/drive/My Drive/Colab Notebooks/deliverables/df_features.pkl')

作成されたDFは以下のようになります。
テクニカル指標は、DF何行分かの価格データ(終値、高値,,,etc)を
利用して作られるものが多いです。
DFの上部は、テクニカル指標作成に必要な終値等の価格データの数が足りていないため、 NaNとなるものがでてきます。
f:id:jodawithforce:20220106003541p:plain






カテゴリ変数のカラムの保存

後ほど標準化を行います。
カテゴリ変数のカラムには標準化の必要がないため、
カテゴリ変数のカラムの配列を作成しておきます。
ユニークな値が2つしかないカラムを取得します。

category_col = df.nunique()
category_col = category_col[category_col==2]
category_col.index





目的変数の計算

こちらに関してはかなりボリューミーになってしまいましたので、
別記事にまとめました。

jodawithforce.hatenablog.com

ざっくりいってしまうと、

  • 時間ごとに現在価格から、ATR/2だけ離れたところに指値をおきます。

  • 結果、得られる損益を目的変数とします。

という内容です。
モデルによる学習前の結果は以下のようになりました。

f:id:jodawithforce:20220106004423p:plain
図1:累積リターン

こちらは、赤線が売り、青線が買い、緑が合計のリターンを表しています。
割合で表記されています。
こちらのロジックで1BTCを売買していた場合、
20倍位にはなっていたということになりますね、、、、

matplotlibの日本語化がされていないようで、
警告が出てしまいました。
以下のように修正可能のようですが、 英語にして対応ます。 colab.research.google.com

指値を出す距離や、horizon、最大のポジション数など、
調整可能なパラメーターがあり、この辺が面白そうですが、
ひとまず学習させてどのような値になるか見てみます。






標準化

dfを標準化していきます。
もっとよい方法があるきがしますが、以下のやり方で標準化を行いました。
まずは、保存しておいたカテゴリ変数と、手数料、目的変数のカラムをdfから分離します。
(手数料カラムはすべて0なので、そのまま標準化すると変なことになった記憶があります)
その後、scipyをもちいて、dfを標準化します。
分離しておいたカラムをjoinして完了です。

また、この時点で不要なカラムを削除しておきます。

#モデルの学習とOOS予測値計算
df = pd.read_pickle('/content/drive/My Drive/Colab Notebooks/deliverables/df_y.pkl')
df = df.dropna()

#標準化が必要のないカラムについて、別で保存しておく
not_stz_col = list(category_col.index) + ["fee", "y_buy", "y_sell"]
not_stz = df[not_stz_col]
df = df.drop(not_stz_col, axis=1)

columns = df.columns.values
index = df.index

#DFを標準化
df = scipy.stats.zscore(df, ddof=1)
df = pd.DataFrame(df)
df.columns = columns
df.index = index
df = df.join(not_stz)

#カラムを特徴量と目的変数のみにしておく
drop_col =  ['op', 'hi', 'lo', 'cl',  'buy_price','sell_price', 'buy_fep', 
              'buy_fet', 'sell_fep', 'sell_fet', 'buy_executed', 
              'sell_executed', 'buy_cost', 'sell_cost', 'fee']

df = df.drop(drop_col, axis=1)





学習用データによるモデルの学習

学習用データを使って、モデルを評価します。
学習用データである程度結果が改善したら、後述のテストデータでモデルを評価します。

LightGBMについて

それではモデルの学習をおこなっていきます。
今回利用するモデルは、 LightGBM です。
「Kaggler」の上位6割以上が LightGBM を用いている、
という集計結果もある模様です。

LightGBMについて、別記事にまとめました。

jodawithforce.hatenablog.com


関数の作成

今回の場合、買いと売りで同じような処理を2回行います。
そのため、コードを関数化いたします。

まずは、実際に学習を行う関数です。
特徴量、目的変数、cvのインデックスを入力すると、
予測値、各Foldの特徴量重要度を格納した辞書を返します。

今回の場合、KFoldを利用するので、
forループのたびに特徴量重要度を保存しておく必要があります。

また、optunaを利用しておりますが、optunaは時間がかかります。
import文をimport lightgbm as lgb_oとすれば、通常のlightGBMを利用できます。

テストデータを入れることで、テストデータに対する予測を返します。

def my_cross_val_predict(X, y, cv_index=None, x_test=None, y_test=None):
    
    #変数の用意
    y_pred = np.zeros(len(X))
    y_pred[:] = np.nan

    #テストデータが存在する時
    if x_test is not None:
        y_test_pred = np.zeros(len(x_test))
    else:
        y_test_pred = None

    #固定するパラメータ
    params = {'objective': 'regression',
               'metric': 'rmse',
               'extra_trees': True,
               'random_seed': 0,
             } 

    evaluate_score_list = []
    importance_df = pd.DataFrame({"feature":[], "importance":[]})

    for trn_idx, val_idx in cv_index:
        x_trn, x_val = X.iloc[trn_idx], X.iloc[val_idx]
        y_trn, y_val = y.iloc[trn_idx], y.iloc[val_idx]

        # 訓練データ
        lgb_train = lgb_o.Dataset(x_trn, y_trn)
        # 評価データ
        lgb_eval = lgb_o.Dataset(x_val, y_val, reference=lgb_train)
        
        #モデルの学習
        # モデルの学習
        model = lgb_o.train(params,
                         train_set=lgb_train,
                         valid_sets=lgb_eval,
                         early_stopping_rounds=100,
                         verbose_eval=-1,
                         )


        #Foldごとに重要度の格納
        cols = X.columns.values
        f_importance = np.array(model.feature_importance(importance_type='gain')) # 特徴量重要度の算出
        importance_df = pd.concat((importance_df, pd.DataFrame(
                                        {"feature": cols, "importance": f_importance})))

        y_pred[val_idx] = model.predict(x_val)

        if x_test is not None:
            y_test_pred += model.predict(x_test)
            print(y_test_pred)
    
    #テストデータでの評価
    if y_test_pred is not None:
        print(y_test_pred)
        y_test_pred = y_test_pred / len(cv_index)
        print(y_test_pred)

    ret = {"y_pred": y_pred, 
            "importance_df": importance_df,
            "y_test_pred": y_test_pred} 

    return ret


dfの処理

ここまでで、ベースとなるdf、分析用の関数ができております。
以降、分析を行っていくのですが、
更にdfに様々な処理を加えていきます。
以下に、最もリターン予測の良かった際の処理を記述します。

ざっくり申し上げますと、
約定した時のリターンのみを学習させる、ということになります。
目的変数に0が多いという偏りがあった、かつ、
0は約定していないということなので、 予測がどうであろうがリターンに影響しない
と考えました。目的はリターンを最大化することなので、
以下のように、リターンが0でない行を抽出しました。

#損益が0の行を取り除いたDF(回帰用) 
df_reg_buy = df[df["y_buy"]!=0]
df_reg_sell = df[df["y_sell"]!=0]

#可視化
plt.hist(df_reg_buy["y_buy"], bins=100, color="lightblue", label="y_buy")
plt.xlim(-0.03, 0.03)
plt.title("y_buy")
plt.show()

plt.hist(df_reg_sell["y_sell"], bins=100, color="lightcoral", label="y_sell")
plt.xlim(-0.03, 0.03)
plt.title("y_sell")
plt.show()

可視化した目的変数の分布は以下のようになります。
正規分布に近い気がしますね。
f:id:jodawithforce:20220105064206p:plain f:id:jodawithforce:20220105064218p:plain


CV

dfを分割していきます。
時系列分析を行うわけではないのですが、
実際には過去データから未来を予測するということで、
cvにはKFoldを利用します。
3ヶ月に1回再学習するとして、テストデータを3ヶ月分にしようか迷いましたが、
テストデータが少なすぎる気がしたので、テストサイズは0.2としました。

#学習用と テスト用 に分割
#一応時系列として、過去から未来を予測する。
test_size = 0.2
train_df, test_df = train_test_split(df, test_size=test_size, shuffle=False)

#エントリーしたところのみのDF作成(y_buy!=0  y_sell!=0)
train_df_buy = train_df[train_df["y_buy"]!=0]
train_df_sell = train_df[train_df["y_sell"]!=0]

# KFold法 インデックス番号(0,1,2,3,,,,)で出てくるので注意! ilocをつかう
cv_indicies_buy = list(KFold().split(train_df_buy))
cv_indicies_sell = list(KFold().split(train_df_sell))


特徴量カラム名の配列を取得

特徴量選択を行う前は、すべての特徴量を利用します。
特徴量のカラム名を取得しておきます。

features = df.columns.values
features = features[(features!="y_buy") & (features!="y_sell")]
features


モデルの学習

実際にモデルの学習を行っていきます。
my_cross_val_predict()を使います。
はじめの1回目は、
buy_features = features
sell_features = features
として、すべての特徴量をつかって学習をおこないます。
特徴量選択を行う場合は、カラム名を記述していきます。

print("prediction of y_buy")
#Light_GBMの結果から、特徴量重要度を降順にソートしてリスト化する。
#重要度の低いものから省いて、特徴量選択を行う
buy_features = []

r_buy_reg = my_cross_val_predict(train_df_buy[buy_features], train_df_buy['y_buy'], cv_indicies_buy)

print(print("prediction of y_sell"))
#Light_GBMの結果から、特徴量重要度を降順にソートしてリスト化する。
#重要度の低いものから省いて、特徴量選択を行う
sell_features = []

r_sell_reg = my_cross_val_predict(train_df_sell[sell_features], train_df_sell['y_sell'], cv_indicies_sell)


モデルの評価1

評価指標の確認

まずは以下の関数を作成します。
回帰における観測値と予測値を引数とします。
R2,MSE,RMSEをdfにして出力します。
また、横軸に観測値、縦軸に予測値をとった散布図を出力します。

def calculate_scores(true, pred):
    """全ての評価指標を計算する

    Parameters
    ----------
    true (np.array)       : 実測値
    pred (np.array)       : 予測値

    Returns
    -------
    scores (pd.DataFrame) : 各評価指標を纏めた結果

    """
    mod = LinearRegression()

    plt.plot(true, true, color = 'red', label = 'x=y') # 直線y = x (真値と予測値が同じ場合は直線状に点がプロットされる)
    plt.scatter(true, pred) # 散布図のプロット
    plt.xlabel('true') # x軸ラベル
    plt.ylabel('pred') # y軸ラベル
    plt.title('true vs pred') # グラフタイトル

    df_true = pd.DataFrame(true)
    df_pred = pd.DataFrame(pred)

    mod.fit(df_true, df_pred)
    lin_reg = mod.predict(df_true)
    plt.plot(true, lin_reg, color = '#008000', linewidth=2)
    plt.show()

    scores = pd.DataFrame({'R2': r2_score(true, pred),
                             'MSE': mean_squared_error(true, pred),
                             'RMSE': np.sqrt(mean_squared_error(true, pred))},
                             index = ['scores'])
    return scores

以下のコードで、買いと売りの評価を出力していきます。

evaluations_buy_reg = calculate_scores(train_df_buy["y_buy"], r_buy_reg["y_pred"])
evaluations_sell_reg = calculate_scores(train_df_sell["y_sell"], r_sell_reg["y_pred"])

evaluations_buy_reg["label"] = ["buy_evaluations_validation"]
evaluations_sell_reg["label"] = ["sell_evaluations_validation"]

evaluations_df_validation = pd.concat([evaluations_buy_reg, evaluations_sell_reg])
evaluations_df_validation = evaluations_df_validation.set_index(["label"])
evaluations_df_validation

初期の学習では以下のような結果になりました。


f:id:jodawithforce:20220107131851p:plain
図2 : 予測値と観測値の散布図



f:id:jodawithforce:20220106010333p:plain
↑図3 : **回帰の評価指標**


累積リターンの比較

以下のようなコードを実行し、
予測値を利用しない場合と、予測値を利用したときのリターンを比較します。

#元のDFと、予測値を結合し、y_pred>0のときだけ取引した場合の結果を表示
train_df_buy['y_pred_buy'] = r_buy_reg["y_pred"]
train_df_sell['y_pred_sell'] = r_sell_reg["y_pred"]

train_df['y_pred_buy'] = train_df_buy['y_pred_buy']
train_df['y_pred_sell'] = train_df_sell['y_pred_sell']

train_df = train_df.fillna(0)

train_df['y_buy'].cumsum().plot(label='buy')
train_df['y_sell'].cumsum().plot(label='sell')
(train_df['y_buy'] + train_df['y_sell']).cumsum().plot(label='buy+sell')
plt.title('Cumulative return')
plt.legend(bbox_to_anchor=(1.05, 1))

train_df[train_df['y_pred_buy'] > 0]['y_buy'].cumsum().plot(label='buy_pred')
train_df[train_df['y_pred_sell'] > 0]['y_sell'].cumsum().plot(label='sell_pred')
(train_df['y_buy'] * (train_df['y_pred_buy'] > 0) + train_df['y_sell'] * (train_df['y_pred_sell'] > 0)).cumsum().plot(label='buy+sell_pred')
plt.title('Cumulative return')
plt.legend(bbox_to_anchor=(1.05, 1))
plt.show()

得られたグラフは以下のようになりました。
f:id:jodawithforce:20220106012415p:plain
図4 : 予測の有無によるリターンの比較

予測値を利用しない場合の累積リターン → 15.9
予測値を利用したときの累積リターン → 17.8
となり、およそ2改善しております。
15.9倍と17.8倍では大きく結果が異なるので、 この時点で機械学習の強力さが伺えます。
しかしながら、あくまでバリデーションのため、まだわかりません。


特徴量重要度の出力

まずは関数を作成していきます。
特徴量重要度のDFを引数に指定すると、横向き棒グラフを作成してくれる関数です。
値はソートされます。

def plot_feature_importance(df, title=None): 
    n_features = len(df)                              # 特徴量数(説明変数の個数) 
    df_plot = df.sort_values('importance')            # df_importanceをプロット用に特徴量重要度を昇順ソート 
    f_importance_plot = df_plot['importance'].values  # 特徴量重要度の取得 
    plt.figure(figsize=(6,15))
    plt.barh(range(n_features), f_importance_plot, align='center') 
    cols_plot = df_plot['feature'].values             # 特徴量の取得 
    plt.yticks(np.arange(n_features), cols_plot)      # x軸,y軸の値の設定
    plt.xlabel('Feature importance(%)')                  # x軸のタイトル
    plt.ylabel('Feature')                             # y軸のタイトル
    plt.title(title)
    plt.grid(True)
    plt.show()

my_cross_val_predict()は、特徴量重要度をFoldごとにconcatしたDFを返します。
こちらをgroupbyメソッドにより、特徴量ごとにまとめ、平均します。
plot_feature_importance()により、可視化されます。

def importance_groupby(importance_df, graph_title="feature importance", normalize=True):
    importance_df = importance_df.groupby("feature", as_index=False).mean().sort_values("importance", ascending=False)
    if normalize:
        importance_df["importance"] = (importance_df["importance"] / np.sum(importance_df["importance"])) * 100 #正規化(%)
    plot_feature_importance(importance_df, graph_title)
    
    return importance_df

関数を実行していきます。
np.set_printoptions(threshold=1000)
とすることで、配列を出力した際に、1000要素まで表示されます。
pprint(np.array(importance_buy_reg[importance_buy_reg["importance"]>border]["feature"]))
は、ソートされた特徴量重要度のうち、border以上のものを表示します。
出力された配列を、上述のbuy_featuressell_featuresにコピペで代入するだけで、
特徴量重要度のフィルタリングが可能です。
今回は、border = 1.5としました。

np.set_printoptions(threshold=1000)
border = 1.5

#買いにおける特徴量重要度
importance_buy_reg = importance_groupby(r_buy_reg["importance_df"])
#重要度を降順にソートし、特徴量選択を行う
pprint(np.array(importance_buy_reg[importance_buy_reg["importance"]>border]["feature"]))

#売りにおける特徴量重要度
importance_sell_reg =  importance_groupby(r_sell_reg["importance_df"])
#重要度を降順にソートし、特徴量選択を行う
pprint(np.array(importance_sell_reg[importance_sell_reg["importance"]>border]["feature"]))

得られた結果は以下のようになりました。

f:id:jodawithforce:20220106015235p:plain
↑図 : 5 買いの特徴量重要度を可視化


以下、買いにおける 特徴量重要度>1.5 のカラム名の配列です。

array(['LINEARREG', 'volume', 'TRIX', 'ADOSC', 'BBANDS_upperband',
       'CORREL', 'BETA', 'OBV', 'AD', 'NATR', 'STDDEV',
       'HT_PHASOR_quadrature', 'HT_DCPERIOD', 'HT_TRENDLINE',
       'STOCHF_fastk', 'AVGPRICE', 'ULTOSC', 'MACD_macdsignal', 'ATR',
       'APO', 'HT_PHASOR_inphase', 'STOCHRSI_fastd', 'AROON_aroondown',
       'ADXR', 'HT_TRENDMODE', 'ADX', 'MIDPOINT', 'TEMA',
       'STOCHRSI_fastk'], dtype=object)

f:id:jodawithforce:20220106015518p:plain
↑図 : 6 売りの特徴量重要度を可視化


以下、売りにおける 特徴量重要度>1.5 のカラム名の配列です。

array(['HT_PHASOR_quadrature', 'NATR', 'TRIX', 'HT_PHASOR_inphase',
       'HT_DCPERIOD', 'ADOSC', 'CORREL', 'AVGPRICE', 'STOCHRSI_fastd',
       'volume', 'LINEARREG', 'ULTOSC', 'HT_SINE_sine', 'BOP', 'OBV',
       'STOCH_slowd', 'BBANDS_upperband', 'ADXR', 'BBANDS_lowerband',
       'HT_TRENDLINE', 'ADX', 'PPO', 'MACD_macdsignal', 'BETA', 'PLUS_DM',
       'SAREXT', 'MIDPOINT', 'HT_TRENDMODE'], dtype=object)


モデルの評価2

今回の場合、回帰を使って予測しますが、実際の利用方法は分類です。
そのため、分類の評価指標を使ってモデルを評価していきます。

混同行列の作成と、評価指標の表示

機械学習による予測は、リターン予測がプラスならトレード、そうでないならトレードしない、
という使い方をします。
そのため、要素が0より大きい場合1、 そうでない場合0としてカラムを作成していきます。
その後、混同行列を作成し、ヒートマップを出力します。
また、classification_report()を利用して、正確性、精度、再現性、f値を出力します。

#2値分類したカラムを作成
train_df_buy["y_buy_bin"] = np.where(train_df_buy["y_buy"]>0, 1, 0)
train_df_buy["y_buy_pred_bin"] = np.where(train_df_buy["y_pred_buy"]>0, 1, 0)

print("買い 混同行列")
cm_buy = confusion_matrix(train_df_buy["y_buy_bin"], train_df_buy["y_buy_pred_bin"])
print(cm_buy)
cm_buy = pd.DataFrame(data=cm_buy, columns=['Predict Negative:0', 'Predict Positive:1'], 
                                           index=['Actual Negative:0', 'Actual Positive:1'])
cm_buy = cm_buy.reindex(index=['Actual Positive:1', 'Actual Negative:0'], columns=['Predict Positive:1', 'Predict Negative:0'])
sns.heatmap(cm_buy, annot=True, fmt='d', cmap='YlGnBu')

buy_clf_report = classification_report(train_df_buy["y_buy_bin"], train_df_buy["y_buy_pred_bin"])
print(buy_clf_report)


train_df_sell["y_sell_bin"] = np.where(train_df_sell["y_sell"]>0, 1, 0)
train_df_sell["y_sell_pred_bin"] = np.where(train_df_sell["y_pred_sell"]>0, 1, 0)

print("売り 混同行列")
cm_sell = confusion_matrix(train_df_sell["y_sell_bin"], train_df_sell["y_sell_pred_bin"])

cm_sell = pd.DataFrame(data=cm_sell, columns=['Predict Negative:0', 'Predict Positive:1'], 
                                           index=['Actual Negative:0', 'Actual Positive:1'])
cm_sell = cm_sell.reindex(index=['Actual Positive:1', 'Actual Negative:0'], columns=['Predict Positive:1', 'Predict Negative:0'])
sns.heatmap(cm_sell, annot=True, fmt='d', cmap='YlGnBu')

sell_clf_report = classification_report(train_df_sell["y_sell_bin"], train_df_sell["y_sell_pred_bin"])
print(sell_clf_report)

以下のような出力が得られます。

f:id:jodawithforce:20220106022548p:plain
↑図7 : 買い 混同行列


f:id:jodawithforce:20220106022636p:plain
↑図8: 買い 評価指標(分類)


f:id:jodawithforce:20220106022855p:plain
↑図9 : 売り 混同行列


f:id:jodawithforce:20220106022832p:plain
↑図10 : 売り 評価指標(分類)

混同行列を見ると、1(リターンがプラス)と予測することで、
一見予測がうまくいっているよに見えているだけで、あまりうまく分類できていない模様です。
検定を行っていないのでわかりませんが、f値も50を少し超える程度で、信頼するには心もとないですね。


テストデータによるモデルの評価

学習用データによるCVで、ある程度結果が改善したらテストデータで確認します。
my_cross_val_predict()の引数にテストデータを渡すだけなので、
学習用データの評価の際とほとんどコードは同じになります。

train_df, test_df = train_test_split(df, test_size=test_size, shuffle=False)

#エントリーしたところのみのDF作成(y_buy!=0  y_sell!=0)
train_df_buy = train_df[train_df["y_buy"]!=0]
train_df_sell = train_df[train_df["y_sell"]!=0]

# KFold法 インデックス番号(0,1,2,3,,,,)で出てくるので注意! ilocをつかう
cv_indicies_buy = list(KFold().split(train_df_buy))
cv_indicies_sell = list(KFold().split(train_df_sell))

r_buy_reg = my_cross_val_predict(train_df_buy[buy_features], train_df_buy['y_buy'], cv_indicies_buy, test_df[buy_features], test_df["y_buy"])
r_sell_reg = my_cross_val_predict(train_df_sell[sell_features], train_df_sell['y_sell'], cv_indicies_sell, test_df[sell_features], test_df["y_sell"])

test_df["y_pred_buy"] = r_buy_reg["y_test_pred"]
buy_evals = calculate_scores(test_df["y_buy"], test_df["y_pred_buy"])

test_df["y_pred_sell"] = r_sell_reg["y_test_pred"]
sell_evals = calculate_scores(test_df["y_sell"], test_df["y_pred_sell"])

##evaluation 
buy_evals["label"] = np.array(["buy_evaluations_test"])
sell_evals["label"] = np.array(["sell_evaluations_test"])

evaluations_df_test = pd.concat([buy_evals, sell_evals])
evaluations_df_test = evaluations_df_test.set_index(["label"])

evaluations_df_comparison = pd.concat([evaluations_df_validation, evaluations_df_test])
evaluations_df_comparison

print('毎時刻、y_predがプラスのときだけトレードした場合の累積リターン 比較')
test_df['y_buy'].cumsum().plot(label='buy')
test_df['y_sell'].cumsum().plot(label='sell')
(test_df['y_buy'] + test_df['y_sell']).cumsum().plot(label='buy+sell')

test_df[test_df['y_pred_buy'] > 0]['y_buy'].cumsum().plot(label='buy_use_pred')
test_df[test_df['y_pred_sell'] > 0]['y_sell'].cumsum().plot(label='sell_use_pred')
(test_df['y_buy'] * (test_df['y_pred_buy'] > 0) + test_df['y_sell'] * (test_df['y_pred_sell'] > 0)).cumsum().plot(label='buy+sell_use_pred')

plt.title('Cumulative return Comparison')
plt.legend(bbox_to_anchor=(1.05, 1))
plt.show()

#2値分類したカラムを作成
test_df_buy = test_df[test_df["y_buy"]!=0]  
test_df_buy["y_buy_bin"] = np.where(test_df_buy["y_buy"]>0, 1, 0)  
test_df_buy["y_buy_pred_bin"] = np.where(test_df_buy["y_pred_buy"]>0, 1, 0)  
print("買い 混同行列")  
cm_buy = confusion_matrix(test_df_buy["y_buy_bin"], test_df_buy["y_buy_pred_bin"])
print(cm_buy)
cm_buy = pd.DataFrame(data=cm_buy, columns=['Predict Negative:0', 'Predict Positive:1'], 
                                           index=['Actual Negative:0', 'Actual Positive:1'])
cm_buy = cm_buy.reindex(index=['Actual Positive:1', 'Actual Negative:0'], columns=['Predict Positive:1', 'Predict Negative:0'])
sns.heatmap(cm_buy, annot=True, fmt='d', cmap='YlGnBu')

buy_clf_report = classification_report(test_df_buy["y_buy_bin"], test_df_buy["y_buy_pred_bin"])
print(buy_clf_report)

test_df_sell = test_df[test_df["y_sell"]!=0]
test_df_sell["y_sell_bin"] = np.where(test_df_sell["y_sell"]>0, 1, 0)
test_df_sell["y_sell_pred_bin"] = np.where(test_df_sell["y_pred_sell"]>0, 1, 0)

print("売り 混同行列")
cm_sell = confusion_matrix(test_df_sell["y_sell_bin"], test_df_sell["y_sell_pred_bin"])

cm_sell = pd.DataFrame(data=cm_sell, columns=['Predict Negative:0', 'Predict Positive:1'], 
                                           index=['Actual Negative:0', 'Actual Positive:1'])
cm_sell = cm_sell.reindex(index=['Actual Positive:1', 'Actual Negative:0'], columns=['Predict Positive:1', 'Predict Negative:0'])
sns.heatmap(cm_sell, annot=True, fmt='d', cmap='YlGnBu')

sell_clf_report = classification_report(test_df_sell["y_sell_bin"], test_df_sell["y_sell_pred_bin"])
print(sell_clf_report)

OHLCVを分析するテンプレートコードは以上になります。
他のOHLCVデータを上記のコードに適用したり、目的変数の作成方法を変更したりして利用可能です。






モデルの改善

効率的市場仮説というものがあります。
簡単にいってしまうと、株価等の予測は不可能であるという学説です。
そのような学説があるデータを分析するにあたって、精度の高いモデルを作成することは容易ではないと考えております。
そのため、モデルを改善するために、様々なアプローチをいたしました。

うまくいったものもあればそうでないものもありました(ほとんどはあまり効果がありませんでした)。
概要で述べたように、エッジ保護のため、変更点の一部をダイジェスト的に記述します。
いくつかはブログにしました。

volumeを特徴量に追加

以下に記述しています。

jodawithforce.hatenablog.com


時間をカテゴリ変数に  

以下に記述しています。
どの国のマーケットがオープンしているか(サマータイム考慮)であったり、
時間アノマリーに着目して、何時であるかをカテゴリ変数にしました。
jodawithforce.hatenablog.com


talibを用いてテクニカル指標を追加  

より多数のテクニカル指標を特徴量として利用しました。
テクニカル指標のパラメーター変更
指標の期間等パラメーターを調整しました。

特徴量選択  

得られた特徴量重要度をソートし、特徴量選択を行いました。
モデルにやdfに変更を加えるたびに特徴量選択を行いました。

マルチパラメータ化  

同一指標の異なるパラメーターによって発生する現象があります。
例えば、異なる期間の移動平均線(短期、長期)をりようしたゴールデンクロスデッドクロスなどです。
これらを特徴量とするため、同一テクニカル指標を複数パラメーターで作成しました。

分類モデルを利用  

予測結果が正であれば注文を、負であれば注文をしない、というモデルの利用をします。
回帰で学習を行い、情報の損失が少ないかと考えていました。
しかしながら、最初から分類を行ったほうが、 良い結果が得られる可能性もあると思い、分類モデルで学習をしました。

他の通貨データを特徴量として追加  

通貨同士の相関はよく話題になります。
特徴量として利用できないかやってみました。

アンダーサンプリング 

分類モデルを利用する際、データの偏りが問題になる可能性があったので、
アンダーサンプリングによって調整を行いました。

学習データの利用法の変更 

KFoldによるFoldごとのモデルでバギングするのではなく、
学習データ全体を学習させ、テストデータを予測しました。

optunaによるハイパーパラメーター最適化  

optunaをつかってハイパーパラメーター最適化を行いました。

標準化  

コードには標準化を記載していますが、はじめは標準化を行っていませんでした。







モデル改善前後の比較

モデル改善のため、上記ような思いつくことを試していきました。
納得いくレベルまでモデルが改善いたしましたので、結果を記載していきます。
得られた結果はすべてテストデータでのものとなります。

回帰の評価指標の比較

得られたモデルの評価指標を比較していきます。

true_vs_predの散布図

下記の図11、12をご覧いただきますと、回帰直線がy=xに近づいたように見えます。

f:id:jodawithforce:20220107163728p:plain
図11 : 買い true_vs_pred 散布図 比較


f:id:jodawithforce:20220107163915p:plain
図12 : 売り true_vs_pred 散布図 比較

評価指標(数値)の比較

以下に、各評価指標の数値をまとめた表を掲載します。
買いの初期値においては、予測よりも平均値のほうが予測できている、という状況でしたが、各種数値はすべて改善しております。
買いの決定係数においては、特に顕著な改善が見られます。 (というよりも、初期値がひどすぎるだけかもしれません。)
個人的に、金融商品の価格予測は難易度がかなり高いと考えています。
こちらの結果は、自分にしては、まあ納得のできる改善であるとおもいます。

f:id:jodawithforce:20220107164702p:plain
図13 : 評価指標の数値 比較


分類の評価指標の比較

上述のように、回帰のモデルを利用していますが、実際には分類の利用をします。
テストデータにおける分類の評価を比較していきます。


混同行列 以下に、混同行列を比較した図を掲載します。
左上、右下の色が濃くなるほどよい分析ができています。
視覚的に、モデルが改善されていることがわかります。
また、dfに加えた操作によってデータの総数が変わっております。

f:id:jodawithforce:20220107170615p:plain
図14 : 買い 混同行列ヒートマップ 比較

f:id:jodawithforce:20220107170950p:plain
図15 : 売り 混同行列ヒートマップ 比較



評価指標(数値)の比較  

評価指標を比較していきます。
最も重要なのは、f1_scoreのweighted avg(表の赤いセル)です。
簡単にいってしまうと、重み付けされたf値の平均です。
買い、売り、どちらも初期値は0.5前後(ほぼランダムと同じ)でしたが、
改善後は0.6後半とスコアが伸びています。


f:id:jodawithforce:20220107203307p:plain
図16 : 買い 評価指標 比較

f:id:jodawithforce:20220107203327p:plain
図17 : 売り 評価指標 比較

累積リターン 比較

最後に、最も重要な累積リターンを比較していきます。
大幅な改善が見られます。



f:id:jodawithforce:20220107201840p:plain
図18 : 累積リターン 比較

最後に

今回、機械学習を利用して取引のリターン予測を行いました。
こちらの戦略が実際にこれだけの利益をあげられるか、というのはまた別問題です。
また、リターンが明らかに改善しすぎています。 改善しすぎて、何かしらのミスを疑っています。
問題ないはずですが、特徴量として予測するべき未来のデータが含まれてしまっているかもしれません。
しかしながら、機械学習を学ぶにおいて、特徴量の作成、選択、モデルへのアプローチなど、
大変よい経験になったと思います。
pandasの操作やモデルの改善では、やりたいことができない、というもどかしさを感じました。
執行方法を変えたり、取引所を変えたり、実装してみたり、やりたいことがたくさんあります。
何をするにしてもわからないことだらけですので、今後も勉強を続けていきたいと思います。






参考

GitHub - richmanbtc/crypto_data_fetcher

機械学習で仮想通貨取引をするチュートリアル。python - Qiita

Python:TA-Libでテクニカル分析、Plotlyでローソク足の描画 - Investment Tech Hack

LightGBM 徹底入門 – LightGBMの使い方や仕組み、XGBoostとの違いについて

lightgbm.LGBMRegressor — LightGBM 3.3.1.99 documentation

Python: pandas の DataFrame を scikit-learn で KFold するときの注意点 - CUBE SUGAR CONTAINER

【Python】Kaggleで引っ張りだこ!lightgbmの2種類の使い方!Training APIとScikit-learn API!【lightgbm】│SKJブログ

【Kaggle】imbalanced-learn を使ってアンダーサンプリングをしてみた - Qiita

LightGBMのハイパーパラメータチューニング(Optuna)をしてみた - Qiita