隼時系列本：ノート5 - クッキーの日記

以下の本を読みます。

時系列分析と状態空間モデルの基礎: RとStanで学ぶ理論と実装
馬場真哉

プレアデス出版 2018-02-14
売り上げランキング : 7742

Amazonで詳しく見る by G-Tools

※ 以下、キャラクターが会話します。原作とは関係ありません。上の本からそれる話も多いです。誤りがあればご指摘ください。

前回：ノート4 ／次回：まだ

f:id:cookie-box:20180305231302p:plain:w60

モデルを選択したら残差の自己相関や正規性をチェックするって話だったけど、本に紹介されている他にはどんな方法があるんだろう？なんでこの方法なのかって気になるし…。

f:id:cookie-box:20180305232608p:plain:w60

ちょっと調べてみましょうか。まず自己相関の有無の検定から。主に Wikipedia を参考にしました。
上2つの検定は「 $1 \sim h$ 次の自己相関が全て0」を帰無仮説としていて、このような検定を総称してかばん検定（portmanteau test）というようですね。なぜかばんというのかはよくわからなかったんですが…。

Box–Pierce test （Ljung-Box test 原論文: https://www.researchgate.net/publication/246995234_On_a_Measure_of_Lack_of_Fit_in_Time_Series_Models）	帰無仮説が正しいならば、サンプルサイズ $n \to \infty$ の極限で $(\rho_1, \rho_2, \cdots , \rho_h)$ が平均ゼロ、分散 ${\rm Var}(\rho_k)=(n-k)/(n(n+2))$ 、共分散ゼロの多変量正規分布にしたがうことに着目すると、検定統計量を $Q=n \sum \hat{\rho}_k^2$ とすればこれが近似的にカイ2乗分布にしたがうので、この $Q$ が大きいときに帰無仮説を棄却しよう（何か自己相関があると疑おう）というものです。ただ、これは近似の仕方の関係で、自己相関があるにもかかわらず $Q$ が小さくなってしまう場合があったようなんです。これを改善したのが Ljung-Box test です。
Ljung–Box test （同上）	Box–Pierce test の検定統計量を $Q=n(n+2)\sum(n-k)^{-1} \hat{\rho}_k^2$ に修正したものです。本でメインに紹介されていますね。
Breusch–Godfrey test （Breusch–Godfrey test - Wikipedia）	これは時系列モデルの残差というより、線形回帰の残差に対して自己相関の有無を判定する検定のようなんですが…線形回帰の残差系列 $u_t$ を AR(p) でフィッティングしたとき、1～p次の自己相関が全てゼロならばこのフィッティングの決定係数がカイ2乗分布にしたがうことを利用するようです。
Durbin–Watson test （Durbin–Watson statistic - Wikipedia）	線形回帰の残差に1次の自己相関があるかどうかを検定します。この本では125ページで出てきますね。上の3つの検定が色々な次数の自己相関がゼロであることを一度に検定する一方で、1次の自己相関のみを検定というのは随分限定的なようですが、上の3つの検定が対応できない、誤差の分散が不均一な場合などにも対応できるようです。

次に正規性の検定です。これらの検定の帰無仮説は「データが正規分布にしたがう」ですね。正規性の検定は古くより数多く考案されており、初期には尖度や歪度が着目されていたようです。
Normality test - Wikipedia

D'Agostino's K-squared test	データの尖度と歪度がどれだけ正規分布のそれより逸脱しているかを $K^2$ 統計量としてまとめて検定するものです。 $K^2$ 統計量は近似的にカイ2乗分布にしたがいます。
Jarque–Bera test	本で紹介されていますね。考え方が D'Agostino's K-squared test と同じですが、検定統計量の近似がより洗練されているようです。
Anderson–Darling test （http://www.hep.caltech.edu/~fcp/statistics/hypothesisTest/PoissonConsistency/AndersonDarling1954.pdf）	この検定と Kolmogorov–Smirnov test は、正規性の検定というよりはもっと一般的に、データがある特定の確率分布から生成されたものかを検定します。もしデータがある確率分布 $f(x)$ から生成されているならば、そのデータのaパーセンタイル点と $f(x)$ のaパーセンタイル点は近いだろうというのはもっともなアイデアですよね。昇順にソート済のデータ $\{ x_i \}$ について、「そのデータがデータ中で小さい方から何%か」を横軸に、「その値は $f(x)$ の何パーセンタイル点か」を縦軸にプロットすると、 $(0, 0)$ から $(1, 1)$ まで直線上に並ぶだろうということです。この直線からのずれの2乗を足し合わせたものが Anderson–Darling test の検定統計量に相当しますが、帰無仮説の下でこのずれの分散は $F(x)(1-F(x))$ となり、分布の裾ほど小さくなります。なので、検定統計量は各点の誤差をこの分散の逆数で重み付けしています。この重み付けは、サンプリング誤差を均しているという見方もできますし、分布の裾のずれほど重くみているという見方もできます。 Anderson–Darling test によって正規性を検定するときは、平均／分散が既知／未知かでバリエーションがあることに注意ですね。
Kolmogorov–Smirnov test	Anderson–Darling test に先立つ手法ですね。この検定統計量は単に経験分布の累積分布関数と検定したい分布の累積分布関数の差のL∞-ノルムですね。なので、「2つのデータセットが同じ分布から生成されているか」という検定にもそのまま利用できます。2つの経験累積分布の差に置き換えるだけですので。これは画一的で便利な検定である一方、正規性の検定としての検定力は Anderson–Darling test や Shapiro-Wilk test に劣るそうですので、これらの検定が利用できる場面ではそちらを利用した方がよさそうですね。
Shapiro-Wilk test	これは分布全体の形というよりもそれぞれの値の位置に着目していますね。サンプルサイズが $n$ だったとき、小さい方から $i$ 番目のデータが、正規分布から生成された $n$ 個のサンプルの $i$ 番目に小さいデータの期待値に近いか、といったイメージです。2011年のある研究によると、正規性の検定の中で最もよい検出力をもつとか。検出力がよいということはつまり、データが正規分布から生成されたものでないときに、ちゃんと帰無仮説を棄却できる確率が高いということですね。一方で、データの中にタイ（複数のデータが全く同じ値となる）が多い場合に弱いなど、弱点もあります。小さい方から $i$ 番目のデータの位置、などに着目しているので、タイには弱そうですよね。

それを聞くと自己相関の有無の検定はやっぱり Ljung–Box test がスタンダードなのか。正規性の検定もたくさんあるな。同じことを確認しようとするのに、色々な考え方があるってのが面白いけど面倒かも。

統計的検定では、帰無仮説を棄却できるかどうか、つまり、目の前のデータがあるパラメータをもつ確率分布から生成されていることを否定できるかどうかを判定しようとします。「帰無仮説が正しくなさそうな度合い」を何で測るかによって色々な考え方が生じるでしょう。

でも目の前のデータの背後に確率分布があるっていうのが人間の空想なのに、この確率分布じゃないと否定されても、確率分布にしてみればとばっちりだよな。

え？すみませんが、棄却された確率分布の側がどう感じられているかは僕にはちょっと…。

72ページからは時系列分析を実際にやってみようって話か…といっても今パソコンないしな…。

ノートPCならありますよ。

なんで学校にノートPC持ってきてんの！？

ちょっとした編曲作業に使おうと思いまして。学校でアイデアが浮かぶこともあるかもしれませんし。

いや確かに曲のアイデアは逃がしたくないけどさ…。

でもこれ再現実行しても本の内容なぞるだけになっちゃいますし、とばして第3部を読みましょうか。見せかけの回帰を別の具体例でやってみましょう。

そうそう、第3部に「見せかけの回帰」ってあったけど、何が見せかけなんだ？

見てみるのが早いかと。互いに全く関係ない2つの時系列モデル

$x_t = 0.8 x_{t-1} + 0.1 x_{t-2} + 0.2 \varepsilon _{x, t-1} + \varepsilon _{x, t}$
$y_t = 0.6 y_{t-1} + 0.2 y_{t-2} + 0.1 \varepsilon _{y, t-1} + \varepsilon _{y, t}$

を用意しましょう。たまたまどちらも ARMA(2,1) にしてみました。では、これらのモデルから 100 点ずつデータを生成し $(x_1, \, y_1) , \cdots, (x_{100}, \, y_{100})$ を得て、線形回帰するとどうなるでしょう。

どうなるって…互いに全く関係ないデータに線形回帰しても意味ないだろ。

実際にやってみましょう。データを生成して R の lm() で線形回帰分析します。lm() はデータに線型モデルを仮定すると、モデル係数の最小2乗推定量や残差、決定係数やF値などを計算してくれる関数ですね。以下では $y_t = \alpha + \beta x_t + \varepsilon_t \; \bigl( \, \varepsilon_t$ は互いに独立に $N(0, \sigma^2)$ にしたがう $\bigr)$ を仮定することになります。

set.seed(1234)
x <- arima.sim(n=100, model=list(ar=c(0.8, 0.1), ma=c(0.2)))
y <- arima.sim(n=100, model=list(ar=c(0.6, 0.2), ma=c(0.1)))
plot(as.numeric(x), as.numeric(y))
print(summary(lm(y~x)))

実行結果は以下です。傾きが約 -0.4、切片が約 1 となりました。傾きのp値も切片のp値も小さく、傾きも切片も有意であるという結果になっています。つまり、傾きも切片もゼロとは思えないということですね。

Call:
lm(formula = y ~ x)

Residuals:
     Min      1Q  Median      3Q     Max 
 -3.7131 -0.9341 -0.1749  0.9870  3.7385 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  1.06548    0.16071   6.630 1.84e-09 ***
x           -0.43461    0.08694  -4.999 2.52e-06 ***

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.454 on 98 degrees of freedom
Multiple R-squared:  0.2032,    Adjusted R-squared:  0.1951 
F-statistic: 24.99 on 1 and 98 DF,  p-value: 2.523e-06

え…？ $x$ も $y$ も互いに関係ないモデルから生成されたのに、線型モデルをあてはめると傾きや切片がゼロとは思えないって…どうして？？

これが見せかけの回帰ですね。

えー何が起こってるのか全然わかんない！これって、p値的にはそのモデルで表される線形関係があるはずって結果なんだよな？じゃあもう見せかけの回帰じゃなくてリアル回帰じゃん。

まだ気付きませんか、ハヤト。バンドのリーダーなんだからしっかりしてください。

バンド関係ない！

線形回帰分析の係数のt検定の検定統計量がt分布にしたがうのは、あくまで「残差が互いに独立で平均ゼロの正規分布にしたがう」という仮定の下での話です。 $x_t$ と $y_t$ が自己回帰モデルから生成される場合にこの仮定は成り立つでしょうか。簡単のために、どちらもランダムウォークとして計算してみます。

$x_t = x_{t-1} + \varepsilon_{x, t}$
$y_t = y_{t-1} + \varepsilon_{y, t}$
$y_t = \alpha + \beta x_t + \varepsilon_t \; \Rightarrow \; \varepsilon_t = y_t - \alpha - \beta x_t$
$\varepsilon_t = y_{t-1} + \varepsilon_{y, t} - \alpha - \beta (x_{t-1} + \varepsilon_{x, t}) = \varepsilon_{t-1} + \varepsilon_{y, t} - \beta \varepsilon_{x, t}$

これより、残差 $\varepsilon_t$ もランダムウォークになることがわかります。残差がランダムウォークするのに、最小2乗推定で回帰係数を推定するべくもありません。ましてや、検定統計量がt分布にしたがうはずがないんです。これを考慮せずに lm() で線形回帰分析しても、その結果はまやかしであり、紛れもなく「見せかけ」の回帰なんです。

t検定が利用できる状況が整っていないのにt検定をしたのがいけなかったのか。じゃあどうすればいいんだ？

残差に自己相関があるのにt検定をしてしまうのが悪いんですから、残差に自己相関が残らないように階差系列をとるなどデータの方を加工するか、回帰モデルの方で自己相関を表現できるように修正するか、あるいは残差に自己相関がある可能性を踏まえた回帰分析の方法を導入する必要があるでしょう。先ほどのデータに一般化最小2乗法（GLS）を適用してみましょう。まず本で紹介されている prais パッケージの関数から。

library(prais)
print(prais.winsten(y~x, data.frame(x=x, y=y), iter=1))

結果はこうです。もはや傾きのp値が小さくありませんね。

Call:
lm(formula = fo)

Residuals:
     Min      1Q  Median      3Q     Max 
 -3.3421 -0.5651  0.0661  0.5300  3.0788 

Coefficients:
          Estimate Std. Error t value Pr(>|t|)  
Intercept  0.76130    0.32288   2.358   0.0204 *
x         -0.04359    0.10781  -0.404   0.6869  

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.9879 on 98 degrees of freedom
Multiple R-squared:  0.05423,   Adjusted R-squared:  0.03493 
F-statistic:  2.81 on 2 and 98 DF,  p-value: 0.06508

      Rho Rho.t.statistic Iterations
 0.688919         9.44636          1

他にも GLS は色々なパッケージがあるようなので、少し調べて出てきたのもやってみます。rms パッケージの Gls() ですが、どのような自己相関をもつのか引数 correlation に明示的に与えるようです。この引数を省略すると回帰は有意となってしまいますが、手抜きで corAR1() と指定してみると、回帰は有意ではなくなりました。 Gls function | R Documentation

library(rms)
print(Gls(y~x, data.frame(x=x, y=y), correlation=corAR1()))

次のようなページも見つけました。
一般化最小二乗法と識別問題 - hnami.net_Pukiwiki

ふーん、とりあえず残差に自己相関がある場合の回帰にも対処できる方法があるならよかった。…138ページからの、共和分って何？

簡単な例だと、 $x_t$ と $y_t$ がそれぞれ単位根過程だったとして、 $x_t$ と $y_t$ を適当に線形結合した過程が定常になるなら、 $x_t$ と $y_t$ は共和分しているといいます。

なるほど…だから何？

単位根過程は差分系列をとって定常過程にしてから解析するのがセオリーです。しかし、 $x_t$ と $y_t$ が共和分している場合、先にそれぞれ差分系列をとってしまうと、 $x_t$ と $y_t$ の関係が失われてしまうということですね。下の沖本氏の本の132ページから、共和分を踏まえなければ導出できない経済事象の具体例が3つ挙げられています。
経済・ファイナンスデータの計量時系列分析 (統計ライブラリー) | 沖本竜義 |本 | 通販 | Amazon
つまり、1つ目の例は消費系列と所得系列の関係を調べるのに、これらが共和分していることそのものが示唆になっているのですが、先に差分系列をとってしまうとこの関係が見えなくなってしまいます。共和分している系列どうしに回帰モデルをあてはめると、残差が定常になりますので、共和分が疑われるときはそれを検定することになります。本の142ページをみると、この用途専用の検定をつかわなければならないようですね。

（ノート6があれば）つづく