みなさまこんにちは。
今回は3ファクターモデルを使って重回帰してみます。
長くなるので、理論の詳細は割愛します。

import investpy
from datetime import datetime, timedelta
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
%matplotlib inline
import scipy.stats as stats
import statsmodels.formula.api as smf

今回使う株の指標は、investig.comから持ってくるので、investpyを使います。
pipでインストールできるので、気軽に使えます。

topix_df = jp_index.query("name.str.contains('TOPIX')", engine="python")
print(len(topix_df))
topix_df[["name","symbol"]][:10]

name	symbol
TSE TOPIX17 Wholesale Stock Price	TRAD17.T
TOPIX Banks	IBNKS.T
TSE TOPIX17 Transportation Stock Price	TRAN17.T
TOPIX Air Transportation	IAIRL.T
TSE TOPIX17 Steel Stock Price	STEL17.T
TSE TOPIX17 IT Services Stock Price	ITSV17.T
TSE TOPIX17 Real Estate Stock Price	RLTY17.T
TSE TOPIX17 Power Gas Stock Price	PWGA17.T
TSE TOPIX17 Pharmacetical Stock Price	PHAM17.T
TSE TOPIX17 Machinery Stock Price	MACH17.T
試しに日本の指標を確認しますが、今回使うTopixのものだけを抽出して見てみます。
これは一部抜粋ですが、こんな感じでデータを持ってこれます。

def load_index(lb, idx, country):
#取得するデータの範囲を日付で指定する
    delta = timedelta(days=lb)
    today = datetime.today().date()
    start = datetime.strftime((today - delta), '%d/%m/%Y')

    end = datetime.today().date
    end = datetime.today().strftime('%d/%m/%Y')
#指標データをダウンロードする   
    df = investpy.get_index_historical_data(index=idx, country=country, 
                                               from_date=start, to_date=end)
    return df
lb = (365*10)
df_index = load_index(lb=lb, idx="TOPIX Land Transportation", country="japan")

	Open	High	Low	Close	Volume	Currency
2011-02-28	1041.5	1046.1	1035.9	1044.1	0	JPY
2011-03-01	1048.0	1057.8	1048.0	1054.2	0	JPY
2011-03-02	1051.2	1052.8	1039.2	1039.3	0	JPY
2011-03-03	1044.6	1047.6	1043.8	1044.4	0	JPY
2011-03-04	1051.0	1054.3	1045.1	1045.1	0	JPY
10年分のTOPIXのデータを持ってきますが、分析したい株に関連する指標であるTOPIX Land Transportationを指定します。
日本の指標なので、countryには欲しい指標の国を指定。
指標の名前はシンボルでなく、nameで指定します。

データの取得コマンドの連発は、サーバーの負荷となるのでやめましょう。

#ダウンロードした3ファクターのデータを読み込み
def load_3fact():
    df = pd.read_csv("Japan_3_Factors_Daily.csv")
    df = df.rename(columns={"Unnamed: 0":"Date"})
    df["Date"] = pd.to_datetime(df["Date"], format='%Y%m%d')
    df.index = df["Date"]
    del df["Date"]
    return df
fama = load_3fact()
fama

	Mkt-RF	SMB	HML	RF
1990-07-02	0.85	0.38	-0.06	0.03
1990-07-03	0.07	0.72	0.30	0.03
1990-07-04	1.45	0.52	0.26	0.03
1990-07-05	-0.64	0.48	0.38	0.03
1990-07-06	-0.02	0.40	0.31	0.03
こちらは、3ファクターモデルのデータとなっています。
データは、https://mba.tuck.dartmouth.edu/pages/faculty/ken.french/data_library.html
からダウンロードしました。

#investing.comから個別株のデータをダウンロード
def load_stock(lb, code, country):
    delta = timedelta(days=lb)
    today = datetime.today().date()
    start = datetime.strftime((today - delta), '%d/%m/%Y')

    end = datetime.today().date
    end = datetime.today().strftime('%d/%m/%Y')
    
    df = investpy.get_stock_historical_data(stock=code, country=country, 
                                               from_date=start, to_date=end)
    return df
df_stock = load_stock(lb=lb, code="9020", country="japan")
df_stock = df_stock[["Close"]]
df_stock

	Close
2011-02-28	5690
2011-03-01	5820
2011-03-02	5710
2011-03-03	5790
2011-03-04	5810
Topixを取得した時のように、investing.comからJR東日本のデータを持ってきます。
個別株を指定する際は、銘柄コードを文字列で指定します。
今回は終値だけ使います。

#使用するデータをデータフレーム化
train_data = df_stock.join(fama)
train_data["Topix"] = df_index["Close"]
train_data = train_data[["SMB", "HML", "Topix", "Close"]]
train_data = train_data.rename(columns={"Close":"Stock_Price"})
train_data

	SMB	HML	Topix	Close
2011-02-28	0.64	0.02	1044.1	5690
2011-03-01	-0.39	0.31	1054.2	5820
2011-03-02	0.53	-0.31	1039.3	5710
2011-03-03	0.18	0.20	1044.4	5790
2011-03-04	-0.19	-0.13	1045.1	5810
使うデータを一つのデータフレームにします。
マーケットポートフォリオにはTOPIXの変化率を使用して、SMB、HMLを説明変数とします。
目的変数は個別株のリターンを設定します。

#欠損値の確認
train_data.isnull().sum()
[out]
SMB            36
HML            36
Topix           0
Stock_Price     0

HML、SMBのデータが12月末までしかないので、欠損値が発生しています。

#変化率を算出する
train_data[["Topix", "Stock_Price"]] = train_data[["Topix", "Stock_Price"]].pct_change()
train_data = train_data.dropna()
train_data

pct_changeで変化率、リターンを算出しています。
欠損値がまた発生するので、ここで欠損値を削除します。

stats.probplot(train_data["SMB"], dist="norm", plot=plt)
plt.title("SMB")
plt.show()

QQプロットで可視化します。
HMLとTopixは


正規分布ではなさそうですが、、念のため正規性の検定を行います。

def ks_test(data):
    result = stats.kstest(data, "norm")
    return result
result_smb = ks_test(train_data["SMB"])
print("smb p-value:" ,result_smb[1])

result_hml = ks_test(train_data["HML"])
print("hml p-value:" ,result_hml[1])

result_topix = ks_test(train_data["Topix"])
print("topix p-value:" ,result_hml[1])
[out]
smb p-value: 2.3568451645966648e-54
hml p-value: 7.770824371895942e-58
topix p-value: 7.770824371895942e-58

p値が5%を余裕で下回っているため、帰無仮説が棄却され、正規分布であることが否定されます。
検定や推定などの統計的手法は正規分布を前提としているため、ちょっと困ります。
そこで、中心極限定理を適用することにします。
今回使用したデータは母集団から抽出した標本とし、サイズを大きくしていけば正規分布に近似することにします。
ちょっと強引ですが・・・

#重回帰分析を行う
formula = "Stock_Price ~ SMB + HML + Topix"
result = smf.ols(formula, train_data).fit()
result.summary()
[out]
 OLS Regression Results                            
==============================================================================
Dep. Variable:            Stock_Price   R-squared:                       0.792
Model:                            OLS   Adj. R-squared:                  0.792
Method:                 Least Squares   F-statistic:                     3052.
Date:                Thu, 25 Feb 2021   Prob (F-statistic):               0.00
Time:                        20:38:00   Log-Likelihood:                 8447.1
No. Observations:                2408   AIC:                        -1.689e+04
Df Residuals:                    2404   BIC:                        -1.686e+04
Df Model:                           3                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
Intercept     -0.0001      0.000     -1.010      0.313      -0.000       0.000
SMB            0.0008      0.000      2.812      0.005       0.000       0.001
HML            0.0011      0.000      4.062      0.000       0.001       0.002
Topix          1.0809      0.012     93.540      0.000       1.058       1.104
==============================================================================
Omnibus:                      817.226   Durbin-Watson:                   2.182
Prob(Omnibus):                  0.000   Jarque-Bera (JB):            49433.551
Skew:                          -0.768   Prob(JB):                         0.00
Kurtosis:                      25.143   Cond. No.                         78.2
==============================================================================

重回帰を行った結果、(自由度調整済み)決定係数は約79%です。
100%に近いほど、あてはまりがよい結果となります。
各説明変数に様々な数字が並びますが、P>|t|に注目します。
これはp値なので、この数値が有意水準を下回れば、その変数が目的変数に影響を持つことになります。
そうすると、切片を意味するIntercept以外は日時リターンの予測に影響を持っていることになります。

考察
決定係数を見ると約8割の精度ですが、今回の手法だと将来の株価予測が難しいです。
時系列データを単純に回帰しているだけで、説明変数と目的変数に時間差がないんですよね。
そのため、得られた回帰式を将来の予測に使うのは難しそうです。
指標と個別株のリターンにある程度の時差を与えると回帰の精度が下がりました。
例えばどちらかを始値にして、もう片方を終値にすると時差が生まれますが、精度が落ちます。
将来の予測をする場合は、時間差があっても予測に有効な特徴量を見つけるしかないですね。

まとめ
この3ファクターモデルを使った重回帰分析は、各特徴量の有効性を判断するために用いるそうです。
あとは、長期的な運用の目安として使ったりするのだとか。
将来の予測にも使えるのではと思い付きで試しにやってみましたが、撃沈しました。

最後まで見てくださって、ありがとうございました。