5.5 予測の分布と区間予測

予測の分布

1.7節で議論したように、予測の不確実性は確率分布を使って表現します。この分布を使うことで、当てはめたモデルを使って起こり得る将来の観測値の確率を描写します。点予測はこの分布の平均です。ほとんどの時系列モデルは正規分布に従う予測を生成します。つまり、あり得る将来値の分布は正規分布に従うと想定しています。本章の後半で、正規分布に替わるものを2、3見ていきます。

区間予測

区間予測は、指定した確率で\(y_{t}\)がその範囲内に含まれるはずの区間です。例えば、将来の観測値が正規分布に従うと仮定して、\(h\)期先予測の95%区間予測は、 \[ \hat{y}_{T+h|T} \pm 1.96 \hat\sigma_h \] となります。ここで\(\hat\sigma_h\)\(h\)期先予測の分布の標準偏差の推計値です。

より一般的に、区間予測は以下のように書けます。 \[ \hat{y}_{T+h|T} \pm c \hat\sigma_h \] ここで乗数\(c\)は被覆確率に依存します。本書では通常80%区間と95%区間を計算しますが、どんなパーセントでも使用可能です。表5.1は、正規分布を想定した被覆確率ごとの乗数\(c\)の値を与えてくれます。

表 5.1: 区間予測に用いられる乗数
被覆確率(%) 乗数
50 0.67
55 0.76
60 0.84
65 0.93
70 1.04
75 1.15
80 1.28
85 1.44
90 1.64
95 1.96
96 2.05
97 2.17
98 2.33
99 2.58

区間予測の被覆確率は予測の不確実性を表現しています。点予測だけでは、予測がどれだけ正確か述べられません。区間予測も合わせて生成すれば、各予測がどれだけの不確実性を伴っているか明らかにできます。この理由で、区間予測を伴わない点予測は、ほぼ無価値かもしれません。

1期先の区間予測

1期先を予測する際、予測の分布の標準偏差は残差の標準偏差を使って、以下のように推計できます。 \[\begin{equation} \hat{\sigma} = \sqrt{\frac{1}{T-K}\sum_{t=1}^T e_t^2} \tag{5.1} \end{equation}\] ここで\(K\)は、予測手法の中で推計されたパラメータの数です。

例えば、Google株価データgoogle_2015(図5.8で示しました)のナイーブ予測を考えてみましょう。観測系列の最後の値は758.88ですから、次期の株価の予測は758.88になります。ナイーブ法の残差の標準偏差は、(5.1)式から11.19です。ですから、Google株価の次期の値の95%区間予測は、 \[ 758.88 \pm 1.96(11.19) = [736.9, 780.8] \] 同様に、80%区間予測は、 \[ 758.88 \pm 1.28(11.19) = [744.5, 773.2] \]

になります。乗数の値(1.96と1.28)は表5.1から取っています。

複数期先の区間予測

区間予測に共通する特徴は、予測期間が先に延びるほど通常幅が広くなることです。先を予測すればするほど、予測に伴う不確実性が増え、区間予測の幅が広がるというわけです。つまり、\(\sigma_h\)は通常\(h\)に連れて増加します。(ただ、この性質を持たない非線形モデルもいくつかあります。)

区間予測の生成には、\(\sigma_h\)を推計する必要があります。既に述べたように、1期先予測(\(h=1\))では、(5.1)式が予測の標準偏差\(\sigma_1\)の良い推計値を提供してくれます。複数期先の予測については、より複雑な計算手法が必要になります。その計算は無相関な残差を想定しています。

ベンチマークとなる手法

4つのベンチマークとなる手法では、無相関な残差を想定すれば、予測の標準偏差を算術的に導き出せます。\(\hat{\sigma}_h\)\(h\)期先の予測の分布の標準偏差を表し、\(\hat{\sigma}\)(5.1)式から得られる残差の標準偏差とすると、表5.2に示した式が使えます。これらの式全てで、\(h=1\)\(T\)が大きい場合、大体同じ値の\(\hat{\sigma}\)になることに注意してください。

表 5.2: 4つのベンチマークとなる手法での複数スッテプ先予測の標準偏差。ただし、\(\sigma\)は残差の標準偏差、\(m\)は季節周期、\(k\)\((h-1) /m\)の整数部分(つまり、\(T+h\)時点以前の予測期間に含まれる完全な年の数)。
ベンチマークとなる手法 \(h\)期先予測の標準偏差
平均 \(\hat\sigma_h = \hat\sigma\sqrt{1 + 1/T}\)
ナイーブ \(\hat\sigma_h = \hat\sigma\sqrt{h}\)
季節ナイーブ \(\hat\sigma_h = \hat\sigma\sqrt{k+1}\)
ドリフト \(\hat\sigma_h = \hat\sigma\sqrt{h(1+h/T)}\)

区間予測は、fableパッケージを使っていれば簡単に計算できます。例えば、Google株価にナイーブ法を用いた場合、以下のように出力できます。

google_2015 %>%
  model(NAIVE(Close)) %>%
  forecast(h = 10) %>%
  hilo()
#> # A tsibble: 10 x 7 [1]
#> # Key:       Symbol, .model [1]
#>    Symbol .model      day        Close .mean            `80%`            `95%`
#>    <chr>  <chr>     <dbl>       <dist> <dbl>           <hilo>           <hilo>
#>  1 GOOG   NAIVE(Cl…   253  N(759, 125)  759. [744.5, 773.2]80 [736.9, 780.8]95
#>  2 GOOG   NAIVE(Cl…   254  N(759, 250)  759. [738.6, 779.2]80 [727.9, 789.9]95
#>  3 GOOG   NAIVE(Cl…   255  N(759, 376)  759. [734.0, 783.7]80 [720.9, 796.9]95
#>  4 GOOG   NAIVE(Cl…   256  N(759, 501)  759. [730.2, 787.6]80 [715.0, 802.7]95
#>  5 GOOG   NAIVE(Cl…   257  N(759, 626)  759. [726.8, 790.9]80 [709.8, 807.9]95
#>  6 GOOG   NAIVE(Cl…   258  N(759, 751)  759. [723.8, 794.0]80 [705.2, 812.6]95
#>  7 GOOG   NAIVE(Cl…   259  N(759, 876)  759. [720.9, 796.8]80 [700.9, 816.9]95
#>  8 GOOG   NAIVE(Cl…   260 N(759, 1002)  759. [718.3, 799.4]80 [696.8, 820.9]95
#>  9 GOOG   NAIVE(Cl…   261 N(759, 1127)  759. [715.9, 801.9]80 [693.1, 824.7]95
#> 10 GOOG   NAIVE(Cl…   262 N(759, 1252)  759. [713.5, 804.2]80 [689.5, 828.2]95

hilo()関数は、予測の分布を区間予測に変換します。デフォルトでは、80%と95%の区間予測を出力しますが、level引数で他のパーセントにもできます。

プロットする際は、区間予測は網掛け部分で示し、色の濃さで区間の確率を表します。ここでも、デフォルトで80%と95%の区間予測を示していますが、level引数で他のパーセントにもできます。

google_2015 %>%
  model(NAIVE(Close)) %>%
  forecast(h = 10) %>%
  autoplot(google_2015) +
  labs(title="Googleの日次株価終値", y="米ドル", level = "区間予測")
ナイーブ法に基づくGoogle株価終値の80%と95%の区間予測

図 5.14: ナイーブ法に基づくGoogle株価終値の80%と95%の区間予測

ブートストラップ法で生成した残差による区間予測

残差に正規性を仮定するのが理にかなっていない場合でも、無相関と均一分散さえ仮定できるなら、ブートストラップ法を使うことが一つの代替手段になります。

1期先予測誤差は\(e_t = y_t - \hat{y}_{t|t-1}\)と定義されます。それは、以下のように書き直せます。 \[ y_t = \hat{y}_{t|t-1} + e_t \] ですから、時系列の次の観測値を以下を使ってシミュレーションできます。 \[ y_{T+1} = \hat{y}_{T+1|T} + e_{T+1} \] ここで\(\hat{y}_{T+1|T}\)は1期先の予測で、\(e_{T+1}\)は未知の将来誤差です。将来誤差は過去の誤差と似ていると想定すると、\(e_{T+1}\)を過去に観測した誤差(つまりは、残差)の集まりから標本抽出したもので代替できます。シミュレーションで得た新たな観測値をデータセットに加えるプロセスを繰り返すことで、以下が得られます。 \[ y_{T+2} = \hat{y}_{T+2|T+1} + e_{T+2} \] ここで\(e_{T+2}\)は残差の集合から抽出した値です。このように続けると、時系列の将来値の経路をシミュレーションできます。

これを繰り返し行うと、起こり得る将来値を複数得ることになります。そのいくつかを見るには、generate()関数が使えます。

fit <- google_2015 %>%
  model(NAIVE(Close))
sim <- fit %>% generate(h = 30, times = 5, bootstrap = TRUE)
sim
#> # A tsibble: 150 x 5 [1]
#> # Key:       Symbol, .model, .rep [5]
#>    Symbol .model         day .rep   .sim
#>    <chr>  <chr>        <dbl> <chr> <dbl>
#>  1 GOOG   NAIVE(Close)   253 1      756.
#>  2 GOOG   NAIVE(Close)   254 1      749.
#>  3 GOOG   NAIVE(Close)   255 1      751.
#>  4 GOOG   NAIVE(Close)   256 1      750.
#>  5 GOOG   NAIVE(Close)   257 1      754.
#>  6 GOOG   NAIVE(Close)   258 1      754.
#>  7 GOOG   NAIVE(Close)   259 1      758.
#>  8 GOOG   NAIVE(Close)   260 1      763.
#>  9 GOOG   NAIVE(Close)   261 1      759.
#> 10 GOOG   NAIVE(Close)   262 1      748.
#> # … with 140 more rows

ここでは、次の30営業日の起こり得る標本経路を5つ生成しています。.rep変数はtsibbleの新しいキー変数です。以下のプロットは過去のデータと5つの標本経路を示しています。

google_2015 %>%
  ggplot(aes(x = day)) +
  geom_line(aes(y = Close)) +
  geom_line(aes(y = .sim, colour = as.factor(.rep)),
    data = sim) +
  labs(title="Googleの日次株価終値", y="米ドル" ) +
  guides(colour = "none")
ナイーブ法に基づくGoogle株価終値の、残差からのブートストラップによる5つのシミュレートされた将来の標本経路

図 5.15: ナイーブ法に基づくGoogle株価終値の、残差からのブートストラップによる5つのシミュレートされた将来の標本経路

次に、各予測期間ごとに将来の標本経路の百分位点を算出することで区間予測を計算できます。結果は、ブートストラップによる区間予測と呼びます。「ブートストラップ」という名前は、自助努力(ブートストラップ)で自分自身を引っ張り出すことを指しています。この手順によって、過去のデータだけを使って将来の不確実性を測れるのです。

これは全てforecast()関数に組み込まれているので、直接generate()を呼び出す必要はありません。

fc <- fit %>% forecast(h = 30, bootstrap = TRUE)
fc
#> # A fable: 30 x 5 [1]
#> # Key:     Symbol, .model [1]
#>    Symbol .model         day        Close .mean
#>    <chr>  <chr>        <dbl>       <dist> <dbl>
#>  1 GOOG   NAIVE(Close)   253 sample[5000]  759.
#>  2 GOOG   NAIVE(Close)   254 sample[5000]  759.
#>  3 GOOG   NAIVE(Close)   255 sample[5000]  758.
#>  4 GOOG   NAIVE(Close)   256 sample[5000]  759.
#>  5 GOOG   NAIVE(Close)   257 sample[5000]  759.
#>  6 GOOG   NAIVE(Close)   258 sample[5000]  759.
#>  7 GOOG   NAIVE(Close)   259 sample[5000]  759.
#>  8 GOOG   NAIVE(Close)   260 sample[5000]  759.
#>  9 GOOG   NAIVE(Close)   261 sample[5000]  759.
#> 10 GOOG   NAIVE(Close)   262 sample[5000]  759.
#> # … with 20 more rows

予測の分布が今や5000の標本経路から成るシミュレーションで表されていることに気付いてください。正規分布の想定がないので、区間予測は対称ではありません。.mean列はブートストラップ標本の平均なので、ブートストラップ法ではない方法で得られた結果とは微妙に異なっているかもしれません。

autoplot(fc, google_2015) +
  labs(title="Googleの日次株価終値", y="米ドル", level = "区間予測")
ナイーブ法に基づくGoogle株価終値予測、残差からのブートストラップ法を使用

図 5.16: ナイーブ法に基づくGoogle株価終値予測、残差からのブートストラップ法を使用

forecast()times引数で、標本数は調整できます。例えば、1000個のブートストラップ標本に基づく区間は、以下のようにします。

google_2015 %>%
  model(NAIVE(Close)) %>%
  forecast(h = 10, bootstrap = TRUE, times = 1000) %>%
  hilo()
#> # A tsibble: 10 x 7 [1]
#> # Key:       Symbol, .model [1]
#>    Symbol .model      day        Close .mean            `80%`            `95%`
#>    <chr>  <chr>     <dbl>       <dist> <dbl>           <hilo>           <hilo>
#>  1 GOOG   NAIVE(Cl…   253 sample[1000]  760. [748.2, 770.8]80 [743.9, 777.6]95
#>  2 GOOG   NAIVE(Cl…   254 sample[1000]  760. [743.9, 776.1]80 [734.1, 801.6]95
#>  3 GOOG   NAIVE(Cl…   255 sample[1000]  760. [739.5, 781.7]80 [728.6, 809.0]95
#>  4 GOOG   NAIVE(Cl…   256 sample[1000]  760. [736.7, 784.7]80 [723.4, 813.1]95
#>  5 GOOG   NAIVE(Cl…   257 sample[1000]  760. [734.4, 787.2]80 [719.4, 819.7]95
#>  6 GOOG   NAIVE(Cl…   258 sample[1000]  760. [731.5, 790.2]80 [717.8, 820.3]95
#>  7 GOOG   NAIVE(Cl…   259 sample[1000]  761. [730.4, 793.0]80 [713.0, 826.3]95
#>  8 GOOG   NAIVE(Cl…   260 sample[1000]  761. [726.2, 796.2]80 [706.3, 830.7]95
#>  9 GOOG   NAIVE(Cl…   261 sample[1000]  761. [723.5, 800.2]80 [707.5, 841.0]95
#> 10 GOOG   NAIVE(Cl…   262 sample[1000]  760. [719.2, 801.8]80 [701.9, 841.4]95

このケースでは、正規分布に基づく区間予測と似ています(が、同じではありません)。

下のスライダーを使いブートストラップ標本数(times)を変え、予測の分布にどう影響するか確認してください。