先日書いたこの記事ですが、「トイデータとは言え乱数シードを一つに決めて発生させたランダムウォークに対して実験をしているので、乱数シードを複数通りに変えてみたら結果は変わってくる（再現しない）のではないか？」という指摘を何人かの友人知人から貰いました。正直言って多項式フィッティングには何の思い入れもないのですが、再現性があるかどうかについては単純に気になるところです。
ちなみに、以前沖本本で勉強した際にあった「ランダムウォークには平均回帰性がなく時間と共に不確実性が増すため事実上予測不可能（特に長期予測）」という議論の通りで、本来ならランダムウォークに対して「予測」を行うのはそもそも適切ではありません。ただし、短期予測なら例えば状態空間モデルやBSTS的な方法で多少は精度を改善できるのでは？と考えていたのは事実で、同じことが無根拠でナンセンスな多項式フィッティングでも出来たら面白いかもと思ったというのも先日の記事の狙いの一つでした。
ということで、物は試しということで実験をやってみました。やり方は先日の記事から少し変えていますが大体同じで、
- 乱数シード101から300までの200通りで生成したランダムウォークで実験
- ランダムウォークの長さ（サイズ）は360に固定
- データを359:1に分け、前者をtrain + val、後者をtestとする
- さらにtrain:valは80%:20%とする
- 多項式フィッティングはデフォルト20次で、CVする場合はvalでベストスコアを出した次数のモデルをtrainに対して推定したものを用いてtest期間の値を予測する
- 状態空間モデルによるフィッティングは{bsts}でローカルレベルモデルとローカル線形トレンドモデルをそれぞれ別にtrain + valで推定し、直接test期間の値を短期予測する（予測値はtest期間の各時点で得られた予測分布のmedianをとる）
- 20次多項式・CV付き多項式・BSTSローカルレベル・BSTSローカル線形トレンドの4通りのモデルに対してRMSEを乱数シードごとに記録する（計200個）
- 最後に箱ひげ図でそれぞれのRMSEの分布を比較し、必要に応じてt検定及び順位和検定を行う（RMSEなので小さいほど予測性能に優れている）
ということをやっています。詳細はRコードをご覧いただいた方が早いでしょう。クソコードなので最適でもない上にかなり遅いですが、ご容赦ください。{bsts}についてはill-identifiedさんの記事が最も分かりやすいと思います。ただし{bsts}を選んだ理由は単純で、{dlm}の使い方は忘れてしまった上に{KFAS}は不慣れでよく分からないからです。。。
なお、コード冒頭で大半の設定は変えられるようになっていて、例えば20次多項式だと出力が上下に暴れて不利だというのであれば上限6次ぐらいに制限することも可能です。当然ながらtest期間長はおろか全体のデータの長さを変えることもできます。実験のセットアップやコードに不具合や誤りなどあれば是非ご指摘ください。
library(Metrics) # Just for computing RMSE library(bsts) # Call BSTS n <- 360 # Length of each random walk r_mean <- 0 # Mean of each random walk r_sd <- 20 # SD of each random walk r_int <- 1000 # Intercept of each random walk x <- 1:n # x vector: a basic component of independent variable l_range <- 1:359 # Range of learning dataset t_range <- 360:360 # Range of test dataset dcap <- 20 # Maximum limit of degree: default is 20 q_range <- 101:300 # Range of random seeds: of course you can change this range cv <- rep(0, length(q_range)) # Vector for storing results of polynomial fitting with CV ncv <- rep(0, length(q_range)) # Vector for storing results of polynomial fitting without CV bsll <- rep(0, length(q_range)) # Vector for storing results of bsts with local level bslt <- rep(0, length(q_range)) # Vector for storing results of bsts with local linear trend # Loop for each random seed for (q in q_range){ set.seed(q) # Set a random seed print(paste0(q, 'th loop')) y <- cumsum(rnorm(n, r_mean, r_sd)) + r_int # Generate a random walk # Set up a data frame for each total dataset d <- data.frame(y = y, x1 = x, x2 = x^2, x3 = x^3, x4 = x^4, x5 = x^5, x6 = x^6, x7 = x^7, x8 = x^8, x9 = x^9, x10 = x^10, x11 = x^11, x12 = x^12, x13 = x^13, x14 = x^14, x15 = x^15, x16 = x^16, x17 = x^17, x18 = x^18, x19 = x^19, x20 = x^20 ) d_learn <- d[l_range, ] # Learn dataset, 1:300 d_train <- d_learn[1:(round(0.8 * nrow(d_learn))),] # Training d_val <- d_learn[(round(0.8 * nrow(d_learn)) + 1):nrow(d_learn),] # Validation # CV procedure res <- rep(0, 20) # 20 steps CV for (j in 1:20){ # Polynomial fitting with increasing the degree fit <- lm(y ~., d_train[, 1:(1 + j)]) res[j] <- rmse(d_val$y, predict(fit, d_val[, 1:(1 + j)])) # Record RMSE } idx <- which(res == min(res)) # Choose the best one with the least RMSE fit_best <- lm(y ~ ., d_train[, 1:(2 + idx)]) # Best polynomial fitting model fit_all <- lm(y ~ ., d_train[, 1:(dcap + 1)]) # Full polynomial fitting model # Fit BSTS model with local level + random walk ss1 <- AddLocalLevel(list(), d_learn$y) fit_bs_ll <- bsts(d_learn$y, ss1, niter = 500) # Fit BSTS model with local linear trend ss2 <- AddLocalLevel(list(), d_learn$y) ss2 <- AddLocalLinearTrend(ss2, d_learn$y) fit_bs_lt <- bsts(d_learn$y, ss2, niter = 500) print(paste0(q, 'th loop')) newdata <- d[t_range, ] # Test dataset, 301:360 # Prediction by the best model pred_best <- predict(fit_best, newdata = newdata[, 1:(2 + idx)]) # Prediction by the full model pred_all <- predict(fit_all, newdata = newdata[, 1:(dcap + 1)]) # Prediction by BSTS model with local level + random walk pred_bs_ll <- predict(fit_bs_ll, horizon = length(t_range), burn = 100) # Prediction by BSTS model with local linear trend pred_bs_lt <- predict(fit_bs_lt, horizon = length(t_range), burn = 100) # Store results for each random seed cv[q - q_range[1] + 1] <- rmse(pred_best, newdata$y) ncv[q - q_range[1] + 1] <- rmse(pred_all, newdata$y) bsll[q - q_range[1] + 1] <- rmse(pred_bs_ll$median, newdata$y) bslt[q - q_range[1] + 1] <- rmse(pred_bs_lt$median, newdata$y) } paste0('cv_', dcap, 'deg_polynom_bsts_random_seeds_', q_range[1], '_', q_range[length(q_range)], '.RData') save.image(paste0('cv_', dcap, 'deg_polynom_bsts_random_seeds_', q_range[1], '_', q_range[length(q_range)], '.RData')) # Final evaluation # Boxplot boxplot(cv, ncv, bsll, bslt, names = c('With CV', 'Without CV', 'BSTS local level', 'BSTS local linear trend')) boxplot(cv, bsll, bslt, names = c('With CV', 'BSTS local level', 'BSTS local linear trend')) t.test(cv, ncv) # Welch's two-sample t-test wilcox.test(cv, ncv) # Ranksum test t.test(cv, bsll) t.test(cv, bslt)
今回は完全に数値実験なのですが、どんな感じのモデルフィッティングが行われているかを見るために、適当に乱数シードを4つ抜き出してきてtest期間長を60にした時の結果をプロットしたものを置いておきます。{bsts}でもtestへの予測があまり良くないケースもあることや、最大次数の多項式フィッティングの予測は上下に大きくブレ過ぎてプロット範囲内に描画されないケースがあることが、何となく見て取れます。
では、実験結果を見てみましょう。{bsts}のせいでそこそこ重いコードなので、結構時間がかかると思います。20次を多項式フィッティングの上限次数として、回してみた結果が以下の通りです。
> t.test(cv, ncv) # Welch's two-sample t-test Welch Two Sample t-test data: cv and ncv t = -17.949, df = 199, p-value < 2.2e-16 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: -55175826 -44252124 sample estimates: mean of x mean of y 1100.267 49715075.293 > wilcox.test(cv, ncv) # Ranksum test Wilcoxon rank sum test with continuity correction data: cv and ncv W = 0, p-value < 2.2e-16 alternative hypothesis: true location shift is not equal to 0 > t.test(cv, bsll) Welch Two Sample t-test data: cv and bsll t = 4.1835, df = 199, p-value = 4.303e-05 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: 572.8268 1594.3485 sample estimates: mean of x mean of y 1100.26702 16.67936 > t.test(cv, bslt) Welch Two Sample t-test data: cv and bslt t = 4.1825, df = 199, p-value = 4.32e-05 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: 572.5687 1594.0908 sample estimates: mean of x mean of y 1100.2670 16.9373
ちなみに多項式フィッティングの上限を6次に押さえた場合の結果は以下の通りです。
> t.test(cv, ncv) # Welch's two-sample t-test Welch Two Sample t-test data: cv and ncv t = -10.112, df = 389.23, p-value < 2.2e-16 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: -4125.486 -2782.425 sample estimates: mean of x mean of y 1100.267 4554.223 > wilcox.test(cv, ncv) # Ranksum test Wilcoxon rank sum test with continuity correction data: cv and ncv W = 3632.5, p-value < 2.2e-16 alternative hypothesis: true location shift is not equal to 0 > t.test(cv, bsll) Welch Two Sample t-test data: cv and bsll t = 4.183, df = 199, p-value = 4.313e-05 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: 572.6792 1594.2011 sample estimates: mean of x mean of y 1100.26702 16.82688 > t.test(cv, bslt) Welch Two Sample t-test data: cv and bslt t = 4.1824, df = 199, p-value = 4.323e-05 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: 572.5263 1594.0483 sample estimates: mean of x mean of y 1100.26702 16.97971
当たり前ですが、{bsts}で推定したローカルレベルorローカル線形トレンド状態空間モデルによる予測性能はいずれの多項式フィッティングよりも上回っています（RMSEがより小さくなる）。が、上限20次でも6次でも多項式フィッティングはCVした方がしないよりも予測性能で有意に上回るようです。これの解釈について数理統計が専門の友人の一人にコメントを求めたところ「CVによって直近のところがより当たっていて直近の予測の振る舞いが穏やかなモデルを選んでくれるという感じになるのでは」とのことでした。
効果量が不明なのでサンプルサイズが大き過ぎる可能性はありますが、例えば乱数シード101から115の15個でもあまり結果は変わらなかったということ、test期間の長さを1ではなく元の記事と同様の60に増やしても傾向は同様だった（むしろ{bsts}の予測性能は下がる）ということを、付記しておきます。
なお、おまけですがサンプルを乱数シード101から115の15個とし、test期間の長さを60に増やすと、このような結果になります。
> t.test(cv, ncv) # Welch's two-sample t-test Welch Two Sample t-test data: cv and ncv t = -4.8092, df = 14, p-value = 0.0002777 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: -9096169811 -3485183628 sample estimates: mean of x mean of y 3.40680e+03 6.29068e+09 > wilcox.test(cv, ncv) # Ranksum test Wilcoxon rank sum test data: cv and ncv W = 0, p-value = 1.289e-08 alternative hypothesis: true location shift is not equal to 0 > t.test(cv, bsll) Welch Two Sample t-test data: cv and bsll t = 2.1666, df = 14.002, p-value = 0.04801 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: 33.41021 6599.38457 sample estimates: mean of x mean of y 3406.79965 90.40226 > t.test(cv, bslt) Welch Two Sample t-test data: cv and bslt t = 2.1429, df = 14.005, p-value = 0.05017 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: -2.728891 6563.481958 sample estimates: mean of x mean of y 3406.7997 126.4231
test期間長が60あると、CV付き多項式フィッティングとローカルレベル状態空間モデルとでRMSEに有意差がなくなります（p = 0.04801とギリギリな上に順位和検定であれば有意差が出ますが：なお同条件でtest期間長だけ1に戻すとp = 0.0002で有意となる）。そもそもランダムウォークに対する「短期」予測という課題なのでtest期間が長くなればなるほど状態空間モデルにとっても不都合なのは直感的に分かることですが、思ったよりも予測精度の悪化が大きいなという印象を持ちました。
あとは、状態空間モデルもtrain + valでモデル推定してtrainで学習したモデルから直接val + testの長さ分の短期未来予測をtest期間の分だけ抜き出すべきかとか、逆に多項式の方をtrainで推定してvalでCVしてからベストモデルを改めてtrain + valで学習し直した方が良いのかとか、技術的に気になる点は色々ありますが今回はこの辺で。
追記
以下のようなコメントをいただきました。
ローカルレベルモデルはほぼ横置きに近いですが、試しに横置きでRMSEをシミュレートしたところ、一番結果がよかったです。
— オーレット (@quants1234) 2020年5月11日
結局、ランダムウォークなので次数を増やしたり無駄なことをすると、その分モデルリスクが大きくなるので、この点にも言及した方が良いかと思います。
確かに、結局ローカルレベルモデルの方が良いという話になるのであれば、下手に多項式モデルの次数を増やすよりは次数を低く押さえて「平ら」に近くした方が良い、というのは直感的で理解しやすいと思いました。
