8.6 モデルの推計と選択

ETSモデルの推計

パラメータを推計するのに、2乗誤差和最小化の替わりになるのは、「尤度」最大化です。尤度とは、指定されたモデルからデータ(実績値)が生じる確率です。こうして、大きな尤度は良いモデルと結び付いています。加法誤差モデルでは、(誤差の正規分布を想定した)尤度最大化は2乗誤差和最小化と同じ結果になります。しかし、乗法誤差モデルでは、違った結果になります。この節では、尤度最大化を使って、平滑化パラメータ\(\alpha\)\(\beta\)\(\gamma\)\(\phi\)と初期状態\(\ell_0\)\(b_0\)\(s_0,s_{-1},\dots,s_{-m+1}\)を推計します。

平滑化パラメータが取り得る値は制約されています。伝統的に、方程式を加重平均と解釈できるよう、パラメータは0と1の間に制約されています。つまり、\(0< \alpha,\beta^*,\gamma^*,\phi<1\)です。状態空間モデルでは、\(\beta=\alpha\beta^*\)\(\gamma=(1-\alpha)\gamma^*\)と設定しました。ですから、伝統的制約は\(0< \alpha <1\)\(0 < \beta < \alpha\)\(0< \gamma < 1-\alpha\)に翻訳されます。実務上、減衰パラメータ\(\phi\)は、モデル推計での数値的困難を避けるため、さらに制約されるのが普通です。fableパッケージでは、\(0.8<\phi<0.98\)に制約しています。

パラメータを見るもう一つのやり方は、状態空間モデルの数学的性質の考慮を通じてです。パラメータを制約する理由は、遠い過去の観測値が現在の予測に影響を保持するのを防ぐためです。これは、パラメータへのいくつかの許容性制約につながります。この制約は、伝統的制約の範囲よりも(いつもではありませんが)通常緩くなっています。 (Hyndman et al., 2008, p. Ch10) 例えば、ETS(A,N,N)モデルでは、パラメータの伝統的制約の範囲は\(0< \alpha <1\)ですが、許容性制約の範囲は\(0< \alpha <2\)です。ETS(A,A,N)モデルでは、パラメータの伝統的制約の範囲は\(0<\alpha<1\)\(0<\beta<\alpha\)ですが、許容性制約の範囲は\(0<\alpha<2\)\(0<\beta<4-2\alpha\)です。

モデル選択

ETSの統計的フレームワークの大きな利点は、モデル選択に情報規準が使えることです。7.5節で紹介したAIC, AIC\(_{\text{c}}\)、BICを、ここでは所与の時系列に最適なETSモデルはどれかを決めるのに使えます。

ETSモデルでは、赤池情報量規準(AIC)は以下のように定義されます。 \[ \text{AIC} = -2\log(L) + 2k \] ただし、\(L\)はモデルの尤度、\(k\)は推計されるパラメータと初期状態の総数(残差の分散を含む)です。

標本規模が小さい際のバイアスを修正した赤池情報量規準(AIC\(_\text{c}\))は以下のように定義されます。 \[ \text{AIC}_{\text{c}} = \text{AIC} + \frac{2k(k+1)}{T-k-1} \] そして、ベイズ情報量規準(BIC)は以下のように定義されます。 \[ \text{BIC} = \text{AIC} + k[\log(T)-2] \]

(誤差、トレンド、季節性)の組み合わせのうち、3つは数値的困難につながりかねません。具体的には、そうした不安定化を引き起こしかねないモデルは、ETS(A,N,M)、ETS(A,A,M)、ETS(A,A\(_d\),M)です。状態方程式に、ゼロに近くなり得る値での割り算があるためです。モデル選択の際に、これらの特殊な組み合わせは考えることは通常ありません。

乗法誤差を持つモデルは、データが必ず正の場合に役立ちますが、データにゼロや負値がある場合は数値的に安定しません。ですから、時系列が正とは限らない場合は、乗法誤差モデルは考えません。そのケースでは、完全に加法なモデル6つだけを適用します。

事例: オーストラリア国内、休暇旅行

では、オーストラリア国内、休暇旅行の宿泊客数の2016–2019年を予測するため、ETSの統計的フレームワークを使いましょう。ETS()関数にAICcを最小にするモデルを選択させます。

aus_holidays <- tourism %>%
  filter(Purpose == "Holiday") %>%
  summarise(Trips = sum(Trips)/1e3)
fit <- aus_holidays %>%
  model(ETS(Trips))
report(fit)
#> Series: Trips 
#> Model: ETS(M,N,A) 
#>   Smoothing parameters:
#>     alpha = 0.3484 
#>     gamma = 1e-04 
#> 
#>   Initial states:
#>   l[0]    s[0]   s[-1]   s[-2] s[-3]
#>  9.727 -0.5376 -0.6884 -0.2934 1.519
#> 
#>   sigma^2:  0.0022
#> 
#>   AIC  AICc   BIC 
#> 226.2 227.8 242.9

選択されたモデルはETS(M,N,A)でした。 \[\begin{align*} y_{t} &= (\ell_{t-1}+s_{t-m})(1 + \varepsilon_t)\\ \ell_t &= \ell_{t-1} + \alpha(\ell_{t-1}+s_{t-m})\varepsilon_t\\ s_t &= s_{t-m} + \gamma(\ell_{t-1}+s_{t-m}) \varepsilon_t \end{align*}\]

パラメータの推計値は、\(\hat\alpha= 0.3484\)\(\hat\gamma=0.0001\)です。出力には、初期状態の推計値\(\ell_0\)\(s_{0}\)\(s_{-1}\)\(s_{-2}\)\(s_{-3}\)も出ています。これらを、表8.3で示した、加法季節性を持つHolt-Winters法で得た値と、比べてみてください。

8.10は時を経るにつれた状態の推移を、図8.12はモデルから生成された点予測と区間予測を、示しています。\(\gamma\)の値が小さいことは、季節成分が時を経てもほとんど変わらないことを示しています。

components(fit) %>%
  autoplot() +
  labs(title = "ETS(M,N,A) components")
推計された状態の推移を表すグラフ

図 8.10: 推計された状態の推移を表すグラフ

このモデルは乗法誤差を持っているため、イノベーション誤差は通常の誤差(つまり、訓練データの1期先誤差)と同じではありません。イノベーション誤差は\(\hat{\varepsilon}_t\)で与えられ、通常の誤差は\(y_t - \hat{y}_{t|t-1}\)と定義されます。augment()関数を使えば、両方とも得ることができます。それらを図8.11にプロットしました。

ETS(M,N,A)モデルからの残差と1期先予測の誤差

図 8.11: ETS(M,N,A)モデルからの残差と1期先予測の誤差

参考文献

Hyndman, R. J., Koehler, A. B., Ord, J. K., & Snyder, R. D. (2008). Forecasting with exponential smoothing: The state space approach. Springer-Verlag. http://www.exponentialsmoothing.net