単回帰分析を利用して、時価総額を求めてみる #Python3

あまり、実用的ではありませんが単回帰分析を利用して予測を行ってみたいと思います。

単回帰分析の目的とは

単回帰分析では、以下の数式（モデル）で予測できると仮定して、
予想値と実測値との差分が一番小さくなるような変数a、変数bの値を求めます。

$$ 予測値 = a * 入力値x + b $$

比例の関係であることが前提ですので、相関関係が高い値同士でのみ単回帰分析が可能です。
例えば、営業利益と時価総額、部屋の広さと家賃等です。

中心化

データの中心化というテクニックがあります。中心化とは、その名の通りデータを中心に集めることを言います。具体的には、データを全て、平均値で引いた状態にします。
こうすることによって、どんなメリットがあるのでしょうか。

もともと、「$ 予測値 = a * 入力値x + b $」で予測していたモデルを「$ 予測値 = a * 入力値x $」で考えることができるようになります。

式変形

単回帰分析では、予測値と実測値の差分が一番小さくなるような変数aの値を求めます。
それは、以下の式の値が最小化するような a の値を求めることを言います。

\sum_{n=1}^{N} (実測値_n - 予測値_n)^2

この式では、実測値と予測値の差分の２乗の合計を求めています。
なぜ２乗するのかというと、２つの理由があります。

実測値と予測値の差分がマイナスの可能性もあること。
２乗することによって式が二次関数で表せるようになり最小となる値を微分することによって簡単に求めることができるようになる。

では、上記の式を変形していきます。
$ 予測値_n $は、$ a * 入力値x_n $ とも表示できます。
そのため以下のようにも表せます。

\sum_{n=1}^{N} (実測値_n - a * 入力値x_n)^2

また、展開するとこうなります。

\sum_{n=1}^N (実測値_n^2 - 2 * 実測値_n* a * 入力値x_n   + a^2 * 入力値x_n^2)

また、シグマの足し算は展開できるので

\sum_{n=1}^N (実測値_n^2) - 2 * \sum_{n=1}^N (実測値_n * 入力値x_n)* a  + a^2  * \sum_{n=1}^N (入力値x_n^2)

この式を微分したものが0になるようなaを求めると最適解となります。

\frac{\partial}{\partial a}(\sum_{n=1}^N (実測値_n^2) - 2 * \sum_{n=1}^N (実測値_n * 入力値x_n)* a  + a^2  * \sum_{n=1}^N (入力値x_n^2)) = 0

aで微分するので、aが含まれない個所は、0になります。
aをaで微分すると1、$ a^2 $をaで微分すると2aなので

- 2 * \sum_{n=1}^N (実測値_n * 入力値x_n)  + 2 * a  * \sum_{n=1}^N (入力値x_n^2) = 0

a  = \frac{\sum_{n=1}^N (実測値_n * 入力値x_n)}{\sum_{n=1}^N (入力値x_n^2)}

となります。
これは、あくまでもデータが中心化されていることが前提なので注意してください。

実装をしてみる

では、今回は、株価（時価総額）を予測するための変数aを求めたいと思います。

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

データの読み込み

今回使用するデータは、時価総額があるCSVと経常利益があるCSVが別々なのでそれぞれ別に読み込みます。

df1 = pd.read_csv(filepath_or_buffer="~\japan-all-stock-data.csv", encoding="ms932", engine="python", sep=",")
df2 = pd.read_csv(filepath_or_buffer="~\japan-all-stock-financial-results.csv", encoding="ms932", engine="python", sep=",")

df1.head(3)

	SC	名称	市場	業種	時価総額(百万円)	発行済株式数	配当利回り	1株配当	PER（予想）	PBR（実績）	EPS（予想）	BPS（実績）	最低購入額	単元株	高値日付	年初来高値	安値日付	年初来安値
0	1	日経225（日経平均株価）	東証	株価指数	-	-	-	-	-	-	-	-	-	-	-	-	-	-
1	2	TOPIX（東証株価指数）	東証	株価指数	-	-	-	-	-	-	-	-	-	-	-	-	-	-
2	1301	極洋	東証一部	水産・農林	42948	10928283	1.53	60.00	15.29	1.49	257.01	2644.98	393000	100	2018/1/10	4460	2017/1/18	2681

df2.head(3)

	SC	名称	決算期	決算発表日(本決算)	売上高(百万円)	営業利益(百万円)	経常利益(百万円)	当期利益(百万円)	総資産(百万円)	自己資本(百万円)	資本金(百万円)	有利子負債(百万円)	自己資本比率	ROE	ROA	発行済株式数
0	2788	アップルインターナショナル	2017/12	2018/02/19	13634	325	335	204	8230	7357	4322	1476	89.40	2.86	2.58	13841400
1	5101	横浜ゴム	2017/12	2018/02/19	668049	51933	52887	35217	929029	383925	38909	315916	41.30	9.61	3.84	169549081
2	6278	ユニオンツール	2017/12	2018/02/19	23188	3698	3718	2655	57605	52440	2998	-	91.00	5.22	4.81	20788590

データの結合

まず、別々に読み込んだデータフレームを結合します。
結合するためには、データフレームのmerge関数を利用します。
merge関数では、SQLの外部結合みたいな指定もできます。今回は、SQLでいう内部結合で結合キーはSCになります。

df_merged = pd.merge(df1, df2, how='inner', on='SC')
df_merged.head(3)

	SC	名称_x	市場	業種	時価総額(百万円)	発行済株式数_x	配当利回り	1株配当	PER（予想）	PBR（実績）	...	経常利益(百万円)	当期利益(百万円)	総資産(百万円)	自己資本(百万円)	資本金(百万円)	有利子負債(百万円)	自己資本比率	ROE	ROA	発行済株式数_y
0	1301	極洋	東証一部	水産・農林	42948	10928283	1.53	60.00	15.29	1.49	...	3709	2422	97391	24976	5664	50919	25.60	10.19	2.52	10928283
1	1332	日本水産	東証一部	水産・農林	171524	312430277	1.46	8.00	8.55	1.22	...	24884	14216	451876	120973	30685	207749	26.80	13.17	3.17	312430277
2	1333	マルハニチロ	東証一部	水産・農林	169292	52656910	1.24	40.00	10.92	1.46	...	27874	15446	501303	100664	20000	272208	20.10	16.62	3.13	52656910

欠損値の削除

次は、欠損値を削除していきます。また、時価総額に比例する変数xを何にするか決めていないので、まずは、「経常利益(百万円)」、「営業利益(百万円)」、「当期利益(百万円)」の3つを候補にしてこの3つのカラムで欠損値がないデータを対象にしようと思います。
まずは、欠損値がある行を含む行数を確認します。

len(df_merged)

# 欠損値を削除していきます。
# 対象のデータでは、欠損値が「-」で登録されているので「-」以外の行を対象にします。
df_removed = df_merged[df_merged['経常利益(百万円)'] != '-']
df_removed = df_removed[df_removed['営業利益(百万円)'] != '-']
df_removed = df_removed[df_removed['当期利益(百万円)'] != '-']
len(df_removed)

これで約100行を削除できました。

型の変換

カラムの型を確認します。

df_removed.dtypes [['経常利益(百万円)', '営業利益(百万円)', '当期利益(百万円)', '時価総額(百万円)']]

経常利益(百万円)     int32
営業利益(百万円)     int32
当期利益(百万円)     int32
時価総額(百万円)    object
dtype: object

型を確認してみると、objectとなっています。
これは、欠損値として「-」が登録されていたので数値と文字が混じっていると判断され何でもという意味のobjectになっています。
これを、数値型に変換します。

# 各カラムを数値に変換します。
df_removed = df_removed.astype({'SC':'object', '経常利益(百万円)':'int', '営業利益(百万円)':'int', '当期利益(百万円)':'int', '時価総額(百万円)':'int'})
df_removed.dtypes [['経常利益(百万円)', '営業利益(百万円)', '当期利益(百万円)', '時価総額(百万円)']]

経常利益(百万円)    int32
営業利益(百万円)    int32
当期利益(百万円)    int32
時価総額(百万円)    int32
dtype: object

相関関係を見てみる

相関関係をみて、時価総額を求めるために最も関連のあるカラムを選択します。
相関関係を判定するために相関係数を使用すると良いです。1に近いほど関係性があり、0に近いほど関連がありません。
DataFrameでは、corrというメソッドを使用するとみることができます。

df_target = df_removed [['SC', '名称_x', '市場', '業種', '経常利益(百万円)', '営業利益(百万円)', '当期利益(百万円)', '時価総額(百万円)']]
df_target.corr()

	経常利益(百万円)	営業利益(百万円)	当期利益(百万円)	時価総額(百万円)
経常利益(百万円)	1.000000	0.989025	0.893595	0.908569
営業利益(百万円)	0.989025	1.000000	0.895535	0.906487
当期利益(百万円)	0.893595	0.895535	1.000000	0.842074
時価総額(百万円)	0.908569	0.906487	0.842074	1.000000

この値を見てみると時価総額(百万円)に対して、経常利益(百万円)が一番相関関係が高いので経常利益(百万円)から時価総額を予測するための変数aを求めたいと思います。

# センタリングを実施
# mean関数を使用することで平均値を求めることができます。
df_target.mean()

SC             5732.019755
経常利益(百万円)     14613.265442
営業利益(百万円)     14226.430161
当期利益(百万円)      9506.916249
時価総額(百万円)    183635.331107
dtype: float64

# dataFrame同士を引き算するだけで、各行づつ自動で計算してくれます。
x = df_target['経常利益(百万円)']
y = df_target['時価総額(百万円)']
xc = x - x.mean()
yc = y - y.mean()

a  = \frac{\sum_{n=1}^N (実測値_n * 入力値x_n)}{\sum_{n=1}^N (入力値x_n^2)}

なので、aは以下で求められます。

a = (yc * xc).sum()/(xc * xc).sum()
a

9.3824176098557999

# matplotlibを使用して、計算した結果を表示します。
# 実際のデータの位置を散布図で表示します。
plt.scatter(x, y, label='actual')
# 変数aに基づき、予測値の一次関数のグラフを表示します。
plt.plot(x, a*x, label='predict', color='red')
# 判例をグラフにプロット
plt.legend()
plt.show()

このように赤線が予測のラインです。

以下の計算で実際の時価総額を予測できます。
業種を絞って、計算してみるとか、もっと工夫しないと全然精度の高い予測はできませんが、ひとまずこんな感じで予測できました。

predict = a * (df_target['経常利益(百万円)'] - df_target['経常利益(百万円)'].mean()) + df_target['時価総額(百万円)'].mean()

次は、重回帰分析を勉強していきたいと思います。