回帰モデルでファンダメンタルデータから理論株価を分析する

売上や配当などのファンダメンタルズデータから理論株価を分析します。
今回はSBI証券から手に入るデータを例にとり、回帰モデルとしてニューラルネットワークを作成し理論株価を計算します。
分析にはpythonを使いました。

データの入手

売上高などのファンダメンタルデータを説明変数にし、株価を目的変数にするため、それらのデータを入手する必要があります。
今回はSBI証券の国内株式のスクリーニング機能からファンダメンタルデータをダウンロードします。スクリーニング機能の利用にはSBI証券へのログインが必要です。

具体的な手順としてはまず、スクリーニング機能の絞り込み項目から必要な項目を選択します。この機能では選択した項目が全て揃っている銘柄のみが抽出されます。そのため、たくさんの項目を選択すると抽出される銘柄が極々小数になってしまいます。

検索後は結果の右上に表示される「CSVダウンロード」を押すことで検索結果がダウンロードされます。

ダウンロードされたファイル"screener_result.csv"を適当な場所に移動すればデータの準備は完了です。今回はソースコードと同じ場所に保存しました。
本記事では2024/2/10(土)にダウンロードしたデータでの分析結果を示します。

データの読み込み

pythonを用いたデータ分析に移ります。
まずはSBIからダウンロードしたスクリーニングの結果を読み込みます。

import pandas as pd
import numpy as np

screener_result = pd.read_csv('screener_result.csv')
screener_result

特徴量エンジニアリング

読み込んだデータを回帰モデルで学習しやすいデータに成形していきます。

# 特徴量エンジニアリング用にコピー
data = screener_result.copy()

# 数値データ列のリネームと数値化
rename_map = {'現在値':'株価',
              '売上高(百万円)':'売上高',
              '経常利益(百万円)':'経常利益',
              '当期純利益(百万円)':'当期純利益',
              'EPS(1株当たり当期利益)(円)':'EPS',
              'BPS(1株当たり純資産)(円)':'純資産',
              '配当利回り(%)':'配当利回り'}
## リネーム
data = data.rename(columns = rename_map)

読み込んだデータの型を確認します。

ほとんどの列が文字列になってしまっているので数値に変換します。文字列中の桁区切りのカンマも消します。

## 数値化
def convert_to_numeric_with_comma(col):
    return pd.to_numeric(col.astype(str).str.replace(',', ''), errors='coerce')

data[list(rename_map.values())] = data[list(rename_map.values())].apply(convert_to_numeric_with_comma, axis=0)

目的変数が「株価=一株あたりの時価総額」なので売り上げ高などのファンダメンタルデータを株式数で除して一株あたりの値に変換します。
SBIのスクリーニングでは株式数は得られないため、純利益をEPSで除して株式数を計算する必要があります。

# 一株あたりの数値に変換
data['株式数'] = data['当期純利益'].div(data['EPS'], axis=0)*1000000
data[['売上高', '経常利益', '当期純利益']] = data[['売上高', '経常利益', '当期純利益']].div(data['株式数'], axis=0)*1000000
data['配当'] = data['配当利回り']*data['株価']/100

data = data[['コード', '銘柄名', '株価', '売上高', '経常利益', '当期純利益', '純資産', '配当']]

# 欠損値をdrop
data = data.dropna(axis=0, how='any')

data = data.sort_values('コード').reset_index(drop=True)

特徴量の分布を確認します。プロットするコードは省略します

金額のデータにありがちな裾の長い分布になっていますね。このままだと学習時の勾配が過大・過少になり学習しにくいので対数変換を行います。符号付のデータなので、符号付き対数変換を適用します。

cols = ['株価', '売上高', '経常利益', '当期純利益', '純資産', '配当'] #'株価', 

# 符号付き対数に変換
data[cols] = np.sign(data[cols]) * np.log1p(np.abs(data[cols]))

（複数の）正規分布のような分布に変換できました。
御覧の通り、今回は目的変数である株価にも対数変換を適用しました。そのため、この後行うモデル学習時においてモデルの損失関数の設計を工夫する必要があります。

最後にデータを説明変数と目的変数に分けます。

X = data[['売上高', '経常利益', '当期純利益', '純資産', '配当']].values
y = data[['コード', '株価']].values

目的変数としては株価だけが必要ですが、out-of-foldな学習・予測をした後の結果の取り扱いを楽にするため、対応する銘柄コード情報も残します。

ニューラルネットワークモデルによる回帰分析

今回はkeras、tensorflowを用いてニューラルネットワークによる回帰モデルを作成します。回帰モデルには色々な選択肢がありますが、説明変数と目的変数の間の緩い単調性と非線形な関係性を表現する上でニューラルネットワークが適当だと判断しています。データがわかりやすいテーブル形式なので回帰ブースティング木も思いつきますが、単調性の表現が難しいため私は使いこなせていません。

モデルの定義

まずはモデルの定義を行います。
学習時に何度かモデルを作成するため関数化しておきます。

from keras.models import Model
from keras.layers.core import Dense
from keras.layers import Input
from keras import regularizers

def get_model():
    inputs = Input(shape=[5, ], name='input')
    z = Dense(48, activation='relu', kernel_regularizer=regularizers.l2(0.001))(inputs)
    z = Dense(48, activation='relu', kernel_regularizer=regularizers.l2(0.001))(z)
    outputs = Dense(1, activation='linear', kernel_regularizer=regularizers.l2(0.001), name='outputs')(z)
    return Model(inputs, outputs, name='model')

中間層の活性化関数は単調性を表現するためにreluを用いました。eluやseluを用いると予測値が小さくなる傾向が見られましたが、これは重みの初期値やepoch数などとの組み合わせで変わる気がしています。

損失関数の定義

次に損失関数の定義を行います。
今回は対数変換をするまえの目的変数に対してMSPEを適用したいと思います。そのため、対数変換後の目的変数に合わせた損失関数を定義します。と言っても難しいことはなく、MSPE内のyをexp(y)-1に変えるだけです。

import keras.backend as K
import tensorflow.compat.v2 as tf

def mspe_exp(y_true, y_pred):
    y_true_exp = K.exp(K.maximum(y_true, K.epsilon()))-1.
    y_pred_exp = K.exp(K.maximum(y_pred, K.epsilon()))-1.
    percentage_error = (y_true_exp - y_pred_exp) / y_true_exp
    squared_percentage_error = tf.square(percentage_error) 
    return tf.reduce_mean(squared_percentage_error)

web上の記事を見ると自作の損失関数を作成する際はkeras.backendのみを用いて実装している例が多いです。
ただ、keras.backendのみを用いると学習の途中でlossが突然nanになってしまい学習が不可能になることが度々あります。
ちょっと調べてみると同じような問題を抱えている人が多いようです。

Nan when training model wit RMSE/ RMSLE loss function

色々試した結果、私はkeras内の実装と同じように損失関数をtensorflow.compat.v2を用いて作成することで、この問題を回避できています。理由はわかりません。

モデルの学習・予測

モデルの学習と学習したモデルによる予想株価の計算をします。この際、out-of-foldを行うことによって、全データの一部を予測用データにして、その他を学習データにします。実装内容は完全にクロスバリデーションですが目的はモデルの評価するためではなく、ある銘柄の株価をほかの銘柄のデータで予測するためです。

from sklearn.model_selection import KFold
from keras.optimizers import Adam

result_list = [[],[]] # 各foldの結果を保存するリスト
kf = KFold(n_splits=5, shuffle=True)  # 分割条件を設定
for fold, (train_idx, test_idx) in enumerate(kf.split(y)):
    # 学習データと検証データに分ける
    X_train, X_test = X[train_idx], X[test_idx]
    y_train, y_test = y[train_idx][:, 1].astype(float), y[test_idx][:, 1].astype(float)

    # モデルの作成
    model = get_model()
    model.compile(    
        loss=mspe_exp, # 自作した損失関数
        optimizer=Adam(learning_rate=0.0008, amsgrad=True), 
    )

    # 学習
    temp_history = model.fit(
        X_train, y_train,
        validation_data=(X_test, y_test),
        batch_size=16,
        shuffle=True,
        epochs=100,
        verbose=0,
    )
    
    # 予測
    y_test_pred = model.predict(X_test)

    result_list[0].append(y[test_idx][:, 0])
    result_list[1].append(y_test_pred)
    
    print(f"fold: {fold}, loss: {temp_history.history['loss'][-1]}, val_loss: {temp_history.history['val_loss'][-1]} ")
    
    del model

なお、モデル構造を含むハイパーパラメータは最適化済みです。めんどくさかったです。

コードは省略しますが適当なCallbackを設定しておいて学習過程を確認しました。

最後にkfoldした結果をDataFrameにまとめます。
予測値を対数変換前の状態に戻すため、このタイミングでexp(y)-1しています。

idx = []
for idx_list in result_list[0]:
    idx.extend(idx_list)
idx = np.array(idx).ravel()

value = []
for value_list in result_list[1]:
    value.extend(value_list)
value = np.array(value).ravel()

pred = pd.DataFrame(np.array([idx, np.expm1(value)]).T, columns = ['コード', '理論株価'])
pred['理論株価'] = pred['理論株価'].astype(float)

結果の確認

予測した理論株価と銘柄の情報をマージします。

result = screener_result.merge(pred)
result['現在値'] = pd.to_numeric(result['現在値'].astype(str).str.replace(',', ''), errors='coerce')

結果をplotします。両軸ともに対数表示です。
コードは省略します。

下図は適当な点のデータを確認してみたもの。plotlyだと結果の確認が直感的でいいですよね。

以上によって各銘柄の現在の株価と理論株価の比較ができようになりました。分析結果は割安株や割高株の発掘などに活用します。

発展

今回用いたデータはファンダメンタルデータとして代表的なものではありますが極々小規模なものとなっています。より高い精度の分析を行おうとするのであれば、データやモデルを見直す必要があります。以下でいくつかアイデアを列挙します。大部分は私が実際に行っているものです。

データの種類を増やす
今回は売上などの代表的なものを説明変数に利用しましたが、キャッシュフローや優待利回りなどを追加することが思いつきます。また、すでに入手したデータを用いてROAや売上高経常利益率などを計算し追加するという選択肢もあります。
割安・割高の基準は業種ごとに異なるため、業種をダミー変数にして用いるのも手です。業種と言わず事業内容を分析して低次元ベクトル表現にしたものを入力することも思いつきます。
ちょっとテクニカル寄りになりますが信用残情報も利用が考えられます。

複数年度のデータを用いる
今回は直近の決算データのみを利用しましたが、過去〇年のデータ用いたり会社予想のデータを用いたりするなど時系列方向にデータを拡張することが考えられます。入手先としてはTDNETやEDINET、最近ではJ-Quants APIといったものが選択肢と挙げられます。特にJ-Quants APIは有料ではあるものの、整形済みで確度の高いデータが利用できます（配当データがもっと安価に手に入れば…）。
また、データを時系列方向に拡張した場合はモデルもそれに合わせたものが必要になります。可変長入力に対応するのであれば、ニューラルネットワーク系ではLSTMなどの利用が思いつきます。

初期値の影響を低減する
回帰モデルにニューラルネットワークを用いたため、予測値は重みの初期値の影響を受けます。初期値の違いが割安か割高かの分水嶺になったり、予想結果として外れ値のような値を出力してしまったりする可能性があります。こういった初期値の影響をなるべく除外するためには、学習から予測のプロセスを複数回実行することが思いつきます。例えば同じ学習データ・予測データで初期値を変えたモデルを30回次々学習すれば30個の予測値が得られます。最終的な予測値に30個の中央値を使ったり、平均と分散を計算し予測値の確からしさを定量的に計算したりすることで、モデルの初期値の影響を低減することができます。

モデルの改良
モデル構造やハイパーパラメータを変えることで精度が向上するかもしれません。

終わりに

いかがでしたでしょうか。
間違えている部分も有るかもしれませんが、皆様のヒントになる情報がひとつでも残せていたら幸いです。

お願い

会社の予想配当利回りが簡単に入手できる方法があればこっそり教えてください。予想配当は入手しやすいのですが、株式分割が予定されている銘柄だと配当と分割のタイミングを考慮してデータを整形しないといけないので。長らくトライしていますが、完璧な整形は私には難しくギブアップ気味です…。