時系列予測システムProphetを使った日経平均株価の予測能力を評価する

以前の記事においてProphetを使ってWikipediaのアクセス数日経平均株価の予測を行いました。
それらの記事では、モデルが実データに追従できていることを確認し、そのモデルから予測結果を導きました。
しかし予測するだけでは実務に役立てることはできません。
なぜならばその予測が100%正解であることはありえないからです。
だから我々は予測結果を使う前にその予測がどの程度正解に近いのか、すなわち予測精度はいかほどなのか、ということを合わせなければ恐ろしくて結果を使うことはできないのです。
よって我々に必要なのはProphetを使った予測性能の評価方法を確立することです。

【理論】

Prophet原論文で提示されている予測性能の評価方法を下敷きにします。
納得いかないところを変えていますので変数の表記が微妙に変わっていることをご容赦ください。

ある時点tからh日先を予測したいとします。
このとき予測したい実現値をy(h+t)と表します。
また予測に使えるデータが過去を遡ってT日分(すなわちt-Tからt-1の間)を使って求めた予測値を\hat y(h+t|T)と表します。
このとき予測誤差を変数t,T,hを用いて


\displaystyle e(t,T,h)=y(h+t)-\hat y(h+t|T)

と表すことができます。

さらに論文では相対誤差の絶対値を用いる方法が提示されています。よって誤差を


\displaystyle \epsilon(t,T,h)=\frac{|e(t,T,h)|}{y(h+t)}

として扱うこととします。

我々が知りたいのはh日先の予測性能です。
これを評価するためにまずTを固定します。
次にたくさんの異なる時点t_1,t_2,\ldots,t_Nから予測を行い、予測誤差\epsilon(t_i,T,h)を計算します。
最後に得られた予測誤差の平均値


\displaystyle \xi_T(h)=\frac{1}{n}\sum_{i=1}^N\epsilon(t_i,T,h)

をもって予測誤差の期待値とします。

これがProphet論文で提案されたSimulated Historical Forecasts(SHFs)と呼ばれる予測能力の評価手法です。

【実装】

ProphetライブラリにはSHFsは実装されていないため、自分で実装します。
使うデータは日経平均株価データです。

#R
library(dplyr)
library(prophet)
nikkei <- read.csv("data/nikkei_long.csv", header=TRUE)
nikkei <- nikkei %>% select(date, end) %>% rename(y = end, ds = date)  %>% arrange(ds)
SHF <- function(input, n.td = 365 * 3, max.horizon = 180, step = 90, log = FALSE){
        if(log){
                dataset <- input %>%
                  filter(y != 0) %>%
                  mutate(y = log(y))
        }else{
                dataset <- input
        }
        dataset <- dataset %>%
          mutate(ds = as.Date(ds)) %>%
          arrange(ds)

        cutoffs <- seq(min(dataset$ds) + n.td, max(dataset$ds), by = step)
        fcst <- data.frame()
        for(i in cutoffs){
                training.data <- dataset %>% filter(ds < i, ds >= i - n.td)
                m <- prophet(training.data)
                future <- make_future_dataframe(m, periods = max.horizon)
                fcst <- rbind(fcst,
                              predict(m, future) %>%
                                filter(ds >= i) %>%
                                mutate(h = seq(1, length(ds)), cutoff = min(ds))
                              )
        }

        if(log){
                fcst$yhat <- exp(fcst$yhat)
                dataset$y <- exp(dataset$y)
        }

        fcst <- fcst %>%
          inner_join(dataset, by = "ds") %>%
          select(ds, yhat, y, h, cutoff) %>%
          mutate(error = yhat - y) %>%
          mutate(relative.error = error/y)
        fcst.horizon <- fcst %>%
          group_by(h) %>%
          summarise(mape = mean(abs(relative.error)), vpe = var(relative.error))
        return(list(forecast = fcst, horizon = fcst.horizon))
}
 res <- SHF(nikkei, log = TRUE)

関数SHFinputに過去データ全体を突っ込むと、日数n.tdの学習データを使って1日先からmax.horizon日先までの予測を行います。
そしてそれらの予測誤差を計算したあと、現在のカレンダーをstep日将来にズラして、また同じように予測を行います。
首尾よくいけば、おおよそ

(inputの日数 – n.td) / step

回のサンプルを取ることができます。
これが論文で提示されているSHFです。
予測結果を図1に示します。

図1:SHFによる日経平均株価の予測性能。横軸hが何日予測するかの日数。縦軸は相対誤差の平均

横軸が何日将来を予測するか、縦軸が相対誤差の平均です。
これを見てまず分かるのは、予測日数45日あたりまでは予測日数に対して予測誤差の変動が非常に少ないということです。
これはProphetによる予測は45日程度までは安定した予測性能が期待できることを意味します。

次に予測日数100日を過ぎたあたりから誤差の成長が大きなクラスタと、誤差の成長が小さなクラスタとに分かれているということも分かります。
グラフをよく観察すると5日連続で上のクラスタ、2日下のクラスタ、といったように振動を繰り返しています。
この原因は極端に予測精度が悪いstepが1つあるためではないかと考えられます。
極端に予測精度が悪いstepが土日のときは実績値がなく予測が発生しないため誤差も発生せず、周期的に予測精度が改善するのです。

#R
library(ggplot2)
ggplot(filter(res$horizon, h %in% 150:170), aes(h, mape)) + geom_col()
図2:SHFの誤差を100以上の予測日数に拡大。5日、2日の周期で誤差が大きくなったり小さくなったりする

【SHFの改善(randomized SHF)】

なぜ図2のような問題が起こるかというと、異なる予測の予測誤差の間に相関があるからです。
なぜ誤差間に相関があるかというと、SHFではprophetの予測期間をmax.horizon(180日)で設定し、予測日数1から180までの180個の予測を一つの予測モデルで行うからです。
そしてなぜそのように一つの予測モデルで180回の予測を行うのかといいますと、計算時間を節約したいからです。

説明します。
Prophetの予測は予測モデルを作って予測するのには少し計算時間を必要とします。
しかし一回の予測で複数の日数の予測を行うのに計算時間はほとんど増加しません。
よって一つの予測モデルでたくさん予測したほうが計算時間を節約することができるのです。

その見返りとして上で起こったように誤差間の相関が発生してしまいます。
これを避けるためには予測のたびに予測モデルを生成して予測することです。
これをrandomized SHF(rSHF)と呼ぶことにしましょう。
rSHFは普通のSHFよりもかなり計算時間がかかります。

#R
fcst_fixed_horizon <- function(input, horizon = 50, n.trial = 2, n.training = 365 * 3, step = 10){
        if(nrow(input) - n.training < n.trial){
                print("Error: the observations of input data is smaller than training data.")
                return(NULL)
        }

        min.date <- min(input$ds) + n.training + 1
        max.date <- max(input$ds) - horizon
        learning.seq <- sample(input$ds[input$ds > min.date & input$ds < max.date], n.trial)
        out <- data.frame()
        for(di in learning.seq){
                td <- input %>% filter(ds < di, ds >= di - n.training)
                m <- prophet(td)
                future <- make_future_dataframe(m, periods = horizon)
                forecast <- predict(m, future) %>% tail(horizon) %>% mutate(h = seq(1, horizon))
                out <- rbind(out, tail(forecast, step)) #予測した日に対応する実績が無いかもしれないので幅をもたせる
        }
        return(out)
}
rSHF <- function(input, step = 7, min.fcst = 1, max.fcst = 180, n.trial = 2, n.training = 365 * 3, log = FALSE){
        if(log){
                dataset <- input %>%
                  filter(y != 0) %>%
                  mutate(y = log(y))
        }else{
                dataset <- input %>% mutate(ds = as.Date(ds))
        }
        dataset <- dataset %>%
          mutate(ds = as.Date(ds)) %>%
          arrange(ds)

        fcst <- data.frame()
        for(i in seq(min.fcst, max.fcst, by = step)){
                fcst <- rbind(fcst, fcst_fixed_horizon(dataset, horizon = i, n.trial = n.trial, n.training = n.training, step = step))
        }

        if(log){
                fcst$yhat <- exp(fcst$yhat)
                dataset$y <- exp(dataset$y)
        }
        fcst <- fcst %>%
          inner_join(dataset, by = "ds") %>%
          select(ds, yhat, y, h) %>%
          mutate(error = yhat - y) %>%
          mutate(relative.error = error/y)
        fcst.horizon <- fcst %>%
          group_by(h) %>%
          summarise(mape = mean(abs(relative.error)), vpe = var(relative.error))

        return(list(forecast = fcst, horizon = fcst.horizon))
}
res <- rSHF(nikkei, log = TRUE)
ggplot(res$horizon, aes(h, mape)) + geom_point()
図3:randomized SHFによる日経平均株価の予測性能。横軸hが何日予測するかの日数。縦軸は相対誤差。SHFと比較して周期性が無くなった

ランダムサンプリングの結果、hに対する誤差の周期性が観察できなくなりました。
さらに図1で観察されていた予測日数45日以降で予測誤差が安定しなくなる、すなわち予測が悪化するという現象が顕著に確認できるようになりました。
この図3では60日くらいまでは予測が安定していると言えそうです。

【結論】

Prophetを使って日経平均株価を予測したとき、50日程度までは予測することができる。
そのときの50日後の期待予測誤差は10%程度である。

免責事項:このページで著者が主張していることは主観を多く含んでいるため、主張内容を活用したとしても利益を保証するものではありません。

関連記事

時系列予測システムProphetを使った日経平均株価の予測能力を評価する” への3件のコメント

  1. 素人質問ですみません。
    予測する際、prophetで95%信頼区間と思われる薄い青の区間の中に濃い青の線が存在しています。初めは、単純な平均だと思っていたのですがどうやらそうでは、無いようです。株価で大きな転換点を迎える場合などは、異様に太くなります。また、株価で相対誤差の1%の値はどのように計算すれば良いのでしょうか。
    長くなってしまいました。もしよろしければお教えください。

    1. ご質問ありがとうございます。

      >予測する際、prophetで95%信頼区間と思われる薄い青の区間の中に
      まず薄い青の区間について、これは確かに信頼区間なのですが、伝統的な統計学が言うところのいわゆる95%信頼区間ではありません。
      なぜならprophetは統計モデルとしてベイズモデルを用いており、予測も解析的に解くのではなくシミュレーションで行っているため、古典的な意味での信頼区間というものは得られないからです。
      よって薄い青色の帯を直感的に解釈するならば、モデルにしたがってシミュレーションした結果、得られたサンプルの一定割合が青色の帯の中に収まったということです。
      なお、prophetでその割合はデフォルトでは80%に設定されており、私のデモンストレーションもすべて80%の結果です。
      80%という数値は関数prophetの引数interval.widthで設定できます。

      >濃い青の線
      基本的にはシミュレーションで得られたサンプルの平均値を推定値であり、濃い青の線だと思います。
      シミュレーションの結果であり、また古典的な統計のように上下に等しく分布するとは限らないため、薄い青の線の平均値とは必ずしも一致しません。

      不安定な株価の動きを反映して信頼区間が異様に太くなったのならば推定値は信用するべきではないでしょうね。
      なお、上記の話はまだ内部仕様を確認したわけではないので推測です。
      いつか信頼区間についても記事を書きたいと思います。

      >株価で相対誤差の1%の値はどのように計算すれば良いのでしょうか
      相対誤差は予測誤差を株価そのもので割り算することで計算できます。
      たとえば予測誤差が100円でその時の株価が10000円ならば1%の誤差です。
      本記事の【理論】の二つ目の数式を参照ください。
      ご質問の意図とズレていたらすいません。

  2. ご回答ありがとうございました。
    大変勉強になりました!

コメントを残す