重回帰分析を使って、株価の予測を実装してみる。
やりたいこと
ノジマの株価予測
欲しい結果
- 実装手順
- ライブラリのインポート
- 株価データの取得
- 曜日情報を追加
- グラフの描画
- データのコピー
- 移動平均データの追加
- 特徴量の描画
- 目的変数の作成
- 必要なカラムを抽出
- データの分割
- モデル作成・精度検証のためのライブラリのインポート
- 学習、予測データの作成
- 精度確認
- 実測値と予測値の比較
実装コード
# ライブラリのインポート
from datetime import datetime, timedelta
import numpy as np
import pandas as pd
from pandas_datareader import data
import matplotlib.pyplot as plt
%matplotlib inline
# ワーニングを非表示にする設定(任意)
import warnings
warnings.simplefilter('ignore')
# 最大表示行数の指定(任意:ここでは10行を指定)
pd.set_option('display.max_rows', 10)
# pandas_datareaderを使って、2018年始から2021年末までの日経平均株価データの取得
start = '2018-01-01'
end = '2022-06-03'
data_master = data.DataReader('7419.JP', 'stooq', start, end)
data_master
data_master = data_master[::-1]
data_master
# 曜日情報を追加(0:月曜日〜4:金曜日)
data_master['weekday'] = data_master.index.weekday
data_master
plt.figure(figsize=(10,6))
plt.plot(data_master['Close'], label='Close', color='orange')
plt.xlabel('Date')
plt.ylabel('JPY')
plt.legend()
plt.show()
# data_techinicalにデータをコピー
data_technical = data_master.copy()
# 移動平均を追加
SMA1 = 5 #短期5日
SMA2 = 10 #中期10日
SMA3 = 15 #長期15日
data_technical['SMA1'] = data_technical['Close'].rolling(SMA1).mean() #短期移動平均の算出
data_technical['SMA2'] = data_technical['Close'].rolling(SMA2).mean() #中期移動平均の算出
data_technical['SMA3'] = data_technical['Close'].rolling(SMA3).mean() #長期移動平均の算出
# 特徴量を描画して確認
plt.figure(figsize=(10, 6))
plt.plot(d2ata_technical['Close'], label='Close', color='orange')
plt.plot(data_technical['SMA1'], label='SMA1', color='red')
plt.plot(data_technical['SMA2'], label='SMA2', color='blue')
plt.plot(data_technical['SMA3'], label='SMA3', color='green')
plt.xlabel('Date')
plt.ylabel('JPY')
plt.legend()
plt.show()
# 特徴量を描画して確認(x軸の拡大)
plt.figure(figsize=(10, 6))
plt.plot(data_technical['Close'], label='Close', color='orange')
plt.plot(data_technical['SMA1'], label='SMA1', color='red')
plt.plot(data_technical['SMA2'], label='SMA2', color='blue')
plt.plot(data_technical['SMA3'], label='SMA3', color='green')
plt.xlabel('Date')
plt.ylabel('JPY')
plt.legend()
xmin = datetime(2018,1,1)
xmax = datetime(2022,6,3)
plt.xlim([xmin,xmax])
plt.show()
# OpenとCloseの差分を実体Bodyとして計算
data_technical['Body'] = data_technical['Open'] - data_technical['Close']
# 前日終値との差分Close_diffを計算
data_technical['Close_diff'] = data_technical['Close'].diff(1)
# 目的変数となる翌日の終値Close_nextの追加
data_technical['Close_next'] = data_technical['Close'].shift(-1)
data_technical
data_technical = data_technical.dropna(how='any')
data_technical
# 木曜日のデータを抜き出す
data_technical = data_technical[data_technical['weekday'] == 3]
data_technical
# 必要なカラムを抽出
data_technical = data_technical[['High', 'Low', 'Open', 'Close', 'Body',
'Close_diff', 'SMA1', 'SMA2', 'SMA3', 'Close_next']]
data_technical
# 2018年〜2021年を学習用データとする
train = data_technical['2018-01-01' : '2021-12-31']
train
# 2022年をテストデータとする
test = data_technical['2022-01-01' :]
test
# 学習用データとテストデータそれぞれを説明変数と目的変数に分離する
X_train = train.drop(columns=['Close_next']) #学習用データ説明変数
y_train = train['Close_next'] #学習用データ目的変数
X_test = test.drop(columns=['Close_next']) #テストデータ説明変数
y_test = test['Close_next'] #テストデータ目的変数
# 線形回帰モデルのLinearRegressionをインポート
from sklearn.linear_model import LinearRegression
# 時系列分割のためTimeSeriesSplitのインポート
from sklearn.model_selection import TimeSeriesSplit
# 予測精度検証のためMSEをインポート
from sklearn.metrics import mean_squared_error as mse
# 時系列分割交差検証
valid_scores = []
tscv = TimeSeriesSplit(n_splits=4)
for fold, (train_indices, valid_indices) in enumerate(tscv.split(X_train)):
X_train_cv, X_valid_cv = X_train.iloc[train_indices], X_train.iloc[valid_indices]
y_train_cv, y_valid_cv = y_train.iloc[train_indices], y_train.iloc[valid_indices]
# 線形回帰モデルのインスタンス化
model = LinearRegression()
# モデル学習
model.fit(X_train_cv, y_train_cv)
# 予測
y_valid_pred = model.predict(X_valid_cv)
# 予測精度(RMSE)の算出
score = np.sqrt(mse(y_valid_cv, y_valid_pred))
# 予測精度スコアをリストに格納
valid_scores.append(score)
print(f'valid_scores: {valid_scores}')
cv_score = np.mean(valid_scores)
print(f'CV score: {cv_score}')
model = LinearRegression()
model.fit(X_train, y_train)
y_pred = model.predict(X_test)
score = np.sqrt(mse(y_test, y_pred))
print(f'RMSE: {score}')
# 実際のデータと予測データをデータフレームにまとめる
df_result = test[['Close_next']]
df_result['Close_pred'] = y_pred
df_result
# 実際のデータと予測データの比較グラフ作成
plt.figure(figsize=(10, 6))
plt.plot(df_result[['Close_next', 'Close_pred']])
plt.plot(df_result['Close_next'], label='Close_next', color='orange')
plt.plot(df_result['Close_pred'], label='Close_pred', color='blue')
plt.xlabel('Date')
plt.ylabel('JPY')
xmin = df_result.index.min()
xmax = df_result.index.max()
plt.legend()
plt.show()
# 誤差の算出
df_result['diff'] = df_result['Close_pred'] - df_result['Close_next']
# 誤差のグラフ化
plt.figure(figsize=(10, 6))
plt.plot(df_result[['diff']])
plt.plot(df_result['diff'], label='diff', color='red')
plt.xlabel('Date')
plt.ylabel('error')
plt.hlines(0, xmin, xmax, color='gray', linestyle='--')
plt.hlines(750, xmin, xmax, color='gray', linestyle=':')
plt.hlines(-750, xmin, xmax, color='gray', linestyle=':')
plt.legend()
plt.show()
ぼく的ポイント
data_technical['SMA1'] = data_technical['Close'].rolling(SMA1).mean() #短期移動平均の算出
内容:短期移動平均の列をコピーしたデータに追加
valid_scores = []
tscv = TimeSeriesSplit(n_splits=4)
内容:効果検証結果を格納するための空のリストを作成
↓
TimeSeriesSplitをインスタンス化
得られた結果
valid_scores: [50.3596830884276, 51.869664757064214, 61.24507532582339, 87.13981033440787]
CV score: 62.653558376430766
結果からわかること
- グラフを見る限り、実測値と予測値に大きな差異はない
- 2000~3000円の株価に対して、60円ぐらいの誤差がある
まとめ
手順としてはシンプルだなと思った。
しかし、データの加工、整形などを自分1人ではまだとてもできないので、数をこなして精度をあげていきたい。機械学習のように、自分自身も学習を繰り返して精度の高いアウトプットを出せるようにしていきたい。