【R言語と時系列解析】Prophetパッケージを利用した、時系列データの未来予測: ブログ利用者数の推移予測とか、ビットコインや日経225の価格トレンド推移予測とか (original) (raw)

はじめに
ブログ利用ユーザー数の推移を未来予測する
日経225インデックスの将来推移を予測してみたところ
- 2012年以降のデータで、prophetで株価トレンドを予測する
ビットコイン(BTC)の価格推移を予測してみたところ
まとめ
補足
参考資料

はじめに

Prophetパッケージは、Meta社（当時、Facebook社）が2017年頃に開発した時系列予測のオープンソースのライブラリです。 R実装版もあり、今回そのパッケージを使用しました。

実際、私も実務で経験しましたが、、時系列データの予測アルゴリズムの構築は途轍もなく困難を伴います。その行為はまさしく、途方もなく、雲を掴む感じですね。

Prophetパッケージは、そのような時系列モデル構築の難易度を下げて、ユーザーをアシストしてくれます。また、多くの場合、モデルの細かいチューニングをしなくても、デフォルト設定でも案外良好なアウトプットが得られるのも魅力的なところです。

この記事では基本的な関数紹介が中心ではなくて、より実践的な解析内容を中心に紹介します。実データとしては、このブログ「京橋のバイオインフォマティシャンの日常」の利用ユーザー数のデータをもとにして、将来のユーザー数の推移予測をやってみました。2022年11月現在、ブログのユーザー数は月で7000-8000人あたりにいます。

また、それに加えて、日経225やビットコインの推移予測の実施事例も紹介したいと思います。

実行環境

macOS Big Sur (バージョン11.5.2) MacBook Air (M1, 2020) チップ Apple M1 メモリ 16GB

Prophetのデータフォーマット

まずは、Prophetのデータ形式についてです。

Prophetのインプット・データには、ds列とy列という、2列のデータを持つ、データフレームを準備します。

ds列 (=datestamp列) は時系列の列で、日付ならYYYY-MM-DD、タイムスタンプなら YYYY-MM-DD HH:MM:SS というスタイルになります。y列は時系列のそれぞれに対応する数値データで、予測したい測定値の過去データとなります。

データフォーマットも、非常にシンプルですね。

ブログ利用ユーザー数の推移を未来予測する

実行には、prophetパッケージのprophet関数を使います。

prophet関数の最初の引数dfには、過去データのデータフレームを指定します。このデータフレームは、過去データの日付を含むds列、及び過去データの数値データを含むy列を持つ必要があります。

それでは、実際のこのブログの利用ユーザー数のデータを使って、予測を実行してみます。この時系列データは「2020/4/1」か「ら2022/11/12」まであります。

また、過去1日間の訪問ユーザー数(y_1day列)、過去7日間の訪問ユーザー数(y_7day列)、過去14日間の訪問ユーザー数(y_14day列)、過去28日間の訪問ユーザー数(y_28day列)の列があります。

過去14日間の訪問ユーザー数(y_14day列)を使って、ブログ利用ユーザー数の推移を予測する

まずは、過去14日間の訪問ユーザー数(y_14day列)のデータに対して、prophetを実行しました。

ここで、増加し続ける周期性を示すデータの予測を行う場合には、「seasonality.mode = "multiplicative"」を指定します。

また、周期性の誤差を考慮するには、mcmc.samplesというパラメータを設定することで、ベイズ推論を行えます。結構計算時間がかかるので、今回は、デフォルトのMAP推定にしています。

それでは、実行していきます。

url <- "https://gist.githubusercontent.com/kumeS/c677d0b3ceb56f74cd109aae6acfd1a8/raw/6cf600582b8ba5bc3812d2a600de959c01548eca/Blog_Users_20200401-20221112.csv"

df <- read.csv(url, header = T)

head(df)

df_14day <- data.frame(ds=as.Date(df$ds), y=as.numeric(sub(",", "", df$y_14day)), row.names=1:nrow(df))

head(df_14day)

m <- prophet(df=df_14day, growth = "linear", n.changepoints = 100, changepoint.range = 1, changepoint.prior.scale=0.05, seasonality.mode = "multiplicative", interval.width = 0.95)

class(m)

prophet関数を実行して、mというモデル変数を出力します。

次に、make_future_dataframe関数を使って、予測する期間のデータフレームを生成します。 freqの引数としては、day、week、month、quarter、yearなどを選択できます。

さらには、モデルと予測期間のデータフレームをpredict関数に渡して、 plot関数などを使って予測結果を可視化します。

future <- make_future_dataframe(m, periods = 365, freq = "day")

forecast <- predict(m, future)

tail(forecast[c('ds', 'yhat', 'yhat_lower', 'yhat_upper')])

plot(m, forecast)

Sys.setlocale("LC_TIME", "en_US")

prophet_plot_components(m, forecast)

Calc_residuals(df_14day$y, forecast$yhat[1:nrow(df_14day)])

予測結果のプロット

青線の予測結果から、2週間ベースで、2022年末に4500人で、2023年明けに5000人を超えることが期待されますね。また、2023年6月ごろには6000人を目指せそう結果になってますね。

prophet_plot_componentsのプロット結果

トレンドは、2021年の後半から上昇トレンドにありますね。水曜日と金曜日に多い傾向にありますね。

他のプロット方法について

add_changepoints_to_plot関数で、潜在的変化点を重ね合わせできます。また、plot_forecast_componentでは、特定の成分をプロットします。

dyplot.prophetでは、dygraphsパッケージによるインタラクティブなプロットは可能です。

plot(m, forecast) + add_changepoints_to_plot(m)

plot_forecast_component(m, forecast, name="yhat")

plot_forecast_component(m, forecast, name="trend")

dyplot.prophet(m, forecast)

重要な変化点の重ね合わせプロット

特定の成分 yhat をプロットする

特定の成分 trend をプロットする

https://kumes.github.io/Blog/R_prophet/Fig_03_4.html

dyplot.prophet関数は、インタラクティブにプロットできて、使い勝手も良いですね。

過去28日間の訪問ユーザー数(y_28day列)を使って、ブログ利用ユーザー数の推移を未来予測する

次に、過去28日間の訪問ユーザー数(y_28day列)を使って、prophetを実行しました。

ここでは、少し実行内容を省略しています。

url <- "https://gist.githubusercontent.com/kumeS/c677d0b3ceb56f74cd109aae6acfd1a8/raw/6cf600582b8ba5bc3812d2a600de959c01548eca/Blog_Users_20200401-20221112.csv"

df <- read.csv(url, header = T)

df_28day <- data.frame(ds=as.Date(df$ds), y=as.numeric(sub(",", "", df$y_28day)), row.names = 1:nrow(df))

m <- prophet(df=df_28day, growth = "linear", n.changepoints = 100, changepoint.range = 1, changepoint.prior.scale=0.05, interval.width = 0.95, seasonality.mode = "multiplicative")

future <- make_future_dataframe(m, periods = 365, freq = "day")

forecast <- predict(m, future)

plot(m, forecast) + add_changepoints_to_plot(m)

dyplot.prophet(m, forecast)

Calc_residuals(df_28day$y, forecast$yhat[1:nrow(df_28day)])

重要な変化点の重ね合わせプロット

28日の利用者数の予測結果は、14日の利用者数の予測結果の2倍の値を示します。

来年1月初めに、利用者数が激減しているけども、何か対策はあるのかな。。

https://kumes.github.io/Blog/R_prophet/Fig_04_2.html

何方にせよ、2023年もアップ・トレンドを期待したいです。

https://www.macrotrends.net/2593/nikkei-225-index-historical-chart-data

最近、www.macrotrends.netで、日経225インデックスの67年間のチャートデータを見つけました。

ちょうどいいデータなので、今回使ってみることにしました。

まずは、Gistに成形した日経225インデックスのチャートデータをアップしていますので、それをダウンロードします。

url <- "https://gist.githubusercontent.com/kumeS/0b901711f87e5a6854a6b354559b11ff/raw/d018a8700d150d43a528a347ee0ffb73cee322c8/nikkei-225-index-67-year-historical-chart-data.csv"

df <- read.table(url, sep=",", header = TRUE)

df0 <- df rownames(df0) <- df0$ds dygraph(df0) %>% dyRangeSelector()

https://kumes.github.io/Blog/R_prophet/Fig_06.html

2012年以降のデータで、prophetで株価トレンドを予測する

日経225の2012年以降のデータを使用って、prophetで株価トレンドを予測してみました。

changepoint.prior.scaleは、潜在的変化点の自動選択の柔軟性を調整するパラメータです。大きな値を設定すると、多くの変化点が選択され、小さな値を設定すると、変化点は選択され難くなります。今回、そのスケール値を0.85に変えています。

また、make_future_dataframe関数では、24週間分のデータフレームを作成しています。

df2012 <- df[df$ds > "2012-01-01",]

head(df2012)

m <- prophet(df=df2012, growth = "linear", n.changepoints = 50, changepoint.range = 1, changepoint.prior.scale=0.85, interval.width = 0.8, seasonality.mode = "multiplicative")

future <- make_future_dataframe(m, periods = 24, freq="week")

forecast <- predict(m, future)

dyplot.prophet(m, forecast)

https://kumes.github.io/Blog/R_prophet/Fig_07.html

まぁ、フィットが悪いんだけど、年明けから日経225が下がったら面白いね。ふと思うと、アップトレンドをではないので、seasonality.mode = "additive"のほうがフィットは良いのかも・・・

ビットコイン(BTC)の価格推移を予測してみたところ

ここでは、getBTC関数を使って、 2010年7月18日から今日(現在、2022年11月20日)までのBTCの日足データを取得します。

それで、ビットコインの取得データは、2017年以降のデータを使用します。

source("https://gist.githubusercontent.com/kumeS/c466e97d42f6facbb8ba5a22ed361bf4/raw/07728e5e4810e53aa8bdcc411ac51603943761df/getBTC.R")

BTC <- getBTC() BTC0 <- BTC[,c("Date", "Close")] colnames(BTC0) <- c("ds", "y")

BTCdata <- BTC0[BTC0$ds > "2020-01-01",]

やっぱり、仮想通貨のように値動きの激しいデータは、パラメータ・チューニングに苦戦しますね。それでもかなりマシですけども。。。

アップトレンドで無いデータということで、seasonality.modeの引数はadditiveにしています。

m <- prophet(df=BTCdata, growth = "linear", n.changepoints = 250, changepoint.range = 1, changepoint.prior.scale=0.5, interval.width = 0.8, seasonality.mode = "additive")

future <- make_future_dataframe(m, periods = 30, freq="day")

forecast <- predict(m, future)

dyplot.prophet(m, forecast)

https://kumes.github.io/Blog/R_prophet/Fig_08.html

色々とパラメータ調整したからかも知れませんが、、この結果は、興味深いですね。

ビットコインが10000ドルあたりで反転するのか、どうかですね。

まとめ

現在、ブログのユーザー数は月で7000-8000人を推移してますが、今後どう推移してくるのか？、Prophetでの予想通りに増えてくるのか？、1万人の大台を超えてくるのか、楽しみですね。今後、実際のデータとの比較もできるとさらに面白そうだ。。。

あとは、弱気相場も大詰めで、ビットコインの底値がどのあたりになるかも気になりますね。時系列予測がどこまで使えるか、気になりますよね。

補足

Calc_residuals関数

このCalc_residuals関数では、 MSE(平均二乗誤差)、MAE(平均絶対誤差)、RMSE(二乗平均平方根誤差)、MAPE(平均絶対パーセント誤差)、およびR2(決定係数)を計算します。

Calc_residuals <- function(original, predicted, digits=4){

d0 <- original - predicted

d1 <- mean((d0)^2)

d2 <- mean(abs(d0))

d3 <- sqrt(mean((d0)^2))

d4a <- abs((d0/original)) d4b <- mean(d4a[!is.infinite(d4a)])*100

d5 <- round(1-(sum((d0)^2)/sum((original-mean(original))^2)), digits)

return(data.frame(mse=round(d1, digits=digits), mae=round(d2, digits=digits), rmse=round(d3, digits=digits), MAPE=round(d4b, digits=digits), R2=round(d5, digits=digits), row.names = 1)) }

prophet関数のオプション

引数	概要
df	(オプション) 履歴を含むデータフレーム。ds (日付型) と y (時系列) の列を持たなければならない。成長がロジスティックである場合、dfは、各dsの容量を指定する列capも持っていなければならない。モデルをフィットさせるにはfit.prophet(m, df)を使ってください。
growth	文字列(String) 'linear', 'logistic', または 'flat' で、線形、ロジスティック、またはフラットトレンドを指定します。
changepoints	潜在的なチェンジポイントを含む日付のベクトル。指定しない場合、変化点候補は自動的に選択されます。
n.changepoints	含めるべき潜在的な変化点の数。入力 'changepoints' が与えられている場合は使用されません。changepoints' が与えられない場合、 df$ds の最初の 'changepoint.range' の割合から一様に n.changepoints 潜在的な変化点が選択されます。
changepoint.range	トレンド変化点が推定される履歴の割合。デフォルトでは、最初の 80 'changepoints'が指定された場合、0.8 となります。
yearly.seasonality	年間の季節性をフィットさせる。auto', TRUE, FALSE, あるいは生成するフーリエ項の数を指定することができます。
weekly.seasonality	週ごとの季節性をフィットさせる。auto', TRUE, FALSE, あるいは生成するフーリエ項の数を指定することができます。
daily.seasonality	日次の季節性をフィットさせる。auto', TRUE, FALSE, あるいは生成するフーリエ項の数を指定することができます。
holidays	holiday (文字) および ds (日付型) の列と、休日として含める日付を指定する lower_window および upper_window の列を持つデータフレームです。また、オプションでprior_scaleカラムを持ち、各祝日の事前スケールを指定することができます。
seasonality.mode	'additive' (default) or 'multiplicative'.
seasonality.prior.scale	季節性モデルの強さを調節するパラメータ。大きな値を設定すると、より大きな季節変動に対応できるようになり、小さな値を設定すると、季節性が弱くなる。個々の季節性に対して add_seasonality を用いて指定することができる。
holidays.prior.scale	休日入力で上書きされない限り、休日成分モデルの強さを調節するパラメータ。
changepoint.prior.scale	チェンジポイントの自動選択の柔軟性を調整するパラメータ。大きな値を設定すると多くのチェンジポイントが選択され、小さな値を設定するとチェンジポイントはほとんど選択されません。
mcmc.samples	整数(Integer)、 0より大きい場合，指定された数のMCMCサンプルで完全なベイズ推定を行います．0の場合, MAP推定を行います.
interval.width	数値(Numeric)、予測に提供される不確実性区間の幅。mcmc.samples=0の場合、これは外挿された生成モデルのMAP推定値を使用したトレンドの不確実性のみとなる。mcmc.samples>0の場合、これは季節性の不確実性を含むすべてのモデルパラメータにわたって統合される。
uncertainty.samples	不確実性区間を推定するために使用されるシミュレーション描画の数。この値を0またはFalseに設定すると、不確実性推定が無効となり、計算が高速化されます。
fit	true/false, FALSE の場合、モデルは初期化されるが、フィットしない。

make_future_dataframe関数のオプション

引数	概要
m	Prophet モデルオブジェクト
periods	将来予測を行う期間数。
freq	'day', 'week', 'month', 'quarter', 'year', 1(1 sec), 60(1 minute) or 3600(1 hour).
include_history	true/false; 予測用のデータフレームに過去の日付を含めます。