Kaggleに挑戦！【Jane Street Market Prediction】②

　こんにちは。システムトレーダーの卵ことKenKenです。これまでの記事では、KaggleのAPIを用いたデータ取得方法やKaggle上のNotebookからダウンロード無しでデータセットを利用する方法について書きました。
　今回は、”Jane Street Market Prediction”のデータ特性について簡単に見ていきたいと思います。最近、知ったんですかこういう作業のことをEDA(探索的データ解析)と言うそうです。

※本記事は、公開しているNotebookをまとめたものとなっております。

データの読込みとメモリの節約

import numpy as np # linear algebra
import pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)
pd.set_option('display.max_columns', 140) #最大表示列数の指定
import matplotlib.pyplot as plt

# dfの各列の型を設定しメモリ軽減
def reduce_mem_usage(df, verbose=True):
    numerics = ['int16', 'int32', 'int64', 'float16', 'float32', 'float64']
    start_mem = df.memory_usage().sum() / 1024**2
    for col in df.columns:
        col_type = df[col].dtypes
        if col_type in numerics:
            c_min = df[col].min()
            c_max = df[col].max()
            if str(col_type)[:3] == 'int':
                if c_min > np.iinfo(np.int8).min and c_max < np.iinfo(np.int8).max:
                    df[col] = df[col].astype(np.int8)
                elif c_min > np.iinfo(np.int16).min and c_max < np.iinfo(np.int16).max:
                    df[col] = df[col].astype(np.int16)
                elif c_min > np.iinfo(np.int32).min and c_max < np.iinfo(np.int32).max:
                    df[col] = df[col].astype(np.int32)
                elif c_min > np.iinfo(np.int64).min and c_max < np.iinfo(np.int64).max:
                    df[col] = df[col].astype(np.int64)
            else:
                if c_min > np.finfo(np.float16).min and c_max < np.finfo(np.float16).max:
                    df[col] = df[col].astype(np.float16)
                elif c_min > np.finfo(np.float32).min and c_max < np.finfo(np.float32).max:
                    df[col] = df[col].astype(np.float32)
                else:
                    df[col] = df[col].astype(np.float64)
    end_mem = df.memory_usage().sum() / 1024**2
    if verbose: print('Mem. usage decreased to {:5.2f} Mb ({:.1f}% reduction)'.format(end_mem, 100 * (start_mem - end_mem) / start_mem))
    return df

csvから読み込んだデータは、データ型を指定しないとint64, float64になります。大量のデータを扱う場合は、それが原因でメモリ不足エラーになります。それを知らずに、今回は大変苦労(イライラ)しました。。。上記の関数を使用することでそれぞれのデータ列に対して最適なデータ型を設定することができます。

# before
train = pd.read_csv('/kaggle/input/jane-street-market-prediction/train.csv')
train.info()

# after
train = reduce_mem_usage(train)
train.info()

# dfをpickleで保存
train.to_pickle('/kaggle/working/JSMP_train.pickle')

　メモリ使用量を2.5GB⇒631.49MB(74.9%減)と大幅に節約することができました。csvからデータを読み込むとかなり時間がかかるので、pickleで保存しておきます。

各特徴量の欠損値について

　分析を始める前に、各特徴量の欠損値の数について見てみる。

null_num_train = {feature:num for feature, num in zip(train.columns, train.isnull().sum())}
print('max null ratio:', max(null_num_train.values()) /train.index.size)
print('null num:{}'.format(len([i for i in null_num_train.values() if i != 0])))

plt.figure()
plt.bar(null_num_train.keys(), null_num_train.values())

　欠損値は、88個の特徴量に含まれている。欠損値の最大数は、395,535個でfeature_17,feature_18, feature_27, feature_28の4つが欠損値が最も多い。トレーニングデータの各変数は欠損値を含めて2,390,491個あり、欠損値は各変数において多くても約16.5％である。

　テストデータの欠損値について見てみる。同様の手順で欠損値について見たら、テストデータにも80個の特徴量に欠損値が含まれているようである(最大で1変数内に約13.5%)。

　学習データ、テストデータにおいて欠損値を持つ特徴量について確認したところ同様であった。テストデータにおいても欠損値が多いことから、欠損値の補完をどう対処するかが問題となる。テストデータに欠損値が含まれていなければ、トレーニングデータで欠損値を含むものを削除するだけでよかったがそう簡単にはいかないか。。。

欠損値の補完

　欠損値はどのような条件で発生しているのか見てみる。ある時点まではデータが取得できないために欠損値が発生している場合やデータ全体の中で点々と欠損値がある場合など考えられる。前者であれば、ある時点までのデータをごっそり削るのもありだと思われる。後者の場合は、平均値や前の値で埋めるなどいろいろ考える必要がある。

　欠損値が最も多かった4変数のうち、”feature_17″について欠損値の推移を簡単に見てみる。データ量が多い為、今回は最初と最後の10％に絞った。

# 欠損値の推移(欠損値が最も多いfeature_17)をプロットしてみる
# 欠損値:1 else:0として全体の最初から10％と最後の10%
train['feature_17'].map(lambda x:1 if x != x else 0)[:int(train.index.size*0.1)].plot()
plt.legend(loc='lower right')
plt.xlabel('0~10%')

# 90%~
train['feature_17'].map(lambda x:1 if x != x else 0)[int(train.index.size*0.9):].plot()
plt.legend(loc='lower right')
plt.xlabel('90~100%')

　欠損値は、期間を通じてまばらのようだ。欠損値の補完は、単純にゼロ埋めや前後の値やその平均で補完する、全期間の平均、一定期間の平均を用いるなどやり方はいろいろある。今回のデータは時系列で期間幅も短いので前の値で埋め合わせをすることにする。金融データということで、各指数の公表タイミングが異なり欠損値となってしまっている可能性がある。つまりある時点のデータが欠損値の場合、そのデータは公表前であり、市場では公表前の前回の値が参照されているということはよくあることだと思われる。

train = pd.read_csv('/kaggle/input/jane-street-market-prediction/train.csv')
# 欠損値の補完(前の値で補完する)
train.fillna(method = 'ffill', inplace=True) 
train.dropna(inplace=True)
# メモリ対策
train = reduce_mem_usage(train)
train.to_pickle('/kaggle/working/train_without_null.pickle')

fillna()を使う際は、dfの型が同じでなければならないらしい。これまでのままfillna()を使うと「TypeError: No matching signature found」が出てしまう。そのため、再度、train.csvを読み込む処理を織り込んでいる。

前の値で欠損値を補完しているため、欠損値でない値が出るまでは、全て欠損値となっている。補完後の各変数の欠損値の数は、多くて479個であり、欠損値を含む場合は削除することにした。（つまり最初から479レコードまでは削除される）

　ここまでで、最低限の前処理は完了した。テストデータの欠損値処理は、同様に前の値で補完することにする。したがって、テストデータにおいては最初の数レコードは評価対象外とする必要が出てくる。

データ全体について

　トレーニングデータについて簡単にプロットしてみて特徴を見ていく。

# リターンの分布と特徴
# 全部, weight=0, weight=0以外でヒストグラムをみてみる
plt.figure()
train['resp'].hist(bins=25, label='all_data', color='lightblue')
train.loc[train['weight']==0]['resp'].hist(bins=25, label='weight=0')
train.loc[train['weight']!=0]['resp'].hist(bins=25, label='weight!=0', histtype='step')
plt.legend()

まずは、リターンの分布についてヒストグラムにしてみた。weight=0の数が少ないことがわかる(全体の17%程度)。weight=0は、ボラティリティがweight!=0に比べて小さいこともわかる。分布の形状は、weight=0はややプラスより、weight!=0はややマイナスよりとなっている。それぞれの統計量を見てみる。

# 統計量を算出
import statsmodels.stats.api as sms
print('weight=0')
print(train.loc[train['weight']==0]['resp'].describe())
JB = sms.jarque_bera(train.loc[train['weight']==0]['resp']) # JB検定
print('歪度:{0:.4f}, 尖度:{1:.4f}'.format(JB[2], JB[3]))

print('\n weight!=0')
print(train.loc[train['weight']!=0]['resp'].describe())
JB = sms.jarque_bera(train.loc[train['weight']!=0]['resp']) # JB検定
print('歪度:{0:.4f}, 尖度:{1:.4f}'.format(JB[2], JB[3]))

　ヒストグラムを見た時とは若干結果が異なった。weight=0は平均は若干プラスだが、歪度はマイナス（マイナス寄り）となっている。一方、weight!=0は平均は若干プラスで歪度はプラス（プラスより）である。やはり、見るだけでなくちゃんと数値として確認することは大事。

weightとresp(リターン)の関係

　ここまでの結果を踏まえると、weightとリターンには何かしらの特徴がありそう。これらを2変数として散布図を書いてみる。

# weightとリターンの関係
plt.scatter(train['weight'], train['resp'])
plt.xlabel('weight')
plt.ylabel('resp')

なかなか面白い散布図が描かれた。weightが大きくなるに従い、respのばらつきは小さくなっている。weight=0では、高いリターンは得られる可能性はあるものの、同時に大きな損失を被る可能性もある。weightはリスク管理として機能しそうである。ただ、weightを大きくするとリターンも得られないのでその辺の調整が難しい。逆に、weightがある一定値より大きい場合はリターンはほぼゼロなので取引をしないという戦略を加えると良いのかも。そうなると狙いとしては、0<weight<75くらいでresp>0の部分をうまく取り出すことか？

resp(リターン)と各特徴量(feature_X)の関係

　次に、respと130個の特徴量の相関を見てみるために、まずは素直にヒートマップとして表示してみる。

# 相関ヒートマップ
import seaborn as sns
sns.heatmap(train.drop(['date', 'resp_1', 'resp_2', 'resp_3', 'resp_4', 'ts_id'],axis=1).corr())

　2行目がrespだが軸から表示が消えてしまっている。すべての値を表示するには、xticklabels, yticklabelsを指定することで可能だが、変数が多すぎるので結局読めなかった（グラフサイズを大きくすれば可能かもしれないが）。2行目を大雑把に見てみると、ところどころで正・負の相関がある特徴量が見られる。また、変数間で相関があるものがいくつも見られる。したがって、次元圧縮してみるのもよいかもしれない。次元圧縮としては、主成分分析やオートエンコーダ、クラスタリングあたりを試したい。

　ここで、respと相関が高い10変数を抽出してみる。

# 相関が高いものを抽出
cols = train.drop(['date', 'resp_1', 'resp_2', 'resp_3', 'resp_4', 'ts_id'],axis=1).corr().nlargest(10, 'resp')['resp'].index
cm = np.corrcoef(train[cols].values.T)

plt.subplots(figsize = (12,10))
sns.heatmap(cm, vmax=.8, linewidths=0.01, annot=True, cmap='viridis',
            xticklabels=cols.values, yticklabels=cols.values, annot_kws={'size':14});

　相関は高くても0.046であり、どの変数もrespとはほぼ無相関である。相関が低い理由としては、全期間のデータを対象としていることや外れ値の存在などが考えられる。外れ値については、最初に見ておくべきだった。。。

　期間を短くしたときの相関係数について見てみる。今回は、前半10%について見てみた。

cols = train.drop(['date', 'resp_1', 'resp_2', 'resp_3', 'resp_4', 'ts_id'],axis=1)[:int(len(train.index)*0.1)].corr().nlargest(10, 'resp')['resp'].index
cm = np.corrcoef(train[cols].values.T)

plt.subplots(figsize = (12,10))
sns.heatmap(cm, vmax=.8, linewidths=0.01, annot=True, cmap='viridis',
            xticklabels=cols.values, yticklabels=cols.values, annot_kws={'size':14});

　respとの相関について見てみると、抽出された変数は多少異なるものの相関係数のサイズに変化は見られなかった。

外れ値について確認

　外れ値をむやみに欠損値にしたりある値で丸めるのもどうなのだろうか。外れ値だからこそ特徴を持つということはないだろうか。そう考えると、外れ値処理の有無についても見ていく必要がありそうだ。

plt.figure(figsize=(20, 100))

cnt = 0
x = train.drop(['date', 'resp_1', 'resp_2', 'resp_3', 'resp_4', 'ts_id'],axis=1).copy()
for i in range(29):
    for j in range(4):
        plt.subplot(29, 4, cnt+1)
        sns.boxplot(x=x.iloc[:, cnt])
        plt.xlabel(x.columns[cnt], fontsize=10)
        gc.collect() # メモリ解放
        cnt += 1

四分位に基づくと多くの変数で外れ値が散見された。また、全体としては左右対称か正方向に偏った形状をしていた。ここで、四分位範囲に基づく外れ値判定（参考サイト）を行い、外れ値の有無がresp(リターン）に影響があるか調べてみる。

x = train.drop(['date', 'weight', 'resp_1', 'resp_2', 'resp_3', 'resp_4', 'ts_id'],axis=1).copy()
corr_outlier = []
corr_inner = []

for col_x in x.columns[1:]:
    df_x = x[['resp', col_x]]
    Q1 = df_x[col_x].quantile(.25)
    Q3 = df_x[col_x].quantile(.75)
    IQR = Q3 - Q1
    th1 = Q1 - 1.5 * IQR
    th2 = Q3 + 1.5 * IQR
    df_outlier = df_x[(df_x[col_x] < th1) | (df_x[col_x] > th2)]
    df_inner = df_x[(df_x[col_x] > th1) &amp; (df_x[col_x] < th2)]
    corr_outlier.append(df_outlier.corr().loc['resp',col_x])
    corr_inner.append(df_inner.corr().loc['resp',col_x])
    gc.collect()

plt.figure(figsize=(12, 10))
plt.bar(x.columns[1:], corr_outlier, color='red', label='outlier')
plt.bar(x.columns[1:], corr_inner, color='lightblue', label='inner')
plt.legend()