この記事は、以下の@icoxfog417さんによる問題提起に合わせたちょっとした実験をまとめたものです。
時系列予測の問題において、機械学習のモデルより既存の統計モデル(ARMAモデルなど)の方が予測精度において優良な結果が出るという研究。データへの適合=予測精度の向上ではないことも実験で示している。機械学習の研究では統計モデルとの比較も入れるべきという提言をしている。 https://t.co/jboGhYSX6E
— piqcy (@icoxfog417) September 16, 2019
この点について僕はこんなコメントをしたのですが。
だいぶ以前から「一般的な時系列データ予測の問題は単位根過程や季節調整など非定常過程との戦いなので、本質的に定常過程を想定する機械学習手法での予測は計量時系列分析など非定常過程も考慮した古典的なモデルによる予測には及ばない」と言い続けてきたけど、やっぱりそうなのかなという感想
— TJO (@TJO_datasci) September 16, 2019
この手の「一般的な時系列データには機械学習は向かない」論は過去何度か断片的にしてきてはいるものの、そう言えばまとめて論じたことはなかったなと思ったのでした。そこで、今回の記事では何故一般的な時系列データを機械学習で扱うのが難しいのかについて、簡単にまとめてみようと思います。
なお、いつもながらですが記事中で用いているR / Stan実装が途方も無いクソコードである点何卒ご容赦くださいm(_ _)m 加えて、理論的なポイントについてもやはり抜け漏れや理解不足などあるかと思いますので、お気付きの方はバンバンご指摘下さると有難いです(多分もっとしっかりとした確率過程の理論的背景なども踏まえた論文とかどこかにあるのではないかと思うのですが、僕には手が出ないレベルなのでここでは割愛します)。
また、僕自身はDeep Learningでもよく時系列データに対して使われるRNN系統の手法は不案内で、正直言って適切な実装が出来る自信がありませんので、今回はDeep系の手法は割愛し、代わりに枯れた理論&実装で知られるランダムフォレストを機械学習サイドの代表例として取り上げることとしました。ただしもしかしたら以下に指摘した問題点に対してロバストに良い結果を返すRNN系統の手法があるかもしれませんので、予めお断りしておきます。
- データの準備
- ベイジアン動的線形モデルとランダムフォレストでそれぞれモデリングして予測(forecast)してみる
- ランダムフォレストに曜日の情報を与えてやり直してみる
- 機械学習系の手法は定常過程を前提としていて、非定常過程に弱い
- 余談
- 追記
- その1
- その2
- その3
データの準備
共変量(説明変数)付きの時系列データセット"sample_marketing.csv"を僕のGitHubに用意してあるので、こちらをDLしておいてください。
想定としてはいわゆるMedia Mix Model (MMM)を模したもので、1列目に目的変数である"cv"、2列目以降に広告出稿ボリュームを表す"ad1"以下3つの説明変数が入っています。適当にプロットしてみるとこうなります。
par(mfrow = c(4, 1)) par(mar = c(1, 4, 1, 2) + 0.1) plot(d$cv, type = 'l', lwd = 2, col = 2) plot(d$ad1, type = 'l', lwd = 1, col = 1) plot(d$ad2, type = 'l', lwd = 1, col = 1) plot(d$ad3, type = 'l', lwd = 1, col = 1)
cvは見ているとちょっとうねった感じのトレンドと若干の周期性があるようにも見えます。一方、adは期間によっては値が入っていたり入っていなかったりするようです。
ベイジアン動的線形モデルとランダムフォレストでそれぞれモデリングして予測(forecast)してみる
この手の時系列回帰モデル(動的線形モデル)を組む際に考慮することは、基本的には以下の3パートに分けるという点です。
- 共変量(説明変数)による回帰パート
- 季節調整パート
- 長期トレンドパート
1つ目は言うまでもなく回帰モデルです。2つ目は24時間・7日(曜日)・12ヶ月(年)といった主にカレンダーに依存する周期変動です。そして3つ目は単位根過程やローカルレベルトレンドや二階差分トレンドなど、各時点間で関連しながら非線形に発展していくトレンドです。今回は、季節調整には7日周期成分を、トレンドには二階差分トレンドを当てはめてみることにします。つまり以下のようなモデルです。
式だけ見ていると面倒な印象がありますが、要は7日周期成分は「7日間の和が平均0&標準偏差の正規分布に従う」、二階差分トレンドは「現時点と前の時点とのトレンド成分の差が平均0&標準偏差の正規分布に従う」、ということです。
まず、ベイジアン動的線形モデルを推定してみます。方法は幾らでもありますが、古典的な時系列モデルらしさを味わうためにあえてStanで組んでみます。StanファイルをGitHubに上げてありますが、以下からコピペしても大丈夫です。
data { int<lower=0> N; int<lower=0> M; matrix[N, M] X; real<lower=0> y[N]; } parameters { real season[N]; real trend[N]; real s_trend; real s_season; real s_q; vector<lower=0>[M] beta; } model { real q[N]; real cum_trend[N]; for (i in 7:N) season[i]~normal(-season[i-1]-season[i-2]-season[i-3]-season[i-4]-season[i-5]-season[i-6],s_season); for (i in 3:N) trend[i]~normal(2*trend[i-1]-trend[i-2],s_trend); cum_trend[1]=trend[1]; for (i in 2:N) cum_trend[i]=cum_trend[i-1]+trend[i]; for (i in 1:N) q[i]=y[i]-season[i]-cum_trend[i]; for (i in 1:N) q[i]~normal(dot_product(X[i], beta),s_q); }
これをkickするRコードもGitHubに上げてあります。ご覧になればお分かりかと思いますが、100点のデータのうち91点を学習データ&9点をテストデータとしてholdout splitしてあります。つまり、91点でモデルを学習させて後ろの9点のcvを予測(forecast)*1させてその精度を見る、という検証のやり方です。
またこのRコードでは、同時にランダムフォレストで全く同じ学習&テストデータに対する機械学習モデル推定&予測を行なっています。
とりあえず以下のところまで回してみましょう。
# Data preparation d <- read.csv('sample_marketing.csv') d_train <- d[1:91, ] d_test <- d[92:100, ] dvar <- d_train[, -1] dat <- list(N = nrow(d_train), M = ncol(dvar), y = d_train$cv, X = dvar) # Time series modeling with Bayesian library(rstan) rstan_options(auto_write = TRUE) options(mc.cores = parallel::detectCores()) fit <- stan(file = 'sample_marketing_holdout.stan', data = dat, iter = 1000, chains = 4) # Extract samples / parameters, and estimate fitted values slength <- 2000 fit.smp <- extract(fit) tmp <- density(fit.smp$s_trend) s_trend <- tmp$x[tmp$y == max(tmp$y)] tmp <- density(fit.smp$s_season) s_season <- tmp$x[tmp$y == max(tmp$y)] tmp <- density(fit.smp$s_q) s_q <- tmp$x[tmp$y == max(tmp$y)] beta <- rep(0, ncol(dvar)) for (i in 1:ncol(dvar)) { tmp <- density(fit.smp$beta[(slength * (i - 1) + 1) : (slength * i)]) beta[i] <- tmp$x[tmp$y == max(tmp$y)] } trend <- rep(0, nrow(dvar)) for (i in 1:nrow(dvar)) { tmp <- density(fit.smp$trend[, i]) trend[i] <- tmp$x[tmp$y == max(tmp$y)] } season <- rep(0, nrow(dvar)) for (i in 1:nrow(dvar)) { tmp <- density(fit.smp$season[,i]) season[i] <- tmp$x[tmp$y == max(tmp$y)] } beta_prod<-rep(0, nrow(dvar)) for (i in 1:ncol(dvar)) { beta_prod <- beta_prod + dvar[, i] * beta[i] } pred <- beta_prod + cumsum(trend) + season matplot(cbind(d_train$cv, pred), type = 'l', lty = 1, lwd = c(2, 3), ylab = '', xlim = c(0, 100), ylim = c(180, 980)) legend('topleft', legend = c('Data', 'Fitted'), lty = 1, lwd = c(2, 3), col = c(1, 2)) # Validation fc_dvar <- d_test[, -1] val_trend <- rep(0, nrow(fc_dvar)) for (i in 1:2) { val_trend[i] <- trend[nrow(d_train) - 3 + i] } for (i in 3:nrow(fc_dvar)){ val_trend[i] <- 2 * val_trend[i - 1] - val_trend[i - 2] + rnorm(1, 0, s_trend) } val_season <- season[(length(season) - 14 + 1):(length(season) - 14 + nrow(fc_dvar))] + rnorm(nrow(fc_dvar), 0, s_season) fc_beta_prod <- rep(0, nrow(fc_dvar)) for (i in 1:ncol(fc_dvar)) { fc_beta_prod <- fc_beta_prod + fc_dvar[, i] * beta[i] } validate <- cumsum(trend)[nrow(d_train)] + fc_beta_prod + cumsum(val_trend) + val_season # Machine learning: Random Forest library(randomForest) d_train.rf <- randomForest(cv~., d_train) pred.rf <- predict(d_train.rf, newdata = d_train[, -1]) validate.rf <- predict(d_train.rf, newdata = d_test[, -1]) # Plotting w/o days matplot(cbind(d_train$cv, pred, pred.rf), type = 'l', lty = 1, lwd = c(2, 3, 3), ylab = '', xlim = c(0, 100), ylim = c(150, 1050)) segments(91, 0, 91, 1100, lty = 3, col = 1, lwd = 2) par(new = T) matplot(cbind(d_test$cv, validate, validate.rf), type = 'l', lty = 1, col = c('blue', 'purple', '#008000'), lwd = c(2, 3, 3), ylab = '', xlim = c(-91, 9), ylim = c(150, 1050), axes = F) legend('topleft', legend = c('Data', 'Fitted TS', 'Fitted RF', 'Actual', 'Forecast TS', 'Forecast RF'), lty = 1, lwd = rep(c(2, 3, 3), 2), col = c(1, 2, 3, 'purple', 'blue', '#008000'), ncol = 2, cex = 0.75)
ベイジアン動的線形モデルでは綺麗に学習データにフィットした上に予測結果もテストデータに割合良く合っているのに対して、ランダムフォレストはべったりと水平に伸びるだけの有様になってしまっています。それもそのはずで、前者は季節調整を初めから含んでいる*2のに対して、ランダムフォレストにはその情報は与えられていないからです。
ランダムフォレストに曜日の情報を与えてやり直してみる
上記の結果のままだと、ランダムフォレストに曜日の情報が与えられていないことになって不公平です。ということで、曜日を表すダミー変数を追加したデータでランダムフォレストだけやり直してみましょう。
# RF w/ days d_day <- cbind(d, rep(c('Mon', 'Tue', 'Wed', 'Thr', 'Fri', 'Sat', 'Sun'), 15)[1:nrow(d)]) names(d_day)[5] <- 'day' day_train <- d_day[1:91, ] day_test <- d_day[92:100, ] day_train.rf <- randomForest(cv~., day_train) pred.rfday <- predict(day_train.rf, newdata = day_train[, -1]) validate.rfday <- predict(day_train.rf, newdata = day_test[, -1]) # Plotting w/ days matplot(cbind(d_train$cv, pred, pred.rfday), type = 'l', lty = 1, lwd = c(2, 3, 3), ylab = '', xlim = c(0, 100), ylim = c(150, 1050)) segments(91, 0, 91, 1100, lty = 3, col = 1, lwd = 2) par(new = T) matplot(cbind(d_test$cv, validate, validate.rfday), type = 'l', lty = 1, col = c('blue', 'purple', '#008000'), lwd = c(2, 3, 3), ylab = '', xlim = c(-91, 9), ylim = c(150, 1050), axes = F) legend('topleft', legend = c('Data', 'Fitted TS', 'Fitted RF w day', 'Actual', 'Forecast TS', 'Forecast RF w day'), lty = 1, lwd = rep(c(2, 3, 3), 2), col = c(1, 2, 3, 'blue', 'purple', '#008000'), ncol = 2, cex = 0.75)
確かに曜日変動=季節調整の部分は反映されるようになりましたが、当然ながらトレンドの部分は反映されておらず*3、相変わらず平坦なフィット&予測のままです。最後のholdout部分のforecastもベイジアン動的線形モデルに比べると少し外れ方が大きいようです。そこでRMSEを計算してみましょう。
library(Metrics) rmse(d_test$cv, validate) # Bayesian TS rmse(d_test$cv, validate.rf) # RF w/o days rmse(d_test$cv, validate.rfday) # RF w/ days
乱数シード依存ではありますが、RMSEを比べてみると70.7308, 88.6625, 88.78605で、やはりベイジアン動的線形モデルが一番良いという結果になります。
機械学習系の手法は定常過程を前提としていて、非定常過程に弱い
では、何故こんなことになるのでしょうか? 例えば、オリジナルデータをholdout splitする前の目的変数cvの前半50点と後半50点を見比べるとこうなっています。
hist(d$cv[1:50], col = 'red', xlim = c(0, 1000), ylim = c(0, 20), main = '', xlab = '') par(new = T) hist(d$cv[51:100], col = 'blue', xlim = c(0, 1000), ylim = c(0, 20), main = '1:50 (red) vs. 51:100 (blue)', xlab = '')
明らかに目的変数の分布が大きく上下にスライドしていることが分かります。このようなデータを機械学習系の手法で取り扱うのは困難であろうと思われます。
これは、6年も前の話ですがこのブログで計量時系列分析について一連のシリーズ記事を書いた際に、一番最初に挙げていたポイントを思い起こせば容易に理解できると思います。ちょっと引用してみましょう。
まず、弱定常性は過程の期待値と自己共分散が時間を通じて一定である、というものです。これは想像しやすいと思います。
任意のtとkに対して、
が成立する場合、過程は弱定常といわれる
一方、強定常性は同時分布が不変である、というものです。
任意のtとkに対して、の同時分布が同一である場合、過程は強定常といわれる
基本的には、機械学習系の回帰モデルで扱うデータは最低でも弱定常、うっかりすると強定常を仮定しているのではないでしょうか。即ち「学習データもテストデータも原則として同一の分布から生成されたものと仮定する」というわけです。これはランダムフォレストのような枯れた手法でも、Deep Learning系の〇〇Netのような最新手法でも、大きく変わらないはずです。実は5年前に出版した拙著でも、以下のような図でこの概念についての説明を試みたことがあります。
要は「機械学習は『学習データとテストデータを通じて全体として目的変数の分布が不変である(もしくは変化があっても既知である)』ことを前提としている」のですが、これに対して強い自己相関(系列相関)=顕著なトレンドを伴う時系列データだとそもそも目的変数の分布自体が不変どころかどんどんスライドしていってしまうので、どのようにモデルを推定しても全く合わないというわけです。
そこで、上のプロットにベイジアンモデリングで推定した長期トレンドをオレンジ色で重ね書きしてみるとこうなります。
matplot(cbind(d_train$cv, pred, pred.rfday, cumsum(trend)), type = 'l', lty = 1, lwd = c(2, 3, 3, 5), ylab = '', xlim = c(0, 100), ylim = c(150, 1050), col = c(1, 2, 3, 'orange')) segments(91, 0, 91, 1100, lty = 3, col = 1, lwd = 2) par(new = T) matplot(cbind(d_test$cv, validate, validate.rfday, cumsum(val_trend) + cumsum(trend)[nrow(d_train)]), type = 'l', lty = 1, col = c('blue', 'purple', '#008000', 'orange'), lwd = c(2, 3, 3, 5), ylab = '', xlim = c(-91, 9), ylim = c(150, 1050), axes = F) legend('topleft', legend = c('Data', 'Fitted TS', 'Fitted RF w day', 'Actual', 'Forecast TS', 'Forecast RF w day', 'Trend'), lty = 1, lwd = rep(c(2, 3, 3, 5), 2), col = c(1, 2, 3, 'blue', 'purple', '#008000', 'orange'), ncol = 2, cex = 0.75)
当然ですが、ランダムフォレストはこの長期トレンドを捉えることが全く出来ていません。このようなデータに対して、機械学習系の手法ではうまく扱えない理由はこれで大体説明がつくと考えられます。
ちなみに元の指摘では「ARMA*4のような定常過程モデルであっても機械学習系の手法よりも時系列予測では有利」という趣旨のことが言われているようですが、正直これは僕の浅学では良く分かりませんでした。もしかしたら、これも結局は時系列データの中で確率的に振る舞う部分をどれだけ織り込めるかという定式化の部分の差なのかなと思ったのですが、確信がないので余計なコメントは差し控えることとします。
とは言え、この辺の事情がつまるところ「株価を機械学習で予想させようとしても上手くいかない」理由でもあることは間違いないと思います。ちなみに単一の株価を単変量モデリングして予測するのは計量時系列モデルでも基本的には無理、というか出来ていたら考案した人は既に大金持ちになっているはずだと付記しておきます。。。
余談
ところで、@hayashiyusさんからは以下のようなコメントをいただきました。
系列データの生成モデルでは.こうした機械学習モデルの弱点を巧く取り扱っているものが多いと思う.こうしたモデルでは,①変数の時間発展は古典的な予測モデルでモデリングする,②同一時点における識別タスクは機械学習モデルで解く.
— Yusuke HAYASHI 林 祐輔 (@hayashiyus) September 17, 2019
1. https://t.co/9KpNG5GESs
2. https://t.co/31aFsBIyIh
古典的な計量時系列モデルと機械学習モデルとのハイブリッドにするというのは、良案だなと思った次第です。
最後に、今回の僕のR / Stan実装はかなり酷い有様なので、きちんと学びたいという方は以前の書評記事でも絶賛したこちらの一冊をぜひご参考になさってください。
実践Data Scienceシリーズ RとStanではじめる ベイズ統計モデリングによるデータ分析入門
- 作者: 馬場真哉
- 出版社/メーカー: 講談社
- 発売日: 2019/07/10
- メディア: 単行本
- この商品を含むブログを見る
追記
その1
RNN(LSTM)はARIMAの要素を包含しているのではないでしょうか?もしそうならば、単純に機械学習が有効になるだけのデータ不足ということはないでしょうか?
— suama (@suamablogger) September 18, 2019
直感的に考えればLSTMでいけるのではないかという気は実は僕もしているのですが、その割にLSTMで社会・ファイナンス時系列データの予測がうまくいったという話をなかなか聞かないので、詳しい方からのコメントがあるまでは一旦保留とします。
ただし、「トレンドを表す特徴量を加えれば良い」というのは解決策になりにくいように思います。何故なら、学習データには目的変数の増加率などの名目でトレンドを表す説明変数を与えることが出来ますが、そもそも目的変数が未知であるテストデータにはそのようなものは与えようがないからです。これが(あくまでも可能だというレベルですが)計量時系列モデルであれば出来る、という側面を今回の実験では強調しています。
その2
トレンドへの対処として、「見せかけの回帰」を回避するのと同様に目的変数の一階差分を取って差分過程に直すという考え方もあります。
とはいえ、差分過程はそもそも様々な情報が「落ちて」しまうのと、また共和分関係にある変数がある場合にはこれはこれで(VARモデルに限らず)機械学習モデルに悪さをする可能性がなきにしもあらずという気もするので、ベストな解決策かというと確信がないです。。。