メーカー勤務・文系の私がProphetで株価予測してみた（Optuna）

# はじめに
はじめまして。
文系大学出身、メーカー勤務で営業をしている私が今回時系列データの分析に挑戦してみました。
初心者で時系列データ分析やProphet気になるなって方、是非読んでみてください:relaxed:

# 目次
* [きっかけ](#きっかけ)
* [概要](#概要)
* [実行環境](#実行環境)
* [Prophetとは](#prophetとは)
* [Optunaとは](#optunaとは)
* [全体の流れ](#全体の流れ)
* [インストールとインポート](#インストールとインポート)
* [銘柄の確認](#銘柄の確認)
* [データを取得・確認](#データを取得確認)
* [全体像をつかむ](#全体像をつかむ)
* [Prophet(デフォルト値)で予測](#prophetデフォルト値で予測)
* [Optunaでベストパラメーターを導きだす](#optunaでベストパラメーターを導きだす)
* [MAPE(平均絶対パーセント誤差)](#mape平均絶対パーセント誤差)
* [パラメーターごとの比較](#パラメーターごとの比較)
* [最後に](#最後に)
* [参考](#参考)


# きっかけ
そもそも営業の私がなぜPythonを勉強しているかというと…


ネットの広告…

そう、よくある広告につられてここまで来ちゃったんですね:sweat_drops:

データを見ることに割と興味はあったものの基礎知識はゼロ。
営業の未来予測はいつも肌感覚＆多少のPOS分析:see_no_evil:

月一くらいで株で稼げないかなとか思いながら営業生活続けている私ですが…
そもそも株って何買えばいいの？？？

というわけで今回私の気になる銘柄をPythonのProphetで予測することにしました:clap:

株価予測は営業の売上データと同じく時系列データなのでゆくゆくは売上データを分析して需要予測もしたいですね:raised_hands:

# 概要

* Pythonで時系列データを分析
* Prophetを活用
* Optunaでハイパーパラメーターを最適化

# 実行環境
Google colab
Python 3.10.2
Windows 11 Home

# Prophetとは
Facebook社が開発した時系列予測解析用のライブラリがProphetです。
トレンド・周期性・イベント・ノイズをそれぞれでモデル化し足し合わせたモデルです。

***y(t)=g(t)+s(t)+h(t)+ϵ(t)***
* g(t):トレンド
* y(t):周期性
* h(t):イベント
* ϵ(t):ノイズ


ポイントは

1. 柔軟なモデリングができる
2. 予測結果を解釈しやすい
3. ***簡単に予測モデルを作ることができる***:clap:

* [公式ドキュメント](https://facebook.github.io/prophet/)

# Optunaとは
 ***ハイパーパラメーターの最適化を自動化するためのもの***です。
ハイパーパラメーターのチューニングはすごく重要で、初心者には難しいですよね:sweat_drops:
そこを手助けしてくれるのがOptunaです:muscle:

* [公式ドキュメント](https://www.preferred.jp/ja/projects/optuna/)

# 全体の流れ
1. データの取得
2. 全体像をつかむため可視化
3. Prophet（デフォルト値）で予測
4. Optunaでベストパラメーターを導き出す
5. ベストパラメーターで予測し、項目3（デフォルト値）の予測と比較を行う

# インストールとインポート
Prophet & Plotlyインストール
（Plotlyはインタラクティブなグラフを作成するためのライブラリ）
```python
!pip install prophet
!pip install plotly==5.1.0
```
Optunaインストール
```python
!pip install optuna
import optuna
```
yfinanceインストール(Yahoo! Financeから情報を取得するためのAPI)
```python
!pip install yfinance
import yfinance as yf
yf.pdr_override()
```
インポート
```python
import pandas as pd
import numpy as np
from pandas_datareader import data as pdr
import os
from sklearn.metrics import mean_absolute_percentage_error
from prophet import Prophet
from prophet import plot
import matplotlib.pyplot as plt
```

# 銘柄の確認
株価の銘柄については以下より確認できます。
[日本取引所グループ](https://www.jpx.co.jp/markets/statistics-equities/misc/01.html)

今回は特定銘柄の株価予測をしていきたいと思います。
私はいつもお買い物マラソンでお世話になっている楽天グループ（コード：4755）の株についてみていきたいと思います:chart:

# データを取得・確認
Yahoo! financeからデータを取得します。
[yfinance](https://pypi.org/project/yfinance/)

pandas_datareaderから株価のデータを取得します。
銘柄コードの番号の後に『.T』をつけることで日本の銘柄の株価が取得できます。

```python
start = '2021-01-01'
end = '2022-07-31'

df = pdr.get_data_yahoo('4755.T', start=start, end=end)
df.head(10)
```
![df.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2779539/842f8b8a-70c4-d025-16bf-47b74661a307.png)

土日祝は株取引がないのでデータもないことがわかりますね。

次にreset_index()でインデックスを連番に振り直します。
```python
df_ri = df.reset_index()
print(df_ri)
```
![df_ri.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2779539/8246ba11-227d-083b-9386-e6ee4bca3983.png)

データを確認します。
```python
df_ri.info()
print(f'df_shape : {df_ri.shape}\n')
print(f'{df_ri.dtypes}\n')
```
![df_ri_info.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2779539/abc85687-c79a-4bb6-6bc9-be20e5b3a2c3.png)

Dateもdatetime64になっているのでこのまま進めていきます！

# 全体像をつかむ
まずは得られた情報がどのようなものか全体像を見てみることも大切なので、matplotlibで可視化します。
```python
plt.figure(figsize=(18, 10))

plt.plot(df_ri['Date'],df_ri.loc[:, "Open"], label='Open price')
plt.plot(df_ri['Date'],df_ri.loc[:, "High"], label='High price')
plt.plot(df_ri['Date'],df_ri.loc[:, "Low"], label='Low price')
plt.plot(df_ri['Date'],df_ri.loc[:, "Close"], label='Close price')
plt.xlabel('year')
plt.ylabel('stock price')
plt.legend(loc="upper right") #右上に凡例表示
plt.show()
```
![df_ri_plt.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2779539/9fc8f23b-c0d8-62d2-1c80-21bdfad0f921.png)

# Prophet（デフォルト値）で予測
終値に着目してProphetで分析していきたいと思います。

まずはProphetでは以下のカラムを用意します。
* ds ：日付　
* y  ：予測値

今回は以下の通りで進めていきます。
* ds : Date（取引日）
* y  : Close(終値）

というわけで、今回は使用しない不要なカラムは削除します！
```python
df2 = df_ri.drop(['Open', 'High', 'Low', 'Adj Close', 'Volume'], axis=1)
print(df2)
```
![2022-08-13.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2779539/d0962f8b-e9e5-68c7-5a47-8fb1b763feec.png)

trainデータとtestデータに分割します。
```python
df2.columns = ['ds', 'y']
train = df2[df2['ds'] < '2022-05-01']
train
```
![2022-08-13 (1).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2779539/eff70d56-b2b9-e979-88f2-b726fbeb6129.png)


```python
test = df2[df2['ds'] >= '2022-05-01']
test
```
![2022-08-13 (2).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2779539/464dbc69-f26f-9f15-7107-0812ca535807.png)

そしてパラメーターデフォルト値で学習します。
```python
m = Prophet() #インスタンス生成
m.fit(train)
forecast = m.predict(df2)
forecast.head(5)
```
![default.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2779539/56cb51b7-dccf-8106-d366-e2ce76ad5161.png)

* yhat:予測値 
* yhat_lower:予測値下限 
* yhat_upper:予測値上限

上記を可視化します。
```python
plot.plot_plotly(m, forecast, trend=True, changepoints=True)
```
trendとchangepointsはオプションです。
* 黒点がActual
* 赤線がTrend
* 青線がPredicted
* 青帯が予測値の幅（yhat_lower & yhat_upper)
![newplot (9).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2779539/4ba7f502-8cf8-635c-7445-1c115837ebe5.png)


予測値をさらに細かい要素に分け可視化します。
デフォルトでの出力内容は、トレンド、季節性（年）、季節性（週）です。
![newplot (10).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2779539/79479ff7-ec7d-c3b5-5a02-c587ad2bdb58.png)

ここで最後にOptunaでチューニングした際の予測精度との比較のためにMAPEを確認しておきます。
```python
test_forecast = forecast.tail(len(test))

test_mape = np.mean(np.abs((test_forecast.yhat - test.y)/test.y))*100
print(f'Default_param_mape : {test_mape}')
```
ここではDefault_param_mape : 17.60029625501451という数値が出ていますね。
# Optunaでベストパラメーターを導きだす
Prophetではパラメーターのチューニングが予測精度を高める上で重要になります。
パラメーターの最適化を行ってくれるOptunaを使っていきたいと思います:v:

パラメーターは

* growth: トレンドを表す関数（今回は線形モデルです）
* changepoints: トレンドの変化点のリスト
* n_changepoints: トレンドの変化点の数
* changepoint_range: トレンドの変化点を推測する幅
* yearly_seasonality: 年次の季節変動の有無(今回はFalse）
* weekly_seasonalityJ: 週次の季節変動の有無（今回はFalse）
* daily_seasonality: 日次の季節変動の有無（今回はFalse）
* holidays: 休日
* seasonality_mode: 周期性の傾向
* seasonality_prior_scale: 周期性の強さを表すパラメータ
* holidays_prior_scale: 休日の強さを表すパラメータ
* changepoint_prior_scale: トレンドの変化点の強さを表すパラメータ


パラメータの範囲を指定して、Optunaを実装します。
```python
def objective(trial):
    prophet_params = {'changepoint_range':trial.suggest_discrete_uniform('changepoint_range', 0.70, 0.95, 0.01),
                      'n_changepoints': trial.suggest_int('n_changepoints', 20, 40),
                      'changepoint_prior_scale': trial.suggest_discrete_uniform('changepoint_prior_scale', 0.001, 0.5, 0.001),
                      'seasonality_prior_scale': trial.suggest_discrete_uniform('seasonality_prior_scale', 0.01, 25, 0.5),
                      }
    m = Prophet(changepoint_range = prophet_params['changepoint_range'],
                n_changepoints = prophet_params['n_changepoints'],
                changepoint_prior_scale = prophet_params['changepoint_prior_scale'],
                seasonality_prior_scale = prophet_params['seasonality_prior_scale'],
                yearly_seasonality = False, 
                weekly_seasonality = False, 
                daily_seasonality = False, 
                holidays = None,
                growth = 'linear'
                 )
    m.fit(train)
    forecast = m.predict(df2)
    valid_forecast = forecast.tail(len(test))
    
    #MAPE=平均絶対パーセント誤差
    val_mape = np.mean(np.abs((valid_forecast.yhat - test.y) / test.y))*100
    print(f'Valid_mape : {val_mape}')


    return val_mape


study = optuna.create_study(direction="minimize") 
study.optimize(objective, n_trials=300)
 
print(f'Best param : {study.best_params}')
```
ここでdirectionについてですが、今回は予測精度の評価指標として[MAPE](#mape平均絶対パーセント誤差)を使用しています。
MAPEは数値（パーセンテージ）が小さければ小さいほど誤差が少ない（＝予測精度が高い）といえるので今回は`direction="minimize"`としています。

ベストパラメーターの出力です。
![best_param_2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2779539/a955fa69-268e-2154-b8ac-4f71d4f0defb.png)


上記で導き出されたベストパラメーターを用いて再度学習します。

```python
m = Prophet(**study.best_params)
m.fit(train)

forecast2 = m.predict(df2)
forecast2.head(5)
```
![3.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2779539/c225ffb6-44dc-ac5f-19e5-09b249b34250.png)

再度可視化します。
```python
plot.plot_plotly(m, forecast2, trend=True, changepoints=True)
```
![newplot (11).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2779539/8e21cdf4-4b7d-418f-3068-bfa9400c8747.png)


```python
plot.plot_components_plotly(m, forecast2)
```
![newplot (12).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2779539/87d966b8-1374-d2ef-530d-8d1e74b541d9.png)

ベストパラメーターのMAPEを出します。
```python
test_forecast = forecast2.tail(len(test))

test_mape = np.mean(np.abs((test_forecast.yhat - test.y)/test.y))*100
print(f'Best_param_mape : {test_mape}')
```
結果は14.912033396860421でした。


# MAPE(平均絶対パーセント誤差)
予測精度評価指標として用いていたMAPEに関しては以下の通りとなりました。
MAPEは「予測値と正解値との差を、正解値で割った値（＝パーセント誤差）」の絶対値を計算し、その総和をデータ数で割った値（＝平均値）を出力する関数です。

今回の結果は
* Default_param_mape : 17.60029625501451
* Best_param_mape: 14.912033396860421


解釈はDefault_paramに関しては『***平均して17％前後の誤差***』があり、Best_paramは『***平均して15%前後の誤差***』があるということになりますね。

MAPEは
* 10% : 高い予測精度
* 10%-20% : 良好な予測精度
* 20%-50% : 妥当な予測精度
* 50%以上 : 低い予測精度

最初に得たデフォルトパラメーターMAPEの数値『17.60029625501451』に対しベストパラメーターのMAPEは『14.912033396860421』だったので、***Optunaでハイパーパラメーターをチューニングすることで予測精度が多少改善している***といえるのではないでしょうか:clap:

# パラメーターごとの比較
ここでパラメーターごとにどのような変化が起き、予測精度の改善につながったのか見ていきたいと思います。
Prophetのパラメーターのデフォルト値は以下の通りです。
```python
params = {'growth': 'linear',
          'changepoints': None,
          'n_changepoints': 25,
          'changepoint_range': 0.8,
          'yearly_seasonality': 'auto',
          'weekly_seasonality': 'auto',
          'daily_seasonality': 'auto',
          'holidays': None,
          'seasonality_mode': 'additive',
          'seasonality_prior_scale': 10.0,
          'holidays_prior_scale': 10.0,
          'changepoint_prior_scale': 0.05,
          'mcmc_samples': 0,
          'interval_width': 0.80,
          'uncertainty_samples': 1000,
          'stan_backend': None}
```
以下が比較です。
| ハイパーパラメーター名 | デフォルト値 | ベスト値 |
|:-:|:-:|:-:|
| changepoint_range  |  0.8 |  0.7 |
|  n_changepoints |  25 |  24 |
|  changepoint_prior_scale |  0.05 |  0.23 |
|  seasonality_prior_scale |  10 | 19.51  |


 [Prophet公式](https://facebook.github.io/prophet/docs/diagnostics.html#hyperparameter-tuning)によると ***『changepoint_prior_scale』おそらくこれが最も影響のあるであろうパラメーター*** とされています.


* **changepoint_prior_scale**
小さすぎればアンダーフィットとなりモデル化されるべき分散がノイズ処理され、大きすぎればトレンドがオーバーフィットしてしまう可能性があります。

* **seasonality_prior_scale**
値が大きいと季節性が大きな変動に適合し、値が小さいと季節性の大きさが小さくなります。
パラメーターの妥当な範囲は [0.01, 10] とされていますが、今回は19.51ですのでその範囲を超えた調整が行われているようですね。

今回Optunaがこの辺りをしっかり最適化してくれたことが予測精度向上につながった可能性ありというところでしょうか:clap:


結論としては今回はいずれの結果も良好な予測精度の範囲内ではありましたが、***パラメーター最適化によって予測精度が平均2.7%前後の改善に繋がった***ということになりますね:kissing_heart:



# 最後に
今回はProphetを用い予測を行い、Optunaでハイパーパラメーターの最適化を図ることでさらに予測精度を高めることができました。
また、この2つを用いることで簡単に時系列データ分析を行うことが可能となりました。
…とは言え、株価予測は至難の業:sweat_drops:
予測精度を上げるためにはもっと勉強が必要ですね。
今後は時系列解析のRNNやLSTM、その他XGBoostやLightGBMなども学んで知識を深めていこうと思います。
最後までお読みいただきありがとうございました:hugging:

# 参考
* [Pythonで時系列解析・超入門（その4）Prophetモデルで予測する方法](https://www.salesanalytics.co.jp/datascience/datascience091/)
* [【Python機械学習】Prophetを用いて株価予測をやってみました(データサイエンス)](https://prog-you.com/prophet/)
* [「Prophet」とは――Facebook製時系列予測OSSは何が便利なのか](https://atmarkit.itmedia.co.jp/ait/articles/1904/19/news009_2.html)
* [将来予測で用いられるライブラリ「Prophet」](https://www.skillupai.com/blog/tech/prophet/)