5.6 変数変換を用いた予測

モデル作成時によく使われる変数変換のいくつかについては3.1節で議論しました。変換を含むモデルで予測する際は、まず変換したデータで予測を生成します。次に、元の尺度での予測を得るために変換を元に戻す(または逆変換する)必要があります。 (3.1)式で得られるBox-Cox変換の場合、逆変換は以下のようにして得られます。 \[\begin{equation} \tag{5.2} y_{t} = \begin{cases} \exp(w_{t}) & \text{$\lambda=0$の場合}\\ \text{sign}(\lambda w_t+1)|\lambda w_t+1|^{1/\lambda} & \text{それ以外} \end{cases} \end{equation}\]

fableパッケージは、モデル定義の中で変数変換が用いられている場合はいつも、自動的に予測を逆変換します。すると、逆変換された予測の分布は、「変換後は正規な」分布になります。

変換を伴う区間予測

変換を用いている場合、まず変換後の尺度で区間予測を計算し、その端の2点を逆変換して元の尺度の区間予測を得ます。このアプローチによって、区間予測の被覆確率が維持できます。ただ、点予測を中心に対称ではなくなります。

fableパッケージを使っていれば、モデル式の中で変換を使っている限り、予測の逆変換は自動的に行われます。

変数変換することで、点予測はほとんど変わらなくても、区間予測は大きく変わることが時にあります。

バイアス調整

Box-Cox変換のような算術変換を用いることの問題の一つは、逆変換後の点予測が予測の分布の平均でなくなることです。実際、逆変換後の点予測は予測の分布の中位数になるのが(変換後の空間では分布が対称と想定すると)通常です。それで構わないことも多いですが、平均の方が望ましいことも多いです。例えば、さまざまな地域ごとに売上高を予測をして、それを集計して全国予測にしたいとします。中位数は集計できませんが、平均はできます。

Box-Cox変換では、逆変換した平均は、次のようにして(おおよそ)得られます。 \[\begin{equation} \tag{5.3} \hat{y}_{T+h|T} = \begin{cases} \exp(\hat{w}_{T+h|T})\left[1 + \frac{\sigma_h^2}{2}\right] & \text{$\lambda=0$の場合}\\ (\lambda \hat{w}_{T+h|T}+1)^{1/\lambda}\left[1 + \frac{\sigma_h^2(1-\lambda)}{2(\lambda \hat{w}_{T+h|T}+1)^{2}}\right] & \text{それ以外} \end{cases} \end{equation}\] 変換後の尺度において、\(\hat{w}_{T+h|T}\)\(h\)期先予測平均、\(\sigma_h^2\)\(h\)期先予測分散です。予測の分散が大きいほど、平均と中位数との差は大きくなります。

(5.2)式によって単純に逆変換された予測と (5.3)式によって得られる平均との差をバイアスと呼びます。中位数ではなく平均を用いる場合、点予測はバイアス調整されている、と言います。

このバイアス調整でどれだけの差が生じるのか見るために、例えば、\((\lambda=0)\)として対数変換してドリフト法を使い、卵の年平均価格を予測するとしましょう。このケースでは、予測と区間予測を必ず正値にするのに、対数変換が役に立ちます。

prices %>%
  filter(!is.na(eggs)) %>%
  model(RW(log(eggs) ~ drift())) %>%
  forecast(h = 50) %>%
  autoplot(prices %>% filter(!is.na(eggs)),
    level = 80, point_forecast = lst(mean, median)
  ) +
  labs(title = "卵価格",
       y = "米ドル (インフレ調整値、セント表記) ", level = "区間予測",
       linetype = "点予測")
対数変換後データにドリフト法を適用した卵価格の予測。バイアス調整後の予測平均は実線で、予測中位数は破線で表示。

図 5.17: 対数変換後データにドリフト法を適用した卵価格の予測。バイアス調整後の予測平均は実線で、予測中位数は破線で表示。

5.17で、破線は予測の中位数を、実線は予測の平均を示しています。予測の分布が歪んでいるせいで、平均が中位数からどれだけ上方に引き上げられているか、に注目してください。この分が、バイアス調整の追加項分です。

fableパッケージでは、バイアス調整後の予測平均は自動的に計算されます。予測中位数(バイアス調整前の点予測)は、分布を表す列にmedian()関数を使うことで得られます。

prices %>%
  filter(!is.na(eggs)) %>%
  model(RW(log(eggs) ~ drift())) %>%
  forecast(h = 50) %>%
  mutate(.median = median(eggs))
#> # A fable: 50 x 5 [1Y]
#> # Key:     .model [1]
#>    .model                   year             eggs .mean .median
#>    <chr>                   <dbl>           <dist> <dbl>   <dbl>
#>  1 RW(log(eggs) ~ drift())  1994 t(N(4.1, 0.018))  61.8    61.3
#>  2 RW(log(eggs) ~ drift())  1995 t(N(4.1, 0.036))  61.4    60.3
#>  3 RW(log(eggs) ~ drift())  1996 t(N(4.1, 0.054))  61.0    59.3
#>  4 RW(log(eggs) ~ drift())  1997 t(N(4.1, 0.073))  60.5    58.4
#>  5 RW(log(eggs) ~ drift())  1998 t(N(4.1, 0.093))  60.1    57.5
#>  6 RW(log(eggs) ~ drift())  1999    t(N(4, 0.11))  59.7    56.6
#>  7 RW(log(eggs) ~ drift())  2000    t(N(4, 0.13))  59.3    55.7
#>  8 RW(log(eggs) ~ drift())  2001    t(N(4, 0.15))  58.9    54.8
#>  9 RW(log(eggs) ~ drift())  2002    t(N(4, 0.17))  58.6    53.9
#> 10 RW(log(eggs) ~ drift())  2003    t(N(4, 0.19))  58.2    53.0
#> # … with 40 more rows