Pythonで予測｜株価はいくらになる？【機械学習を使って予測する方法を解説：データ取得、データ前処理、モデル作成、モデル評価まで】

# ライブラリのインポート
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 = '2021-12-31'
data_master = data.DataReader('^N225', 'yahoo', start, end)
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(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()
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(2018,12,31)
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年〜2020年を学習用データとする
train = data_technical['2018-01-01' : '2020-12-31']
train

# 2021年をテストデータとする
test = data_technical['2021-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}')

valid_scores: [325.3625074145673, 169.61507596829318, 413.7658021675662, 201.0198662362686]
CV score: 277.4408129466738

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}')

RMSE: 435.1219386786387

先ほどのモデルを使用して、2021年の金曜日の株価を予測してみましょう。ここでは、先ほど作成した学習データとテストデータを使用します。
手順はこのように簡単です。
まず、線形回帰モデル(LinearRegression)のインスタンス化をします。
次に、学習データに対して線形回帰モデルで学習をします。第一引数には2018年〜2020年のデータの説明変数、第二引数には目的変数を指定します。
続けて、テストデータの説明変数を用いて、金曜日の株価を予測します。
最後に、予測した結果のy_predと実際の値であるテストデータの目的変数y_testから、予測精度のRMSEを算出します。

そして結果を表示しましょう。実行します。
RMSEはこのようになりました。交差検証の値よりも悪化したことがわかります。

可視化で予測と実際の値を確認

# 実際のデータと予測データをデータフレームにまとめる
df_result = test[['Close_next']]
df_result['Close_pred'] = y_pred
df_result

予測した値と実際の値がどの様になっているのか、目的変数だけをグラフにして確認してみましょう。
まず、実際のデータと予測データをデータフーレム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のグラフを描画してみます。
実際の値をオレンジ、予測した値を青の線で表示してみましょう。
株価の変動に対して、大まかな傾向は捉えていそうです

# 誤差を算出
df_result['diff'] = df_result['Close_pred'] - df_result['Close_next']

今度は、どの部分の乖離が大きいのかを確認するために、誤差だけをグラフにしてみましょう。
誤差は予測した値と実際の値の差分で計算します。赤い線のグラフにしましょう。また、誤差の目安を750円としてY軸に補助線を追加します。実行します。

# 誤差のグラフ作成
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()

では、誤差が大きかったところに注目して見てみましょう
どの部分も実際の値段が直前に大きく下に変動していることがわかります。
ここでは実際の値段よりも予測値が高い傾向ですが、他のこれらは、実際の変動よりも低い値段で予測していることがわかります。
こうして見ると、今回作成したモデルは、株価が下がる場合に誤差が大きくなる可能性があると考えられます。

つまり、この様な傾向を上手く捉えて予測モデルを改善することで、精度を向上させられると期待できます。
モデルの精度向上の方法としては、

• 説明変数の見直し
• 学習データ期間の見直し
• 予測モデルの見直し

などを行うことが有効です。説明変数については、ファイナンスデータのテクニカル分析を説明した動画の中でいくつか紹介しています。ぜひこちらも参考になさってみてください。

予測モデルの係数と切片を確認

# 予測モデルの係数を確認
coef = pd.DataFrame(model.coef_) # データフレームの作成
coef.index = X_train.columns     # 項目名をインデックスに設定
coef

さて、線形回帰モデルは、このような関数で予測値を表すことが可能です。
予測値 ＝A1×(説明変数1)+A2×(説明変数2)+・・・+An×(説明変数n)+K<
今回作成したモデルに対して、これらの値を確認してみます。
まず、係数は作成したモデルに対してcoef_メソッドで取得できます。ただし、これらがどの特徴量に紐づいているのかわかりにくいため、データフレームにしてみます。説明変数をインデックスにして表示してみましょう。
このように、説明変数に対する係数が表示されました。

# 予測モデルの切片を確認
model.intercept_

409.8866450575515

次に、切片です。intercept_メソッドで取得できます。実行します。

この関数を使って予測をした、ということです。
ここで、1点注意いただきたいのは、各説明変数の係数の大きさが相関の強さを示さない点です。
理由は、学習データの各説明変数ごとの分布が揃っていない為です。

# X_train基本統計量の確認
X_train.describe()

では、説明変数の分布を確認するために、モデルの学習に使用したX_trainの基本統計量を確認しましょう。
標準偏差stdを見ると、小さいもので実体'Body'の162に対して、大きいもので最安値'Low'の1602というように、説明変数間でばらつきに差があることがわかります。
原因は、データのスケールが異なることにあります。より影響の度合いを考慮したい場合は、データを標準化するとよいでしょう。

エンディング

いかがでしたでしょうか？
今回は重回帰分析での予測方法について、ファイナンスデータを使って説明しました。
重回帰分析については、別の動画で詳しく説明します。もう少々お待ちくださいませ。
また、キノコードではわかりやすく見飽きない動画作成を心がけています。
チャンネル登録がまだの方は、新着通知も届きますので、ぜひチャンネル登録をお願いします。それでは次のレッスンでお会いしましょう。