ローカル線形トレンドモデル
2022-06-20
時系列分析と状態空間モデルの基礎：RとStanで学ぶ理論と実装
第5部 9章
実装：変化するトレンドのモデル化
この章で使うパッケージ ——————–
install.packages("KFAS")## Installing package into '/cloud/lib/x86_64-pc-linux-gnu-library/4.1'
## (as 'lib' is unspecified)install.packages("forecast")## Installing package into '/cloud/lib/x86_64-pc-linux-gnu-library/4.1'
## (as 'lib' is unspecified)install.packages("ggplot2")## Installing package into '/cloud/lib/x86_64-pc-linux-gnu-library/4.1'
## (as 'lib' is unspecified)install.packages("ggfortify")## Installing package into '/cloud/lib/x86_64-pc-linux-gnu-library/4.1'
## (as 'lib' is unspecified)library(KFAS)## Please cite KFAS in publications by using:
##
## Jouni Helske (2017). KFAS: Exponential Family State Space Models in R. Journal of Statistical Software, 78(10), 1-39. doi:10.18637/jss.v078.i10.library(forecast)## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoolibrary(ggplot2)
library(ggfortify)## Registered S3 methods overwritten by 'ggfortify':
## method from
## autoplot.Arima forecast
## autoplot.acf forecast
## autoplot.ar forecast
## autoplot.bats forecast
## autoplot.decomposed.ts forecast
## autoplot.ets forecast
## autoplot.forecast forecast
## autoplot.stl forecast
## autoplot.ts forecast
## fitted.ar forecast
## fortify.ts forecast
## residuals.ar forecastlibrary(gridExtra)トレンドと観測値の関係 ——————–
シミュレーションにおけるサンプルサイズ
n_sample <- 450変化しないトレンド
トレンドが0.2：毎時点0.2ずつ増えていくこと
t0 <- 0.2トレンドの累積和が、観測値となる
constant_trend <- cumsum(rep(t0, n_sample))変化するトレンド
t1 <- 0.2
t2 <- 0.4
t3 <- 0
t4 <- -0.2
trend <- c(rep(t1, 100), rep(t2, 100), rep(t3, 100), rep(t4, 150))トレンドの累積和が、観測値となる
change_trend <- cumsum(trend)図示
p1 <- autoplot(ts(constant_trend),
xlab = "Time", main = "変化しないトレンド(トレンド= 0.2)")
p2 <- autoplot(ts(change_trend),
xlab = "Time", main = "変化するトレンド")
grid.arrange(p1, p2)
シミュレーションデータの作成 ——————–
水準の過程誤差を作る
平均０、分散１の正規分布に従う、過程誤差を作る。
set.seed(12)
system_noise <- rnorm(n = n_sample)真の水準値
真の水準値＝前期の水準値＋トレンド＋過程誤差
alpha_true <- numeric(n_sample + 1)
for(i in 1:n_sample){
alpha_true[i + 1] <- alpha_true[i] + trend[i] + system_noise[i]
}cumsum関数を使っても結果は同じ
cumsum(system_noise + trend)## [1] -1.280567595 0.496601877 -0.260142602 -0.980147850 -2.777789947
## [6] -2.850085991 -2.965434703 -3.393689939 -3.300153824 -2.672139022
## [11] -3.249858604 -4.343740902 -4.923307410 -4.711355651 -4.663771889
## [16] -5.167236143 -3.778356987 -3.237844716 -2.530876544 -2.624181693
## [21] -2.200540278 0.006661179 1.218640297 1.116181052 0.290936213
## [26] 0.223551383 0.224445722 0.555568317 0.901368213 1.463432934
## [31] 2.337414098 4.609449866 4.268421217 3.397929059 3.225472326
## [36] 2.940330972 3.415115150 3.135602588 4.133707915 3.329256712
## [41] 3.634240945 2.678248055 3.456382682 2.060757026 1.952253370
## [46] 2.601719292 1.824666009 2.214663869 3.146117226 2.853518115
## [51] 3.010833203 3.098162626 3.754989874 5.975324716 5.124434654
## [56] 6.059086759 6.798336503 5.684063706 5.634024984 6.148229581
## [61] 6.754776274 7.949196874 9.004965306 9.402094223 10.436419260
## [66] 11.483209413 13.637314667 11.688054666 12.859174936 14.204236509
## [71] 13.878835883 14.329155986 14.099749374 14.117229752 14.213919286
## [76] 13.780081082 12.709027295 12.525076901 13.241832704 13.263864160
## [81] 13.468122199 12.394062648 12.391952309 13.756418189 13.933038781
## [86] 15.030195530 15.053471076 16.367180153 16.025291293 15.261892962
## [91] 15.838341361 15.053667564 16.151226967 16.480489500 17.714192481
## [96] 17.571903207 18.224184464 17.729446522 17.690432931 16.883133971
## [101] 19.209852804 19.661283168 21.861806930 21.710349169 22.217408063
## [106] 23.406532107 23.562276790 24.346456754 23.577221240 25.495493666
## [111] 25.968474843 25.840799815 25.716274235 26.863174983 26.080761012
## [116] 25.126052528 25.313767019 25.624285524 26.281318037 27.885604874
## [121] 28.111573797 29.201790588 30.015175170 28.984696123 28.085027021
## [126] 29.888376889 30.069026670 32.601517762 31.168047852 31.457538171
## [131] 32.623248524 33.346753128 33.011374217 34.557091360 34.678853922
## [136] 36.135695245 36.903209683 38.150724297 37.287903831 38.216140499
## [141] 39.705547322 39.406068958 41.227095126 43.066116804 42.170805520
## [146] 41.673745653 40.188969536 40.616467441 40.109050789 38.861489086
## [151] 39.513829125 40.556143390 40.609207282 40.493944051 41.292112974
## [156] 42.645670188 43.305849112 41.655714735 42.343685453 43.105630872
## [161] 41.818889848 43.466059650 42.095866768 42.718702634 42.558182201
## [166] 42.450052639 41.911219550 43.946868935 43.694052046 44.933265665
## [171] 44.721670374 44.749953402 44.634719444 44.692950834 45.776869614
## [176] 46.975212081 46.575407939 47.205038184 49.298586114 48.930115145
## [181] 47.786962422 46.866361516 46.590540011 46.623534585 46.509662404
## [186] 45.950408441 46.228744769 47.294681920 48.933401739 50.517342521
## [191] 52.415704013 52.526599657 54.850286876 55.768799326 56.645390283
## [196] 57.347994030 58.569226739 59.361806000 58.321355076 57.852457592
## [201] 54.866001523 54.815454241 56.114967117 54.763566040 53.813767134
## [206] 55.499868937 55.744164477 54.356084815 55.243204041 55.271376559
## [211] 55.550464044 55.474896121 55.597936568 54.566526227 54.656229475
## [216] 55.978208144 55.670534098 57.149110172 56.679574061 55.813730293
## [221] 56.813011794 55.721143931 56.574772381 55.947695964 56.368380744
## [226] 55.972507448 55.582759327 55.966332711 57.629202993 56.501063735
## [231] 56.369988861 55.953122866 56.605484169 54.957669448 54.281135631
## [236] 52.975667527 52.570287737 53.929384529 53.268135244 52.353553398
## [241] 52.481811398 52.256375401 52.926472632 53.464270257 53.409003327
## [246] 53.911105618 53.270852037 55.073717283 56.072435378 56.950163250
## [251] 56.495299633 56.406905327 56.446567992 56.993637792 56.767950609
## [256] 57.161606892 57.285188820 57.446475453 56.171659173 54.994457621
## [261] 55.121115704 56.862300072 56.276074654 55.830887573 55.743955707
## [266] 56.023786366 56.130802839 56.124787578 57.868902783 57.807533708
## [271] 55.273151516 55.666154056 56.333956065 56.156021033 55.759766787
## [276] 55.239676917 56.118360881 55.579009293 56.943233119 58.609719985
## [281] 58.373254757 58.868843920 57.990807802 57.838813455 58.976094697
## [286] 56.693048883 58.345654194 59.505462694 58.623533947 59.568822020
## [291] 59.666473869 56.728911026 56.708635931 57.098707635 55.155797181
## [296] 54.517640539 54.209974089 54.939692130 53.985217021 54.843252561
## [301] 53.161860480 52.031527599 50.814912636 50.403224349 51.856090776
## [306] 52.381387380 51.949151189 52.333659232 54.549113563 54.369811267
## [311] 54.507157388 54.466093547 55.191731140 55.831290281 54.537365681
## [316] 54.136503256 52.860153134 52.920081283 51.930420359 49.727589943
## [321] 50.119824763 49.510521670 47.851255543 47.801385517 46.835033257
## [326] 47.271993789 47.143790738 47.464466712 47.969143897 49.329769265
## [331] 50.243137830 49.474041855 48.355745221 48.155781989 47.848118501
## [336] 45.693270662 45.179498721 44.614607347 44.145484337 44.465738954
## [341] 43.850268775 42.645723707 42.612450349 42.694527328 43.738049761
## [346] 43.335871794 41.569727975 41.593130648 40.133129745 41.502694496
## [351] 40.997860163 39.306746428 39.339429880 39.705132075 39.178730325
## [356] 39.902167521 38.561909826 37.623731324 36.368267040 36.666488661
## [361] 35.602642698 34.617203417 34.912164703 34.556889932 34.982977791
## [366] 33.844689096 33.613803479 33.632822002 33.567362337 34.213101923
## [371] 33.804192875 33.637593700 34.249916998 36.451003384 36.896916545
## [376] 37.238506616 37.741012695 36.402735094 36.208931742 34.748718463
## [381] 33.297603803 32.658737426 32.628520292 31.802765571 33.233618445
## [386] 31.487575102 31.101752290 31.927547956 31.632422425 31.419956421
## [391] 31.190162488 33.584516534 34.720857328 33.006003614 33.998448455
## [396] 34.910585779 34.643955267 33.638630318 32.671826523 32.222766227
## [401] 32.109891994 30.541707634 28.951624180 29.056891937 29.012263320
## [406] 29.026715384 28.646481695 27.843872194 29.383894186 30.268780088
## [411] 31.377227096 30.217566463 30.498304644 28.732360029 29.731139731
## [416] 30.532240796 30.360255288 29.868522474 30.735494217 30.271268147
## [421] 30.601081767 29.459421128 27.747244404 27.791087846 27.483454683
## [426] 26.978774745 25.827874325 25.850102223 23.715913946 22.503929023
## [431] 22.933293209 22.919459804 23.023835502 23.161152201 23.308444695
## [436] 21.891186299 21.008606017 19.010944900 18.740493137 18.391138037
## [441] 17.724117571 18.177944202 17.365506570 18.327527091 18.226789476
## [446] 18.753619622 19.267089051 19.539036360 19.089967029 19.392886422sum(alpha_true[-1] - cumsum(system_noise + trend))## [1] -8.445566e-11観測誤差を作る
平均０、標準偏差５の正規分布に従う観測誤差を組み込む。
obs_noise <- rnorm(n = n_sample, sd=5)トレンドが変化する売り上げデータ
11を足しているのは負にならないようにするため。
sales <- alpha_true[-1] + obs_noise + 11架空の売り上げデータを図示
autoplot(ts(sales), main = "架空の売り上げデータ")
200時点目までは増えていき、300まで上げどまり、その後減少トレンドに転じるデータを作成することができた。
最後の50期間は、テストのために残しておく
sales_train <- sales[1:400]
sales_test <- sales[401:450]途中の50期間は欠損とする
sales_train[125:175] <- NAKFASによるローカル線形トレンドモデル ——————–
Step1：モデルの構造を決める
SSMtrendのdegree=2:ローカル線形トレンドモデル
build_trend <- SSModel(
H = NA,
sales_train ~ SSMtrend(degree = 2, Q = c(list(NA), list(NA)))
)Step2：パラメタ推定
推定するパラメタは3つ：観測誤差のH、過程誤差の分散Q（水準の変動、トレンドの変動）
fit_trend <- fitSSM(build_trend, inits = c(1, 1, 1))Step3、4：フィルタリング・スムージング
result_trend <- KFS(
fit_trend$model,
filtering = c("state", "mean"),
smoothing = c("state", "mean")
)推定されたパラメタ
観測誤差の分散
fit_trend$model$H## , , 1
##
## [,1]
## [1,] 23.63745観測誤差の分散:23.63
過程誤差の分散
fit_trend$model$Q## , , 1
##
## [,1] [,2]
## [1,] 0.8833319 0.0000000000
## [2,] 0.0000000 0.0004648067水準の変動：0.88
トレンドの変動：0.00046
補足：モデルの行列表現 ——————–
fit_trend$model$T## , , 1
##
## level slope
## level 1 1
## slope 0 1fit_trend$model$R## , , 1
##
## [,1] [,2]
## level 1 0
## slope 0 1fit_trend$model$Z## , , 1
##
## level slope
## [1,] 1 0トレンドの図示 ——————–
levelは推定された売り上げの状態。slopeはトレンド。
推定された平滑化状態はローカルレベルモデルと異なり、2列ある。
head(result_trend$alphahat, n = 3)## Time Series:
## Start = 1
## End = 3
## Frequency = 1
## level slope
## 1 7.592715 0.1979556
## 2 7.573414 0.1980699
## 3 7.817400 0.1981601データの整形
trend_df <- data.frame(
time = 1:length(sales_train),
true_trend = trend[1:length(sales_train)],
estimate_trend = result_trend$alphahat[, "slope"]
)図示
ggplot(data = trend_df, aes(x = time, y = true_trend)) +
labs(title="トレンドの変化") +
geom_line(aes(y = true_trend), size = 1.2, linetype="dashed") +
geom_line(aes(y = estimate_trend), size = 1.2)
点線：真のトレンド
実線：推定されたトレンド
おおよその傾向を捕らえることができているのがわかる。
補間と予測 ——————–
平滑化状態と予測区間interval_trend <- predict(
interval_trend <- predict(
fit_trend$model, interval = "prediction", level = 0.95)将来予測の結果と予測区間
forecast_trend <- predict(
fit_trend$model, interval = "prediction", level = 0.95, n.ahead = 50)過去の状態と予測結果をまとめた
estimate_all <- rbind(interval_trend, forecast_trend)ローカル線形トレンドモデルによる予測の考え方 ——————–
最後のレベルと傾きを取得
last_level <- tail(result_trend$a[, "level"], n = 1)
last_trend <- tail(result_trend$a[, "slope"], n = 1)
last_level # 水準## Time Series:
## Start = 401
## End = 401
## Frequency = 1
## [1] 44.60478last_trend # トレンド## Time Series:
## Start = 401
## End = 401
## Frequency = 1
## [1] -0.1395292トレンド成分を足していくだけの予測
fore <- cumsum(c(last_level, rep(last_trend, 49)))参考 for構文を使った場合
fore2 <- numeric(50)
fore2[1] <- last_level
for(i in 2:50){
fore2[i] <- fore2[i - 1] + last_trend
}
sum(fore2 - fore)## [1] -3.225864e-12パッケージの予測結果と比較
forecast_trend[, "fit"]## Time Series:
## Start = 401
## End = 450
## Frequency = 1
## [1] 44.60478 44.46525 44.32572 44.18619 44.04666 43.90713 43.76760 43.62807
## [9] 43.48854 43.34902 43.20949 43.06996 42.93043 42.79090 42.65137 42.51184
## [17] 42.37231 42.23278 42.09325 41.95372 41.81419 41.67467 41.53514 41.39561
## [25] 41.25608 41.11655 40.97702 40.83749 40.69796 40.55843 40.41890 40.27937
## [33] 40.13984 40.00032 39.86079 39.72126 39.58173 39.44220 39.30267 39.16314
## [41] 39.02361 38.88408 38.74455 38.60502 38.46549 38.32596 38.18644 38.04691
## [49] 37.90738 37.76785fore## [1] 44.60478 44.46525 44.32572 44.18619 44.04666 43.90713 43.76760 43.62807
## [9] 43.48854 43.34902 43.20949 43.06996 42.93043 42.79090 42.65137 42.51184
## [17] 42.37231 42.23278 42.09325 41.95372 41.81419 41.67467 41.53514 41.39561
## [25] 41.25608 41.11655 40.97702 40.83749 40.69796 40.55843 40.41890 40.27937
## [33] 40.13984 40.00032 39.86079 39.72126 39.58173 39.44220 39.30267 39.16314
## [41] 39.02361 38.88408 38.74455 38.60502 38.46549 38.32596 38.18644 38.04691
## [49] 37.90738 37.76785predict関数の結果とforループの結果が一致していることを確認できた。
ローカル線形トレンドモデルは「最新のデータを用いて補正されたトレンド」を用いて将来を予測していることを確認できた。
補間と予測結果の図示 ——————–
データの整形
df <- cbind(
data.frame(sales = sales, time = 1:n_sample),
as.data.frame(estimate_all)
)図示
ggplot(data = df, aes(x = time, y = sales)) +
labs(title="トレンドが変わる売り上げの予測") +
geom_point(alpha = 0.6, size = 0.9) +
geom_line(aes(y = fit), size = 1.2) +
geom_ribbon(aes(ymin = lwr, ymax = upr), alpha = 0.3)
400時点以降が予測結果。減少トレンドを上手く表現できていることがわかった。
ARIMAによる予測 ——————–
モデルの構築
mod_arima <- auto.arima(sales_train)予測
forecast_arima <- forecast(mod_arima, h=50, level = 0.95)#予測結果の図示
autoplot(forecast_arima, main = "ARIMAによる予測",
predict.colour = 1, shadecols = "gray")
右肩上がりのトレンドがあるとみなして予測を出してしまっていることがわかった。300時点以降で見ると減少トレンドであるにもかかわらず、0時点から全データで見ているため、上昇トレンドであると予測してしまっている。
予測精度の比較
accuracy(forecast_trend[, "fit"], sales_test)["Test set", "RMSE"]## [1] 7.47772accuracy(forecast_arima, sales_test)["Test set", "RMSE"]## [1] 13.08098予測精度も大幅に劣ることがわかった。arimaのほうが予測誤差が大きいことがわかった。