あかり描像のブログ

思ったことや学習記録を適当に書いていきます。お気軽にコメントください

トップ > 統計検定 > 【統計検定】2019-11 理工学4 (社会科学3)
最終更新日

【統計検定】2019-11 理工学4 (社会科学3)

統計検定 数学 算数

時系列分析の勉強をしようとちょっと触ってみたら行列の算数に蹂躙されて酷い目に遭ったので戦闘の記録をまとめます。
解答書ではありません、あくまで戦闘の記録です。

英語で検索を掛けた方が情報が見つかるので、記事中にちょくちょく英単語を載せてます。

目次

Notation

  • ベクトルや行列は太字で表しません。
  • ベクトル x や行列 A の転置は x^\prime, A^\prime のようにプライムで表します。
  • 確率変数 X が分布 D に従うことを X\sim D と書きます。
    • iid ... 独立同一分布に従う (independently and identically distributed)
  • \mathbb{E}[\cdot] ... 期待値, 期待値ベクトル
  • \mathbb{V}[\cdot] ... 分散, 分散共分散行列 (covariance matrix)
  • \mathrm{Cov}(\cdot, \cdot) ... 共分散 (covariance)
  • \chi^2 (n) ... 自由度 n のカイ2乗分布 (chi-squared with n degrees of freedom)
    • \chi_\nu^2 (n) ... 自由度 n のカイ2乗分布の上側 100\nu % 点

問題

時系列データ (time series data) \left\{ X_t \right\}_{t\in\mathbb{Z}} が次のように線形な \mathrm{AR}(1) 過程 (autoregressive process) に従うとする.

X_t = \phi X_{t-1} + \epsilon_t \quad (t \in \mathbb{Z})

ここで \epsilon_t \sim N(0, \sigma^2) iid であり, (弱) 定常条件 (weak stationarity) |\phi| < 1 をみたすものとする.

(1)

X = (X_1, \dots, X_n)^\prime の自己共分散行列 (autocovariance) T= \left\{ \tau_{ij} \right\} , 自己相関行列 (autocorrelation) R = \left\{ \rho_{ij} \right\} がそれぞれ次のようになることを示せ。

  • \tau_{ij} = \dfrac{\sigma^2}{1-\phi^2} \phi^{|i-j|}
  • \rho_{ij} = \phi^{|i-j|}

(2)

n 次対称行列 A=\left\{ a_{ij} \right\} を次のように定める。
a_{ij} = \begin{cases} 1 & (i = j = 1, \ i = j = n)\\ 1 + \phi^2 & (i = j = 2, \ \dots, n-1)\\ -\phi & (|i-j|=1)\\ 0 & (|i-j| \geq 2) \end{cases}

この時、

(3)

x = (x_1, \dots, x_n)^\prime とする。2次形式 Q_A = x^\prime A x を x_1, \dots, x_n で書き下し、 A が |\phi| < 1 の時に正定値になることを示せ。すなわち、 \forall x \neq 0 に対し Q_A > 0 を示せ。

(4)

行列 A の (1,1) 成分を 1 から \phi^2 に変えた行列を B とする。2次形式 Q_B = x^\prime B x が非負定値になること、すなわち \forall x に対して Q_B が \phi によらず非負になることを示せ。また、 Q_B = 0 となる x \neq 0 はどのようなベクトルか。

(5)

モデルにおける自己回帰係数 \phi が既知であり、誤差分散 \sigma^2 が未知であるとき、 \sigma^2 の 95 % 信頼区間はどうなるか。

戦記

|\phi| < 1 のことを弱定常性 (weak stationarity) と呼んでいますが、実際の定義は違います。

時系列データ \{ x_t \} が弱定常 (weak stationary), 若しくは共分散定常 (covariance stationary), はては二次定常 (second order stationary) であるとは、 \forall t, k に対して、期待値 \mathbb{E}[x_t] と自己共分散 \mathrm{Cov} (x_t, x_{t-k}) が t に依らないことをいう。言い換えると、次のように書けることをいう。

  • \mathbb{E}[x_t] = \mu < \infty
  • \mathrm{Cov} (x_t, x_{t-k}) = \gamma_{k} < \infty

これから、分散 \mathbb{V}[x_t] = \mathrm{Cov} (x_t, x_t) = \gamma_0 も constant になることは直ちに分かります。ちゃっかり有限であることも要請していますが、英語版 Wikipedia によると 2 次モーメント \mathbb{E}[x_t^2] の有限性のみ要求すればよいらしいです。一体何が起きてるんですかね。

これを仮定すると、 \mathrm{AR}(1) が有限の期待値及び分散を持つためには |\phi| < 1 が必要になることはすぐに出ます。しかし、逆に |\phi| < 1 を仮定してこれが出てくるかどうかは 知りません。どうやら同値らしいので出るとは思うのですが、私はそのフォローの仕方が分かりません。天下の沖本やハミルトンでもこの辺がボカされている感じがします。 特性方程式って何やねん。
分かる方、教えてください m(_ _)m
私のお勉強が進んできたらまたここに書きます。

また、 \mathrm{AR}(1) の式を見ただけでは「初期値ねぇじゃんw」という気持ちになりますが、弱定常性を仮定すれば初期値を考えなくても良いことになります。なるそうです。

ごめんなさい。弱定常性、ガチで分からん。

そんなところなので、問題文では |\phi| < 1 しか与えられていないですが、 \mathbb{E} \left[ X_t \right] や \mathbb{V} \left[ X_t \right] が constant になることや、 \mathrm{Cov} (X_t, X_{t+k}) が k のみに依存することなどは無条件に認めてしまうことにします。でも多分恐らくひょっとしたらそれでも論理は崩れていないことが聊か予想されると考えられるものと存じます。本当にごめんなさい。

(1)

時系列分析と言ったらこれ、というくらい very special one pattern な問題。
とりあえず期待値, 分散, 自己共分散, より一般にユール・ウォーカー方程式. ところによってはスペクトル密度関数 (spectral density) を計算させられることでしょう (2021-06 準1級)。

やな弱定常性はぜーんぶゴミ箱に捨てちゃえ。

※無限に解説書があるので省きます orz...
パワーがある時に LaTeX を打ちます。

(2)

行列の算数です。おもしろい問題なのに捨て問たらしめている元凶。この記事を書こうと思った理由。
色々やり方を考えましたが、結局注意深く場合分けするくらいしかないと思います。個人的にはクロネッカーのデルタ等を用いて機械的にガチャガチャするのがしっくりくるのでその方針で書きます。
行列式の計算も色々やり方があると思いますが、行基本変形で計算するのがシンプルだと思います。

過去問の答えは n=4 の場合だけ書いてあるけど、正直これズルだと思う。


a_{ij} はクロネッカーのデルタを用いて次のように書ける。

\begin{aligned} a_{ij} &= (1+\phi^2) \delta_{i,j} - \phi^2 \delta_{i,j} (\delta_{i, 1} + \delta_{i, n}) - \phi (\delta_{i-1, j} + \delta_{i, j-1}) \end{aligned}

これから、

\begin{aligned} &\left( \frac{1}{\sigma^2} TA \right)_{ij}\\ &= \frac{1}{\sigma^2} \sum_{l=1}^n \tau_{il} a_{lj} = \frac{1}{1-\phi^2} \sum_{l=1}^n \phi^{|i-l|} a_{lj}\\ &= \frac{1}{1-\phi^2} \sum_{l=1}^n \phi^{|i-l|} \left[ (1+\phi^2) \delta_{l,j} - \phi^2 \delta_{l,j} (\delta_{l, 1} + \delta_{l, n}) - \phi (\delta_{l-1, j} + \delta_{l, j-1}) \right] \\ &= \frac{1}{1-\phi^2} \left[ \phi^{|i-j|} \left\{ (1+\phi^2) -\phi^2 (\delta_{j, 1} + \delta_{j, n}) \right\} \right.\\ &\left.\qquad\qquad\qquad - \phi\cdot \phi^{|i-j-1|} (1 - \delta_{j, n}) - \phi\cdot\phi^{|i-j+1|} (1-\delta_{j, 1}) \right] \end{aligned}

最後の = では、 \delta_{l-1, j} を l について和を取ると l = j+1 の項だけが残るが、 j=n の場合はこれを満たし得ないので、 \sum_{l=1}^n \phi^{|i-l|} \delta_{l-1,j} = \phi^{|i-(j+1)|} (1-\delta_{j, n}) であることを用いている。 \delta_{l, j-1} についても同様である。

注意深く場合分けすればこれは \delta_{ij} になる、とでもできれば締まりが良いのだが、あんまり自明ではないので一応確認する。 \phi の肩にある 3 つの絶対値の中身の正負で場合分けする。

  1. i=j の時は、計算により 1 となる。
  2. i > j の時は、 |i-j| = i-j, \ |i-j-1| = i-j-1, \ |i-j+1|=i-j+1 であることから計算する。 \delta_{j, n} の項だけ残るが、 i > j では j = n にはならないので 0 である。
  3. i < j の時も同様に、 |i-j| = j-i, \ |i-j-1| = j-i+1, \ |i-j+1|=j-i-1 から計算する。 \delta_{j, 1} の項だけ残るが、 i < j では j = 1 にはならないので 0 である。

結局、 \delta_{ij} であることが分かったので、 \frac{1}{\sigma^2} T の逆行列は A である。


次に行列式 |A| について、 r_2 + \phi r_1, \ r_3 + \phi r_2, \dots, \ r_n + \phi r_{n-1} なる行基本変形を逐次行うことで

\left| A \right| = \begin{vmatrix} 1 & -\phi & & & O \\ &1 & -\phi & & \\ & & \ddots & \ddots & \\ & & & 1 & -\phi\\ O & & & & 1-\phi^2 \\ \end{vmatrix} = 1-\phi^2.

最後に、
\left| R \right| = \left| \dfrac{1-\phi^2}{\sigma^2} T \right| = \left| (1-\phi^2) A^{-1} \right| = (1-\phi^2)^n \left| A \right|^{-1} = (1-\phi^2)^{n-1}.

(3)

何が起きているのか分かりにくいので、とりあえず n=4 の場合を書き下してみると、何やら平方完成っぽいことができることが見て取れます。

\begin{aligned} Q_A &= x_1^2 + (1+\phi^2) x_2^2 + (1+\phi^2) x_3^2 + x_4^2 - 2\phi (x_1 x_2 + x_2 x_3 + x_3 x_4)\\ &= (x_1 - \phi x_2)^2 + (x_2 -\phi x_3)^2 + (x_3 - \phi x_4)^2 + x_4^2 - \phi^2 x_4^2 \end{aligned}

ということなので、これをあたかも n で考えましたけど？って体で答案を書けば良いです。*1

\begin{aligned} Q_A &= x^\prime A x\\ &= x_1^2 + (1+\phi^2) \sum_{i=2}^{n-1} x_i^2 + x_n^2 - 2\phi \sum_{i=1}^{n-1} x_i x_{i+1}\\ &= \sum_{i=1}^{n-1} (x_i - \phi x_{i+1})^2 + (1-\phi^2) x_n^2 \end{aligned}

これは、 x\neq 0, \ |\phi| < 1 で常に正.

(4)

(3) で要領を得られれば、芋づる式にできると思います。


\begin{aligned} Q_B &= x^\prime B x\\ &= \phi^2 x_1^2 + (1+\phi^2) \sum_{i=2}^{n-1} x_i^2 + x_n^2 -2\phi \sum_{i=1}^{n-1} x_i x_{i+1}\\ &= \sum_{i=1}^{n-1} (\phi x_i - x_{i+1})^2 \geq 0. \end{aligned}

Q_B = 0 となる x\neq 0 は、 x_2 = \phi x_1, \ x_3 = \phi x_2, \dots , x_n = \phi x_{n-1}, すなわち
x^\prime = (1, \phi, \phi^2, \dots, \phi^{n-1}) \times \mathrm{const.}

(5)

一般に、 \Sigma を n 次の正定値実対称行列とした時、

X\sim N(\mu, \Sigma) \Longrightarrow (X - \mu)^\prime \Sigma^{-1} (X - \mu) \sim \chi^2(n)

という事実 (公式本 p.47 の定理2.3.2) を知っていれば瞬殺ですが、私みたいに知らない人は 1 からお勉強です。*2


とどのつまり、ここでやりたいことは多変量正規分布 (multivariate normal distribution) の標準化です。

一般に、多変量正規分布は必ずしも標準化することはできませんが、共分散行列が正定値 (positive definite) であるときは標準化ができます。どんな分布であっても (正規分布でなくても) 共分散行列はほとんど確実に半正定値 (positive semi-definite) になりますが、そこに正定値という条件が加われば標準化ができるということです。(大事なことなので 2 χ 言いました。)
一般に、多変量正規分布は共分散行列が正定値でない場合も定義することができますが、そのような場合では密度関数が存在しないので (逆も然り)、普通は正定値になるように定義します。言い換えれば、普通の人ならば多変量正規分布はいつでも標準化できると思っているということです。*3 実際、本問でも T は正定値になっているので標準化できます。

標準化から信頼区間を求めるまでの論理をまとめると大方以下のようになると思います。

  • T は正則 (regular, non-singular) なので *4、 T の半正定値性と併せて正定値です。*5
    • A は正定値なので、その逆行列 \frac{1}{\sigma^2} T も正定値、という論理でも ok です。
  • T は正定値なので、正則行列 C を用いて T = C C^\prime のように分解できます*6。
  • X\sim N(0, T) なので、 X = CY とすると Y \sim N(0, I) が示せます。*7
  • すると、 X^\prime T^{-1} X = Y^\prime Y = \sum_{i=1}^n Y_i^2 が \chi^2 (n) に従うことが分かります。*8
  • 使える分布が分かったので、あとは観測値から \chi^2 統計量を計算して信頼区間が求まります。


という訳で、答案を作ると大体次のようになると思います。過去問の解答もすっごく簡素ですが、この感じだと 2次形式が \chi^2 (n) に従うことは既知としても良さそうですね。


(2), (3) の結果から T は正定値になるので、2次形式 x^\prime T^{-1} x = x^\prime ( A / \sigma^2 ) x = Q_A / \sigma^2 は \chi^2 (n) に従う。よって、観測値 \{ x_i \}_{i=1}^n から Q_A / \sigma^2 を計算することで、 \sigma^2 の 95% 信頼区間が

\chi_{0.975}^2 (n) \leq \dfrac{Q_A}{\sigma^2} \leq \chi_{0.025}^2 (n)

すなわち

\dfrac{Q_A}{\chi_{0.025}^2 (n)} \leq \sigma^2 \leq \dfrac{Q_A}{\chi_{0.975}^2 (n)}

と求まる。

感想

定常性、何も分からん。

本番ではこれを 30 分でやらんといかんのか。無理やな！！

参考

*1:ちなみに、A は実対称行列であり、かつ (2) で計算した |A| = 1 - Φ² から |Φ| < 1 の時に |A| > 0 であるので、A が正定値であること自体は 2 次形式を計算するまでもなく明らかです。

*2:今回はまったく関係ないですが、平均 0 の多変量正規分布 N(0, Σ) に iid する n 個の確率変数ベクトル X1, X2, ..., Xn ～ N(0, Σ) iid に対し、Σi Xi Xi' が従う分布は (自由度 n の) ウィシャート分布 (Wishart distribution) と呼ばれます。n = 1 かつ Σ = 1 の時は自由度 n の χ² 分布になりますね。この分布は平均が既知である場合の多変量正規分布の精度行列 Σ⁻¹ (精度とは分散の逆数のこと) についての共役事前分布としてよく知られています。なお、共分散行列 Σ は正定値である必要はないことに注意してください。 → PRML, etc.

*3:しかし、正定値になっていない場合もそれなりに重要です。

*4:A / σ² という逆行列が見つかったので

*5:分散共分散行列は分布の種類によらず常に半正定値です。そして、半正定値行列が正定値であるための必要十分条件はそれが正則であることです。 → 参考ブログ

*6:実対称行列 T について、T が正定値であることと、 T = CC' (C: 正則行列) と分解できることは同値です。一般に C は一意には定まりませんが、C も正定値となるようにすると一意に定まり、これを T の主平方根ということがあります (この場合は C は対称行列になるので T = C² になります。)。また、 C は下三角にもできて、この場合をコレスキー分解と呼びます。コレスキー分解も一意ではないですが、C の対角成分が正の実数となるようにすると一意になります。どんな方法であれ、正定値行列 T の平方根なるものは適当に定義することができ、その中の一つを T^{1/2} や \sqrt{T} 等と書くことがあるようです。また、C は正則なので逆行列も存在し、それを T^{-1/2} のように書いたりもするそうです。なお、T が半正定値ではあるが正定値ではない場合も T = CC' のような分解は可能ですが、この場合は C を正則にできないので、 T^{-1/2} なるものが存在せず論理が破綻します。標準化の可否が正定値性と直結していることが何となく分かってきましたね。

*7:一般に、正規分布に従う確率変数の組を束ねてベクトルにしてもそれが多変量正規分布に従うとは限らないので、X ～ N(0, T) を説明するのにも案外難儀します。しかし、いったんこれを知らんふりして認めてしまうと、Y も多変量正規分布に従うことは (多変量正規分布に従う確率変数ベクトルを線形変換しても多変量正規分布に従うことから) 直ちに分かり、その期待値, 共分散行列も簡単な算数で計算できます。なお、X の平均が 0 ではなく一般に μ であるときは、X = CY + μ と変換することになります。

*8:2 級くらいの範囲の話だと思いますが、N(0, 1) に iid する n 個の確率変数の 2 乗和が χ²(n) に従うことは、畳み込み積分+帰納法、もしくは特性関数による方法で確認できます。何ならこれが χ²(n) の定義まであります。