【Probspace: 6位】米国株式市場 将来株価予測の分析まとめ (Python)

1. はじめに

本記事では、「Probspace」さん主催の「米国株式市場 将来株価予測」に参加し、6位に入賞できたので、忘備録としてその解析手法をまとめてみました。

すべての特徴量が株価予測に役立つかは不明ですが、本コンペ以外にも役立ちそうなコードを列挙していきます。

2. コンペ概要

  • 米国証券取引所の株価推移データを用いて、1週間後の株価を予測するアルゴリズムを開発する
  • 2011/11/13~2019/11/17週の計419週間にわたる米国株データ(終値)をもとに、2019/11/24週の終値を予測する
  • 外部データは使用不可である
  • モデルの予測性能は評価関数RMSLE(Root Mean Squared Logarithmic Error)で評価する
  • 上位者は、コンペ終了後ソースコードを公開する必要がある(株価を調査して提出するなどの不正を防ぐため)

3. 解析環境

  • Google Colab Pro
  • Memory: 37.8GB

4. 下準備コード

4-0. ライブラリインポート

import os
import gc
import glob
import warnings
warnings.filterwarnings('ignore')
import re
import time
import pandas as pd
pd.set_option('display.max_columns', 500)
pd.set_option('display.max_rows', 100)
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns

from tqdm import tqdm

import scipy.stats as stats
import statsmodels.api as sm

from statsmodels.tsa import stattools
from statsmodels.tsa.api import VAR
from sklearn.metrics import mean_squared_error as mse
from sklearn.cluster import KMeans
from sklearn.preprocessing import RobustScaler, LabelEncoder, MinMaxScaler
from sklearn.model_selection import KFold

import lightgbm as lgb

4-1. Google Colab上でのMy Drive接続

まずは、コンペデータをGoogle Colab上に格納して、データ分析に使用できるようにMy Driveと接続できるようにします。コードは以下です。

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

4-2. Debugとデータ読込

Debug用class、変数処理に便利です。

class Debug:
    def __init__(self):
        self._debug = 0
        self.today = "20211224"
        self.revision = 1
        self.data_dir = './drive/MyDrive/probspace/input/'
        self.output_dir = './drive/MyDrive/probspace/output/'
        self.n_clusters = 11
        os.makedirs(self.output_dir + self.today, exist_ok=True)
        
debug = Debug()

def read_train():
    train = pd.read_csv(debug.data_dir + 'train_data.csv', encoding='cp932')
    train.iloc[-1, 0] = "2019/11/24"
    train["Date"] = pd.to_datetime(train["Date"])
    return train

train_raw = read_train()
original_col_num = len(train_raw.columns) - 1

4-3. 株価のトレンドとQ-Qプロット確認

株価のヒストグラムとQ-Qプロット(正規分布に近いかどうかを目視する)を表示する関数を作成しました。すべての株式を同時に表示すると処理が重いので、List型のcolsを引数にしています。

def hist_prob_show(df, cols):
    df[cols].plot(subplots=True, kind='hist', bins=100, figsize=(25,15))
    plt.show()

    fig, ax = plt.subplots(1, 5, figsize=(25, 6))
    for i, col in enumerate(cols):
        stats.probplot(df[col], dist='norm', plot=ax[i])
        ax[i].set_title(col)
    plt.show()

def log_diff_return(df):
    df = df.dropna().set_index('Date')
    return np.log1p(df).diff(1).fillna(0)

cols = ["KZIA",	"DAC", "SSI",	"VHC", "MTEM"]
train_df_diff = log_diff_return(train_raw)
hist_prob_show(train_df_diff, cols)

解析結果では、株式のex.log1pとdiff(1)を組み合わせると正規分布に近くなるため、前日との株価の変化率を目的関数として設定しました(図1参照ください)

図1 解析結果の一部

4-4. 標準化関連

株価の変化率の標準化を行う関数です。基本的には、株価の変化率は最大最小値がないため、通常の標準化を実施しました。

train_df_diff_norm = (train_df_diff - train_df_diff.mean()) / train_df_diff.std(ddof=0)

scaler = RobustScaler()
train_df_diff_robust = scaler.fit_transform(train_df_diff)
train_df_diff_robust = pd.DataFrame(train_df_diff_robust, columns=train_df_diff.columns)

scaler = MinMaxScaler(feature_range=(0, 1))
train_df_diff_maxmin = scaler.fit_transform(train_df_diff.values)
train_df_diff_maxmin = pd.DataFrame(train_df_diff_maxmin, columns=train_df_diff.columns)

4-5. Companyデータ解析

Comapanyデータにある会社名から、特徴的な文字列の出現個数を検討しています。図2のようにある程度、同じような文字列を有する会社名が存在するようです。

def read_compnay():
    company = pd.read_csv(debug.data_dir + 'company_list.csv').rename(columns={"Symbol": "id"})
    return company

def word_show(df, top=21):
    sns.set(font_scale=1.5)
    gridspec_kw = dict(width_ratios=[2, 1], height_ratios=[1])
    fig, ax = plt.subplots(1, 2, figsize=(30, 6), gridspec_kw=gridspec_kw)
    sns.barplot(data=df.iloc[0:top], x="Word", y="Count", ax=ax[0])
    ax[0].set_xticklabels(df["Word"].tolist(), rotation=90)
    ax[1].pie(df["Count"].iloc[0:top], counterclock=False, startangle=90, wedgeprops={'linewidth': 1, 'edgecolor':"white"}, labels=df["Word"].iloc[0:top], labeldistance=1.1, textprops={'color': "black"})
    plt.show()

def search_keyword(df, kw):
    tmp_df = df[df["name_clean"].str.contains(kw)]
    n_stocks = len(tmp_df)
    print(f"Target: {kw}, {n_stocks}件")
    tmp_df[kw] = 1
    return tmp_df[["id", kw]]

company = read_compnay()
company['name_clean'] = company["Name"].apply(lambda x: re.sub(r'[^(a-zA-Z\s)]', '', x).lower().strip())
txt_lst = company["name_clean"].values.tolist()
txt_lst = [x.split() for x in txt_lst]
company['text_list'] = txt_lst

word_dict = {}
word_lst = []
for i in range(len(company['text_list'])):
    for word in company['text_list'][i]:
        if word not in word_lst:
            word_lst.append(word)
            word_dict[word] = 1
        else:
            word_dict[word] += 1
word_dict = sorted(word_dict.items(), key=lambda x: x[1], reverse=True)
df_word = pd.DataFrame(word_dict, columns=["Word", "Count"])
word_show(df_word, top=25)

for kw in df_word["Word"].tolist()[0:25]:
    search_df = search_keyword(company, kw)
    company = pd.merge(company, search_df, how="left", on="id")
    company[kw] = company[kw].fillna(0)
    
company["corporate"] = company["corporation"] + company["corp"]
company["corporate"] = company["corporate"].apply(lambda x: x//2 if x==2 else x)
company = company.drop(["corporation", "corp"], axis=1)
図2 Company文字列の解析結果

4-6. RMSE-Cluster解析

株価の変化率のトレンドが似ている株式をクラスタリングで取得します。ただし、”mse_result”を得るには、今回のデータで1時間くらい必要なため、二度目以降の解析では、csvデータに保存したデータを使用するようにしました。

def read_data(filepath):
    df = pd.read_csv(debug.output_dir + filepath).iloc[:, 1:]
    df.index = train_raw.columns[1:]
    df.columns = train_raw.columns[1:]
    return df

def cluster_df(df):
    km = KMeans(n_clusters=debug.n_clusters, init='k-means++', n_init=10, max_iter=300, tol=1e-04, random_state=2021).fit(df.values)
    df_cls = pd.DataFrame(km.labels_, index=df.index, columns=['class'])
    return df_cls

df = log_diff_return(train_raw)

"""
2回目以降はcsvファイルを読み込む
mse_result = [[mse(df.iloc[:, i], df.iloc[:, j], squared=False) for j in range(original_col_num)] for i in range(original_col_num)]
"""

pd.DataFrame(mse_result).to_csv(debug.output_dir + "rmse_result.csv")
df_rmse = read_data("rmse_result.csv")

df_cls = cluster_df(df_rmse).reset_index(drop=False).rename(columns={"index": "id", "class": "class_rmse"})

4-7. VAR解析

VAR関数を利用した解析です。各株式でどの程度のLagが良いのかを参考にしました。

class VARS:
    def create_adf_df(self, data):
        adf_lst = []
        for name in data.columns:
            adf_dict = {}
            adf = stattools.adfuller(data[name].values, regression='ctt')
            adf_dict['id'] = name
            adf_dict['p_value_%'] = adf[1] * 100
            adf_dict['adf_lag'] = adf[2]
            adf_lst.append(adf_dict)
        return pd.DataFrame(adf_lst)
   

var_inst = VARS()
train_var = log_diff_return(train_raw)
adf_df = var_inst.create_adf_df(train_var)
adf_df["adf_lag"].value_counts()

5. 分析コード

5-1. Preprocessing

“create_diff_norm”関数では、株価データが、各株式が列ごとに配置されているため、melt関数で縦積み再配置します。その後、統計処理などで目的関数や説明変数を作成します。

def create_diff_norm(df):
    df = log_return(df)
    df.loc[pd.to_datetime('2019-11-24'), :] = np.nan
    train_date = pd.Series(df.index)
    df.reset_index(drop=True, inplace=True)

    week_len = len(df)
    df = df.T.reset_index().rename(columns={'index': 'id'})
    df = pd.melt(df,
                id_vars='id',
                value_vars=[week for week in range(week_len)],
                var_name='Week',
                value_name='y'
                )
    df['Year'] = df['Week'].map(train_date).dt.year
    df['Month'] = df['Week'].map(train_date).dt.month
    df['Day'] = df['Week'].map(train_date).dt.day
    df['WeekOfYear'] = df['Week'].map(train_date).dt.isocalendar().week.astype(int)

    df["y_mean"] = df[['id', "y"]].groupby(['id'])["y"].transform(lambda x: x.mean())
    df["y_std"] = df[['id', "y"]].groupby(['id'])["y"].transform(lambda x: x.std(ddof=0))
    df["y_norm"] = (df["y"] - df["y_mean"]) / df["y_std"]
    df["y_norm_all"] = (df["y"] - df["y"].mean()) / df["y"].std(ddof=0)

    df["y_diff"] =  df[['id', "y"]].groupby(['id'])["y"].transform(lambda x: x.diff(1))
    df["y_diff_norm_all"] = (df["y_diff"] - df["y_diff"].mean()) / df["y_diff"].std(ddof=0)

    val_col = "y_diff"
    df['y_diff_max'] = df[["id", val_col]].groupby(['id'])[val_col].transform(lambda x: x.max())
    df['y_diff_min'] = df[["id", val_col]].groupby(['id'])[val_col].transform(lambda x: x.min())
    df['y_diff_median'] = df[["id", val_col]].groupby(['id'])[val_col].transform(lambda x: x.median())
    df['y_diff_0.75'] = df[["id", val_col]].groupby(['id'])[val_col].transform(lambda x: x.quantile(0.75))
    df['y_diff_0.25'] = df[["id", val_col]].groupby(['id'])[val_col].transform(lambda x: x.quantile(0.25))
    df['y_diff_kurt'] = df[["id", val_col]].groupby(['id'])[val_col].transform(lambda x: x.kurt())
    df['y_diff_skew'] = df[["id", val_col]].groupby(['id'])[val_col].transform(lambda x: x.skew())
    return df

train = create_diff_norm(train_raw)

“preprocessing”では、各株式のLagや移動平均値、ラベルエンコーディングなどを実施します。

def preprocessing(df, target="y_diff_norm_all", lag_len=19):
    def create_lags(df, target, lags):
        print("Create lags")
        lag_df = pd.DataFrame()
        for lag in tqdm(lags):
            lag_df[f'y_diff_lag{lag}'] = df[['id', target]].groupby(['id'])[target].transform(lambda x: x.shift(lag))
            if lag == 1:
                lag_df[f'y_diff_lag_total_{lag}'] = lag_df[f'y_diff_lag{lag}']
            else:
                lag_df[f'y_diff_lag_total_{lag}'] = lag_df[f'y_diff_lag{lag}'] + lag_df[f'y_diff_lag_total_{lag-1}']
        lag_df = lag_df.drop(['y_diff_lag_total_1'], axis=1)
        return lag_df

    def create_rolls(df, target, lags, rolls):
        def sma(series, roll, lag=1):
            return series.shift(lag).rolling(window=roll, min_periods=1).mean()

        print("Create Rollings")
        roll_df = pd.DataFrame()
        for roll in tqdm(rolls):
            for lag in lags:
                roll_df[f'y_diff_lag{lag}_sma_mean_roll{roll}'] = df[['id', target]].groupby(['id'])[target].transform(sma, roll, lag)
         return roll_df

    df = pd.merge(df, company, on='id', how='left')
    df['enc_Sector'] = LabelEncoder().fit_transform(df['Sector'])
    df['enc_Industry'] = LabelEncoder().fit_transform(df['Industry'])
    df['enc_List'] = LabelEncoder().fit_transform(df['List'])

    df = pd.merge(df, df_cls[["id", "class_rmse"]], how="left", on="id")
    df = pd.merge(df, adf_df[["id", "adf_lag"]], how="left", on="id")

    lags = list(np.arange(1, lag_len, 1))
    lag_df = create_lags(df, target, lags)
    df = pd.concat([df, lag_df], axis=1)

    roll_df = create_rolls(df, target, [1], [4, 9, 13, 26, 52])
    df = pd.concat([df, roll_df], axis=1)

    return df

train = preprocessing(train, target="y_diff_norm_std", lag_len=19)

5-2. LGBM class

LightGBM用のClassを作成。また、各Hold中のRMSLEを算出するようにしました。

class LGBM:
    def __init__(self, feature_column, train, target, df_std):
        self.num_boost_round = 5000
        self.verbose_eval = 1000
        self.n_splits = 5
        self.early_stopping_rounds = 100
        self.kf = KFold(n_splits=self.n_splits, shuffle=True, random_state=2021)
        self.y_vals = []
        self.y_preds = []
        self.boosters = []
        self.eval_results_lst = []
        self.train, self.target, self.df_std = train, target, df_std
               
    def lgbm_train(self, params, reg):
        y_oof = np.zeros(len(self.target))
        for fold, (trn_idx, val_idx) in enumerate(self.kf.split(self.train, self.target)):
            eval_results = {}
            print("="*15 + f' Fold {fold+1} started at {time.ctime()} ' + "="*15)
            lgb_trn = reg.Dataset(self.train[trn_idx], label=self.target[trn_idx], categorical_feature='auto')
            lgb_val = reg.Dataset(self.train[val_idx], label=self.target[val_idx], categorical_feature='auto')
                        
            booster = reg.train(params=params, 
                                train_set=lgb_trn, 
                                valid_sets=[lgb_trn, lgb_val],
                                valid_names=['Train', 'Valid'], 
                                num_boost_round=self.num_boost_round, 
                                early_stopping_rounds=self.early_stopping_rounds, 
                                verbose_eval=self.verbose_eval,
                                evals_result=eval_results,
                                categorical_feature='auto')


            self.eval_results_lst.append(eval_results)
            self.boosters.append(booster)
            self.y_vals.append(booster.predict(self.train[-X_test.shape[0]:], num_iteration=booster.best_iteration))
            self.y_preds.append(booster.predict(X_test.values, num_iteration=booster.best_iteration))
            y_oof[val_idx] = booster.predict(self.train[val_idx], num_iteration=booster.best_iteration)

            print(
                 f'Fold {fold+1} RMSLE:',
                 np.sqrt(np.mean(np.square((y_oof[val_idx] - self.target[val_idx]) * self.df_std[val_idx]))),
                 "\n"
               )

        print(
             'Training RMSLE:',
             np.sqrt(np.mean(np.square((y_oof - self.target) * self.df_std)))
             )

5-3. LGBM train

LGBM実行用のコードです。ハイパーパラメターチューニングはoptunaを用いて実行済みです。

start_week = 1
X_test = train[(train['Week'] == 419)].reset_index(drop=True)
X_train = train[(train_tmp['Week'] < 419) & (train_tmp['Week'] >= start_week)].reset_index(drop=True)
y_train = train[(train_tmp['Week'] < 419) & (train_tmp['Week'] >= start_week)].reset_index(drop=True)["y_diff_norm_all"]

print(f"X_train: {X_train.shape}")
print(f"y_train: {y_train.shape}")
print(f"X_test: {X_test.shape}")

lgbm_inst = LGBM(X_train.columns,
          X_train.values,
                 y_train.values,
                 train["y_diff"].std(ddof=0)
                 )

best_params = {
               'objective': 'regression', 
               'metric': 'rmse', 
               'boosting_type': 'gbdt', 
               'verbose': -1,
               'random_state': 2021
              }
best_params.update({
                    'feature_pre_filter': False, 
                    'lambda_l1': 0.0, 
                    'lambda_l2': 0.0, 
                    'num_leaves': 255, 
                    'feature_fraction': 0.82, 
                    'bagging_fraction': 1.0, 
                    'bagging_freq': 0, 
                    'min_child_samples': 25
                    })
lgbm_inst.lgbm_train(best_params, lgb)

5-4. 学習経過の確認

Validに対してTrainが過学習になっていないか、グラフで評価値の傾向を確認します。

scores = list()
for result in lgbm_inst.eval_results_lst:
    scores.append(result['Valid']['rmse'][-1])

5-5. 予測値と実値の差異を確認

予測値と実値の差異を確認し、必要な説明変数を考察します。ここが一番難しいですね。株価は基本的にはあまり変動しませんが、たまに大きく突然変動することがあるため、その傾向をうままくとらえることが困難でした。

train_tmp = train[["id", "y", "y_diff"]].iloc[-X_test.shape[0] * 2:-X_test.shape[0]]

y_pred = np.mean(lgbm_inst.y_vals, axis=0) * train["y_diff"].std(ddof=0) + train["y_diff"].mean()
train_pred_lgbm = pd.DataFrame(y_pred, index=train_tmp["id"].unique())
train_pred_lgbm.columns = ['y']
train_pred_lgbm.index.name = 'id'
train_pred_lgbm['y_real'] = train_tmp["y_diff"].values

train_pred_lgbm["RMSLE"] = rmse(train_pred_lgbm["y"].values, train_pred_lgbm["y_real"].values)
print(f"418 RMSLE: {np.sqrt(train_pred_lgbm['RMSLE'].mean())}")
print(f"Worst10: {list(train_pred_lgbm.sort_values('RMSLE', ascending=False)['id'][:10])}")

5-6. 提出用csv作成

5-3で実行済みの”lgbm_inst.y_preds”から、提出用のcsvファイルを作成するコードです。

train_tmp = train[["id", "y", "y_diff"]].iloc[-X_test.shape[0] * 2:-X_test.shape[0]]

y_pred = np.mean(lgbm_inst.y_preds, axis=0) * train["y_diff"].std(ddof=0) + train["y_diff"].mean()
y_pred = np.expm1(y_pred + train_tmp["y"].values)
submission_lgbm = pd.DataFrame(y_pred, index=train_tmp["id"])
submission_lgbm.columns = ['y']
submission_lgbm.index.name = 'id'
submission_lgbm = submission_lgbm.where(submission_lgbm['y'] > 0, 0)
submission_lgbm.to_csv(f'{debug.output_dir}{debug.today}/lgbm_{debug.today}_all.csv')

少し説明足らずの個所もありますが、おおよその実行内容は以上となります。

6. 所感

株価予測を経験して以下を得ました。

  • 単純な時系列予測では予測が困難
  • 微小な株価変動はある程度予測できるが、変動率が大きい予測は困難。これは全株価の変動率がほぼ1に近いため、変動率が大きい値は±5σ以上の確率でしか起きないためと推察される
  • 必要な説明変数は、株式によって異なるかもしれないので試してみたい
  • lambda関数を用いたtransformとapplyの違いを理解できた

7. 参考文献、サイト

米国株式市場 将来株価予測 2位解法

LightGBM Base line(LB=0.03781)

Pythonを使ってVARモデルにおける多変量時系列予測モデルの構築

Pandas の transform と apply の基本的な違い

Baran

それでは、よいPythonライフを!

Baran

以下の記事も書いています。よかったらご覧ください。

Baran

Twitterもやってまーす。

コメントを残す

メールアドレスが公開されることはありません。 が付いている欄は必須項目です

ABOUT US

Baran-gizagiza
経歴:浪人→理系大学院卒業→大手製造業に就職(技術職)→アメリカ赴任中 仕事は、研究・設計など上流工程の仕事に携わっています。企業勤務を継続しながら、新しいことにチャレンジしたいと思い、ブログを下記はじめました。 このブログでは、趣味である 筋トレ(健康、ダイエット) AIとデータ(高校数学、プログラミング) 読書(主に自己啓発系) を中心に、人生経験やおすすめ情報の発信しています。