5.4 残差診断

良い予測手法から生じるイノベーション残差には次のような性質があります。

  1. イノベーション残差が無相関。もしイノベーション残差の間に相関があるなら、予測に使えたはずの情報が、残差に漏れたことになる。
  2. イノベーション残差の平均がゼロ。もし平均がゼロでないなら、予測にバイアスがあることになる。

これらの性質を満たさない予測手法には、改善の余地があります。だからと言って、これらの性質を満たす予測手法には、改善の余地がないと言っているわけではありません。同じデータセットについて、いくつかの異なる予測手法が全てこれらの性質を満たしていることは、あり得ることです。これらの性質をチェックすることは、予測手法が利用可能な情報を全て使っているかを見るために大切なことです。しかし、それで予測手法を選択するのは、良いやり方ではありません。

これらの性質のどちらか一方でも満たされない場合、予測手法を修正することでより良い予測にできます。バイアスの調整は簡単です。残差の平均が\(m\)なら、全ての予測から\(m\)を差し引くだけで、バイアス問題は解決です。相関問題の解決はそれよりも難しいので、10章まで先延ばしにしておきます。

これらの必須の性質に加えて、残差が以下のような性質を有していることは、(必要ではありませんが)役に立ちます。

  1. イノベーション残差の分散が均一。これは「均一分散」として知られる。
  2. イノベーション残差が正規分布に従う。

これら2つの性質を満たしていると、区間予測の計算が容易になります。(例えば、5.5節を参照) しかし、これら2つの性質を満たしていない予測手法に、必ず改善余地があるわけではありません。Box-Cox変換を適用することが、これら2つの性質を満たす助けになることがありますが、そうでなければ、イノベーション残差が均一分散で正規分布であることを確保するために、できることは通常ほとんどありません。その場合、区間予測を得るための代替手法が必要になるだけです。残差が非正規分布に従う場合の扱い方は5.5節で示します。

事例: Googleの日次株価終値の予測

5.2節に引き続き、Googleの日次株価終値を例とします。個別銘柄の株価や株価指数の予測手法として、しばしばナイーブ法が最良です。つまり、各予測を単に最後の観測値と等しくする、\(\hat{y}_{t} = y_{t-1}\)とするということです。すると、残差は単に隣り合う観測値間の差分に等しくなります。 \[ e_{t} = y_{t} - \hat{y}_{t} = y_{t} - y_{t-1} \]

以下のグラフは、2015年の営業日におけるGoogleの日次株価終値です。大きなジャンプは2015年7月17日で、この日は第2四半期業績が期待以上に良かったため株価が16%上昇しました。(google_2015オブジェクトは5.2節で作成しました)

autoplot(google_2015, Close) +
  labs(y = "米ドル",
       title = "2015年のGoogleの日次株価終値")
2015年のGoogleの日次株価終値

図 5.9: 2015年のGoogleの日次株価終値

ナイーブ法を用いてこの系列の予測から得た残差を図5.10に示します。一番大きな正の残差は7月の予期せぬ株価上昇の結果です。

aug <- google_2015 %>%
  model(NAIVE(Close)) %>%
  augment()
autoplot(aug, .innov) +
  labs(y = "米ドル",
       title = "ナイーブ法からの残差")
ナイーブ法を用いたGoogle株価予測からの残差

図 5.10: ナイーブ法を用いたGoogle株価予測からの残差

aug %>%
  ggplot(aes(x = .innov)) +
  geom_histogram() +
  labs(title = "残差のヒストグラム")
Google株価にナイーブ法を適用して得た残差のヒストグラム。正規分布にしては、右裾がやや長過ぎる感じ。

図 5.11: Google株価にナイーブ法を適用して得た残差のヒストグラム。正規分布にしては、右裾がやや長過ぎる感じ。

aug %>%
  ACF(.innov) %>%
  autoplot() +
  labs(title = "ナイーブ法からの残差")
Google株価にナイーブ法を適用して得た残差のACF。相関がないことは、良い予測であることを示唆している。

図 5.12: Google株価にナイーブ法を適用して得た残差のACF。相関がないことは、良い予測であることを示唆している。

これらのグラフは、ナイーブ法が全ての利用可能な情報を取り込んだ予測を生成したらしいことを示しています。残差の平均はゼロに近く、残差系列に有意な相関はありません。残差の時間プロットは、残差の変動は、1つの外れ値を除き、このデータの期間中は同程度であることを示しており、従って、残差の分散は均一として扱えます。このことは、残差のヒストグラムからも読み取れます。ヒストグラムからは外れ値を無視したとしても、右裾がやや長過ぎるので、残差が正規分布に従わないことが示唆されます。それ故、この手法からの予測はたぶん相当良いけれど、正規分布を想定して計算する区間予測は不正確かもしれません。

これら残差診断グラフを生成する便利なショートカットは、 gg_tsresiduals()関数を使うことです。残差の時間プロット、ACFプロットとヒストグラムを生成します。

google_2015 %>%
  model(NAIVE(Close)) %>%
  gg_tsresiduals()
Google株価にナイーブ法を適用した場合の残差診断グラフ

図 5.13: Google株価にナイーブ法を適用した場合の残差診断グラフ

自己相関のかばん(portmanteau)検定

ACFプロットを見ることに加え、\(r_k\)値を一つ一つ別々に扱うのではなく、グループとして全体の\(r_k\)値の集合と考えることで、より正式な自己相関検定を行うことができます。

\(r_k\)はラグ次数が\(k\)の自己相関係数であることを思い起こしましょう。ACFプロットからスパイクの一つ一つが青色の破線の範囲内にあるか見ているとき、暗黙に複数の仮説検定を行っていることになります。そして、一つ一つの検定には小さいながら偽陽性の可能性があり得ます。十分な数の検定を行えば、少なくとも1つは偽陽性を示す可能性があり、実際には残差に自己相関は残っていないのに、残っていると結論することになりかねません。

この問題を克服するためには、最初の\(\ell\)個の自己相関係数がホワイトノイズ過程からのものと有意に異なるか、検定します。複数の自己相関係数をまとめた検定はかばん検定と呼ばれ、その名は複数の衣服を運ぶスーツケースやコート掛けを意味するフランス語に由来します。

そうした検定の一つが、以下の統計量に基づくBox-Pierce検定です。 \[ Q = T \sum_{k=1}^\ell r_k^2 \] ここで\(\ell\)は考慮する最大ラグ次数、\(T\)は観測値の数です。もし一つ一つの\(r_k\)が全てゼロに近ければ、\(Q\)値は小さくなります。もしいくつかの\(r_k\)が(正であれ負であれ)大きいなら、\(Q\)値は大きくなります。季節性のないデータには\(\ell=10\)を、季節性のあるデータには\(m\)を季節周期として\(\ell=2m\)を提案します。しかし、\(\ell\)が大きい場合は上手く行きません。ですから、提案した\(\ell\)\(T/5\)より大きくなる場合は、\(\ell=T/5\)としましょう。

関連する(そして、より正確な)検定が、Ljung-Box検定で、以下に基づいています。 \[ Q^* = T(T+2) \sum_{k=1}^\ell (T-k)^{-1}r_k^2 \]

ここでも、大きな値の\(Q^*\)値は、自己相関がホワイトノイズ系列から計算したものと違うことを示唆します。

どれだけ大きいと大き過ぎるのでしょう? もし自己相関がホワイトノイズ系列から来ていたとすると、\(Q\)値と\(Q^*\)値の両方とも自由度が\((\ell - K)\)\(\chi^2\)分布に従います。ここで\(K\)はモデル中のパラメータ数です。(モデルの残差ではなく)原データから計算する場合は、\(K=0\)とします。

Google株価の例では、ナイーブ法にはパラメータがないので、このケースでも\(K=0\)とします。以下のコードで、lag\(=\ell\)dof\(=K\)です。

aug %>% features(.innov, box_pierce, lag = 10, dof = 0)
#> # A tibble: 1 × 4
#>   Symbol .model       bp_stat bp_pvalue
#>   <chr>  <chr>          <dbl>     <dbl>
#> 1 GOOG   NAIVE(Close)    7.74     0.654

aug %>% features(.innov, ljung_box, lag = 10, dof = 0)
#> # A tibble: 1 × 4
#>   Symbol .model       lb_stat lb_pvalue
#>   <chr>  <chr>          <dbl>     <dbl>
#> 1 GOOG   NAIVE(Close)    7.91     0.637

\(Q\)値と\(Q^*\)値のどちらでも、結果は有意ではありません(つまり、\(p\)値が比較的大きい)でした。こうして、残差はホワイトノイズ系列と区別できないと結論付けることができます。

Googleの日次株価終値を予測するのに適しているかもしれない、もう一つの単純なアプローチはドリフト法です。tidy()関数は、このモデルの唯一のパラメータ、ドリフト係数、の推計値を示してくれます。この係数は、過去のデータから観測された平均日次変化幅です。

fit <- google_2015 %>% model(RW(Close ~ drift()))
tidy(fit)
#> # A tibble: 1 × 7
#>   Symbol .model              term  estimate std.error statistic p.value
#>   <chr>  <chr>               <chr>    <dbl>     <dbl>     <dbl>   <dbl>
#> 1 GOOG   RW(Close ~ drift()) b        0.944     0.705      1.34   0.182

Ljung-Box検定を適用する際、推計されたパラメータが1つあったので、\(K=1\)とします。

augment(fit) %>% features(.innov, ljung_box, lag=10, dof=1)
#> # A tibble: 1 × 4
#>   Symbol .model              lb_stat lb_pvalue
#>   <chr>  <chr>                 <dbl>     <dbl>
#> 1 GOOG   RW(Close ~ drift())    7.91     0.543

ナイーブ法と同様、ドリフト法からの残差もホワイトノイズ系列と区別できません。