みなさまこんにちは。
今回は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ファクターモデルを使った重回帰分析は、各特徴量の有効性を判断するために用いるそうです。
あとは、長期的な運用の目安として使ったりするのだとか。
将来の予測にも使えるのではと思い付きで試しにやってみましたが、撃沈しました。
最後まで見てくださって、ありがとうございました。