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

以下の本を読みます。

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

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

Amazonで詳しく見る by G-Tools

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

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

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

えーと、前回やったのは、時系列分析っていうのは時系列モデルを推定することで、時系列データには周期性とかの構造があって、第2部が具体的なモデル推定の手順か（前回の最後にジュンがちょっとしゃべっちゃってたけど…）。最初はまず定常過程の説明からで、平均と分散が時刻によらず一定で、自己相関も時点によらず時間差にのみ依存するのが定常過程で、そうじゃないのが非定常過程ね。でも非定常過程でも定常過程に変換できることがあるんだ。変換できる例が単位根過程？なんか変な名前…。

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

単位根過程は1階差分が定常過程になるような非定常過程のことですね。そのような名前である理由まではこの本には触れられていないようですが、自己相関を特徴付けるある方程式が1を解にもつことに由来していて…その方程式については後の方で出てくるようですね。

ふーん？ …なあジュン、38ページの (2-6) 式って「定数＋ホワイトノイズ」なの？「定数＋ホワイトノイズ－1時点前のホワイトノイズ」にみえるんだけど、「ホワイトノイズ－1時点前のホワイトノイズ」は「ホワイトノイズ」ってこと？

確認すればいいでしょう。「ホワイトノイズ－1時点前のホワイトノイズ」の平均が0なのは自明、分散は ${\rm Var}(\varepsilon_t - \varepsilon_{t-1}) = {\rm E} \bigl( (\varepsilon_t - \varepsilon_{t-1})^2\bigr) = {\rm E}(\varepsilon_{t}^2) + {\rm E}(\varepsilon_{t-1} ^2)= 2 \sigma^2$ なので時刻によらず一定ですね。1次の自己共分散は ${\rm Cov}(\varepsilon_t - \varepsilon_{t-1}, \varepsilon_{t-1} - \varepsilon_{t-2}) = {\rm E} \bigl( (\varepsilon_t - \varepsilon_{t-1})(\varepsilon_{t-1} -\varepsilon_{t-2})\bigr) = {\rm E}(\varepsilon_{t-1}^2 )=\sigma^2$ なので1次のみゼロではありませんが、時点によらず時間差のみに依存するので結局定常過程ですね。ならいいじゃないですか。

でもちょっと気になるじゃん。サポートページにコメントしといた…ら返信もらえた！

何やってるんですか…。

42ページからはARIMAモデル？まずはARモデル…AR(1) だったらある時点の値は1時点前の値で回帰されるってことか。これって、昨日の気温が暑かったから今日も暑い、みたいな表現ができるってことだよな。iid 系列よりこういうモデルの方が時系列モデルって感じがするかも。…あれ？ iid 系列は定常過程だったけど、ARモデルは定常過程なの？ Box-Jenkins 法ではまず定常過程にしてから色々なモデルを適用するんだよな。このARモデルもそんな風にしてから適用するモデルの例ってことは、定常過程ってこと？

44ページに、AR(1) の特別な場合について何て書いてあります？

えっと待って…係数が1のときランダムウォーク？ランダムウォークって AR(1) の特別な場合なんだ。ランダムウォークは非定常過程だから、あれ、ARモデルは非定常過程ってこと？？

とは限りません。一般の AR(1) の期待値、分散、自己共分散を考えてみましょう。AR(1) の一般形は $y_t = c + \phi_1 y_{t-1} + \varepsilon_t$ です。両辺の期待値をとると、 ${\rm E}(y_t) = c + \phi_1 {\rm E}(y_{t-1})$ になり、両辺の分散を取ると ${\rm Var}(y_t) = \phi_1 ^2 {\rm Var}(y_{t-1}) + \sigma ^2$ になりますね。もしこれが定常過程であれば、 ${\rm E}(y_t) = {\rm E}(y_{t-1}) = \mu, \; {\rm Var}(y_t) = {\rm Var}(y_{t-1}) = \gamma$ ですが、 $|\phi_1| < 1$ のとき $\mu = c / (1 - \phi_1), \; \gamma =\sigma^2 / (1 - \phi_1^2)$ はこれを満たしますね。このとき自己共分散はどうでしょうか。 $k$ 次自己共分散は ${\rm Cor}(y_t, y_{t-k}) = {\rm Cor}(\phi_1 y_{t-1} + \varepsilon_t, y_{t-k}) = \phi_1 {\rm Cor}(y_{t-1}, y_{t-k}) = \phi_1^2 {\rm Cor}(y_{t-2}, y_{t-k}) = \cdots = \phi_1 ^k \gamma$ となります。これは時点に依存しませんので、定常過程です。つまり、AR(1) は $|\phi_1| < 1$ ならば定常過程になりえます。逆に $|\phi_1| < 1$ ならば期待値や分散は上の値に収束し、結局必ず定常過程になるのですが。

ARモデルが定常か非定常かは係数によるってことか…なんか面倒だな。定常データのモデリングにつかうつもりなら、もう $|\phi_1| < 1$ ってことに最初から決めておけばよくない？

1次のARモデル AR(1) の場合はそれでよいんですが…2次以上になると係数への制約をそう簡単に書き下せないんです。先ほどの AR(1) の $k$ 次自己共分散は、 $k-1$ 次自己共分散と ${\rm Cor}(y_t, y_{t-k}) = \phi_1 {\rm Cor}(y_{t-1}, y_{t-k})$ という関係が成り立っていましたよね。もっと一般に、p次のARモデル AR(p) の $k$ 次自己相関 $\rho_k$ について、 $\rho_k = \phi_1 \rho_{k-1} + \phi_2 \rho_{k-2} + \cdots + \phi_p \rho_{k-p}$ が成り立つんです。この式は、 $\rho_k$ からはじめて順に $\rho_{k-1}, \, \rho_{k-2}, \, \cdots, \rho_{0}=1$ まで自己相関を紡いでいく漸化式と見做すことができるでしょう。この漸化式の特性多項式は $f(z) = 1- \phi_1 z - \cdots \phi_p z^p$ です。同じ式が本の50ページにありますね。 $\rho_k$ はこの特性多項式の $p$ 個の根を $-k$ 乗したものの線形和で書けるんです。漸化式は数学Bでやりましたよね？数列 $a_n$ の隣接3項間漸化式を特性多項式の根を用いて解くときに、 $a_n$ の一般形に根 $\alpha$ の $n$ 乗があらわれたのを思い出してください。ここでは（添え字の向きがややこしいですが） $\rho_{0} \propto \alpha^k \rho_{k}$ かつ $\rho_0 = 1$ なので $\rho_k \propto \alpha^{-k}$ ということです。ということは、もし根 $\alpha$ の中に絶対値が1以下のものが混じっていたら、 $k$ が大きくなるほど $\alpha ^{-k}$ が振動か爆発し、ARモデルは不安定になってしまいます。というわけで、特性多項式のすべての根の絶対値が1より大きいことが、ARモデルが定常である条件なんです。

…よくわかんないけど、ARモデルが定常になる条件が面倒なのはわかったよ…。

ちなみに、37ページで出てきた「単位根過程」の名前はこのARモデルの特性多項式の「根」に由来します。単位根過程はARモデルで表現したとき特性多項式の根の1つが1になるんです。これは定常過程の1階和分を取って特性多項式を確かめてみるといいと思います。例えば―

あー！それは置いとこう！なんか全然話進んでないし！本に戻ると、ARモデルとは別にMAモデルっていうのもあるって。これは…過去 q 時点のノイズで現時点の値を回帰するってことか。…でも、ARモデルも過去の値を通して過去のノイズが含まれているよな。どう違うんだろう？

…すみません、確かに脇道にそれ過ぎましたね。ARモデルとMAモデルの関係は49ページにありますね。AR(1) は実質 MA(∞) だと。だから、過去の実現値そのものと関係がありそうな時系列データなら、ARモデルが向いているでしょう。無理やり MA(q) モデルで推定すると、q を非常に大きくしなければならないでしょうから。逆にMAモデルをARモデルで表現するときは…この話は50ページにありますね。MAモデルは（同じ期待値、分散、自己相関をもつ過程が複数存在するのですが、その内「反転可能」といわれる1つの過程は）AR(∞) に書き直すことができます。なので、MAモデルから生成されたデータを無理やり AR(p) モデルで推定しても p が大きくなってしまい、よくないでしょう。ARモデルとMAモデルはどちらかがどちらかを兼ねるというわけにはなっていないんですね。

なるほど。47ページにはARモデルとMAモデルのコレログラムがある…これは、これ見てどう思えばいいの？

モデルを推定するのに、対象時系列データの自己相関や偏自己相関のコレログラムを眺めてモデル選択するのが伝統的な手順だったんです。しかしこの本は、そう教えていませんね。そのようなコレログラムの観察は、計算資源が乏しかった頃の「古い同定の手順」だと60ページに。ただ、伝統的な手順にならえば…例えばハヤトが分析したい対象時系列データの標本自己相関をコレログラムにプロットしたところ、2次以降の自己相関がほとんどゼロだったとします。ARモデルとMAモデルのどちらを適用したいですか？

自己相関？って偏自己相関じゃない方だよな？ 47ページだと ACF っていうグラフの方…がすぐにゼロになっているのは…MAモデルの方？

ええ、逆に偏自己相関の方がすぐにゼロになっていたらARモデルを適用、という具合です。MA(q) は q+1 次以降自己相関をもちません。AR(p) は p+1 次以降偏自己相関をもちません。しかし、MA(q) は AR(∞) なので無限に偏自己相関をもち、AR(p) は MA(∞) なので無限の自己相関をもつんですね。

ARモデルとMAモデルってそういう違いがあるのか。それで、両方を組み合わせたのがARMAモデルで、ARIMAモデルというのは、何度か差分を取ってからARMAモデルを適用するってことか。54ページからのSARIMAモデルは、季節成分を考慮したARMAモデル…え、(2-39) 式って何これどうなってんの？？

落ち着いてください。例えば SARIMA(0,0,1)(0,0,1)[12] を考えてみましょう。このとき、(2-39) 式は、 $y_t = \theta(B) \Theta(B) \varepsilon_t = (1+\theta_1 B)(1+\Theta_1 B^{12}) \varepsilon_t = (1 + \theta_1 B + \Theta_1 B^{12} + \theta_1 \Theta_1 B^{13} ) \varepsilon_t$ ですから、ラグ演算子 $B$ を用いずに書けば、 $y_t = \varepsilon_t + \theta_1 \varepsilon_{t-1} + \Theta_1 \varepsilon_{t-12} + \theta_1 \Theta_1 \varepsilon_{t-13}$ ということです。

なんだそういうことか…じゃあ、SARIMA(0,0,1)(0,0,1)[12] は、1時点前のノイズに依存していてその係数が $\theta_1$ で、12時点前のノイズにも依存していてその係数が $\Theta_1$ で、13時点前のノイズにも依存していてその係数が $\theta_1 \Theta_1$ なのか…。なんか13時点前って気持ち悪くない？「1時点前」はさっきのノイズへの依存って感じで、「12時点前」は昨年のこの月のノイズへの依存って感じでわかりやすいけど、「13時点前」ってさあ…。

昨年のこの月の影響を織り込み済みのノイズ「 $\varepsilon_t + \Theta_1 \varepsilon_{t-12}$ 」を新たに「 $\tilde{\varepsilon_t}$ 」とみなしてこれで MA(1) を適用するって考えればいいんじゃないでしょうか。

あー確かに。SARIMAモデルって（月ごとに値がある周期12の時系列だったら）こんな感じなのかな？

必要ならデータ全体を前年同期との（D階）差分にする。
必要ならデータ全体を前月との（d階）差分にする。
同じ月の系列「2000年1月の値、2001年1月の値、2002年1月の値、…」を説明する ARMA(P, Q) モデルを決める（但し、どの月も同一のモデル）。
さっきの ARMA(P, Q) での過去のノイズへの依存も含めたノイズを新しい「現時点のノイズ」とみなし、過去の値への依存も含めた値を新しい「現時点の値」とみなし、ARMA(p, q) を適用する。

でもこれってなんかパラメータ多いけど、同じ期待値や自己相関をもつ過程が1つに決まるのかな…。

SARIMA が実装されているパッケージでは上手くやっているんじゃないでしょうか…。

第2部の4章の最後は、ARIMAXモデル…これはイベントの影響を考慮することができるモデルで…ダミー変数って何？過去の値とか過去のノイズだけじゃなくて、こういうのにも依存していいの？

そうですね、 $y_t$ が各日のドーナツ屋の売上高だったとしたら、よりよく推定するために、売上高以外にも日ごとの情報を色々追加してしまおうということですね。 $x_{k, t}$ は「その日が平日か休日か」「その日が晴れか雨か」「その日に近くでコンサートがあるか」などがありえると思います。しかし、「平日」も「休日」も数字ではありません。そこで仮に「平日＝0」「休日＝1」などと適当な数字を割り当ててしまいます。これがダミー変数です。

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