2018-04-23

隼時系列本：ノート7

本読み

以下の本を読みます。

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

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

Amazonで詳しく見る by G-Tools

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

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

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

後半からはナツキが本読み手伝ってくれるんだな。よろしく。

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

うん。ジュンは学校も仕事も忙しいから、最後まで付き合うのは難しい。だから、後は俺に任せて…。

いや俺も学校も仕事も忙しいけどな！？ …でも、ナツキって前半読んでないんだよな？いきなり後半からで大丈夫なのか？

大丈夫、この本の前半と後半は比較的別々に読めるはず…ハヤト、前半の話をまとめられる？

この本の前半の話は…まず、時系列分析っていうのは、手元の時系列データの将来値の予測をしたりその時系列データについて何かしらの知識を得ることが目的。その目的のために何をするかっていうと、時系列データには各時点での値の背後に確率分布があると仮定して、尤もらしい確率分布列を与えるモデルを特定する。そのモデルの候補として、ARIMAやVAR、GARCHとかが紹介されてた。ただ、このモデルを同定するときの注意点として、予め対象時系列データを定常にしておくことや、あてはめたモデルで説明できない残差の正規性や自己相関の確認が要るんだったな。そのために用いる統計的検定も色々出てきた。

まとめありがとう。後半の「状態空間モデル」も、目的はだいたい一緒だよ（「状態」を明示的に扱うなりの目的も加わってくるかもしれないけど）。でも、目的のために何をするかがもっと一般的になる。時系列データには各時点での値の背後に、直接観測できない「状態」と、状態をその時点まで運んできた「システムモデル」、状態を目に見える観測値の分布に写す「観測モデル」があると考えるよ。もし「状態」が観測される変数そのもので、「観測モデル」が恒等写像だったら、前半までの話と一緒の構図になるけどね。

直接観測できない「状態」があると考えるって？

本の179ページにも魚の数の例があるけど、「状態」を導入する意義の1つは、時系列データって、「観測するときにだけ混ざってくるノイズ」があるかもしれない。例えば、俺たちが持ち回りで部室の気温を毎日記録するとして、ハヤトが記録する日だけはハヤトが温度計に近づくせいで温度計の指示値が少し高めに記録されるかもしれない。でも、そのハヤトのせいで記録値が少し大きくなった分は、翌日以降の気温に影響する成分とは思えないよね。

なんていうか、それはもう俺を記録係から外した方がいいな…。

この本の前半のモデルにも誤差項はあるけど、それはその場限りの誤差ではなくて、その影響はその後ずっと（MAモデルの場合はどこかで打ち切られるけど）続いていくものになっている。そうじゃなくて、現実に何かを測るときはその場限りの誤差ってあるはず。だったら、その場限りの誤差が付加される前の「状態」の推移を考えればいい。それが状態空間モデルの意義。それに加えて、第4部の冒頭には明示的に書いてないと思うんだけど、状態って何も「観測される変数にその場限りのノイズが付加される直前のもの」とは限らなくて、「観測される変数そのものではないけど、その推移を追うと有用な量」だったり「観測できないけど興味がある量」かもしれない。状態空間モデルならそういう設計もできる。例えば、ハルナが毎日食べるドーナツの個数は直接観測できて、ハルナの毎日の空腹度は直接観測できないけど、したがう方程式がわかっているのは後者という場合とか。

だから、ハルナの空腹度がしたがう方程式が研究で明らかにされてたらおかしいだろ！ …でもナツキ、状態空間モデルでは、その状態がしたがう方程式が予めわかってないといけないってこと？前半ではさ、ARIMA だったら ARIMA の次数と係数は AIC なり何なりを指標に計算機にゴリ押しで決めさせるって感じで、「対象時系列はこんなメカニズムをもつんじゃないか」って考えるプロセスがなかった気がするんだけど…。

ARIMA を使ってみる時点である意味「ARIMA で表現できるんじゃないか」っていう考えがあるとも思うけど…状態空間モデルでは、対象時系列に対して「こんなメカニズムをもつんじゃないか」というのがもっと明示的にある場合が多いんだろうね。もちろん、データに対する強い信念なしに、色々な次数の ARIMA を試すような用法もあると思うよ？システムモデルの係数を未知パラメータにして。182ページの最初で「このあたりの手順に関しては」って別の文献を紹介しているね。でもそういうモデル探索よりは、何かこういうシステムモデルだろうという信念が予めある程度あって、直接観測できない「状態」の推移やその分散を捉えたいって向きが強いんじゃないかな。それが180ページの「過去の知見（中略）自由にモデル化できる」「（前略）モデルの解釈が容易」ってことだと思うから。俺も普段から状態空間モデルを使うわけじゃないから、よくわからないけど…。

普段から状態空間モデルを使う高校生はいないだろ！

まあそれで、この後の中心話題は未知パラメータをどうやって推定するかだね。この本では2つの異なるアプローチのそれぞれの代表として、カルマンフィルタとHMC法を大きく取り上げてる。この辺について、ちょっとプロデューサーさんが昔書いた記事を掘り出しておくね。

カルマンフィルタに関する記事は以下など。
- 状態空間モデルのまとめ絵
- R dlm をつかってみるだけ
1時点前までの観測値を得た下での現在の状態の分布は、一時点前の状態を用いて、

だけど、これに現在の観測値が手に入った下での現在の状態はベイズの定理より

になる。もし状態空間モデルのシステムと観測が線形で、システムノイズと観測ノイズがガウシアンなら（線形・ガウス状態空間モデル）、をガウス分布とすると、もガウス分布、さらにもガウス分布、…となって、1つ手前を含む形の式になって、どんどん逐次的に計算していけるんだよね。これがカルマンフィルタだね。
HMC法に関する記事。後者の記事、収束してなさすぎだけどスクリプト合ってるのかなあ…。

もっと具体的な計算手順や、これら以外の状態空間モデルの解き方とかはまた今度にする？

うーん、183～185ページに今ナツキが言ったような2つのアプローチのパラメタ推定があるってあるけどさ、なんでその2つなのかよくわからないんだよな。要するに、真面目に解くか、乱数シミュレーションで解くかって違い？でも、逐次的かどうかも違って…真面目に同時に解くとか、乱数的に逐次的に解くとかはないの？

ないとは言い切れない、けど…真面目に全期間同時に解けることなんてまずないんじゃない？逐次的に解けるような特殊なケースだから真面目に書き下せたってだけで…。あと185ページに、HMC法ではパラメタの推定を同時にやるから逐次的にはしにくいとも書いてある。

そっか、どちらかというと、特殊なケースでは真面目に解けて、副産物的に逐次解法になっている、って感じなのかな…。あとナツキ、209ページの散漫カルマンフィルタってのは？

あまり他では聞かない用語、だけど…状態の分散の初期値をものすごく大きくしておいて、状態の期待値の初期値の影響をすぐ消滅させるテクニックを指しているみたいだね。

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

2018-04-03

隼時系列本：ノート6

本読み

以下の本を読みます。

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

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

144ページからは、VAR モデル？ …これは、複数の時系列の過去の値で複数の時系列を回帰するモデルで、複数時系列を用いることによってよりよい予測ができるようにしたり、時系列化間の相互作用を調べる目的で利用するのか…具体的に、どんな時系列に対して適用するんだろう？

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

例えば春名さんの毎月の収入と、毎月のドーナツ消費量という2つの時系列データがあったとして、春名さんは収入が増えるほどドーナツを消費するかもしれませんが、好物のドーナツを消費することによって仕事に精が出て収入も増加するかもしれません。と考えると、春名さんの収入とドーナツ消費量の間には相互作用があるのではないかと疑われます。まあこの例は本に載っていないのでGranger因果性が存在するかはわかりませんが。

本に載ってたらおかしいだろ！ハルナのドーナツ消費行動が！

また沖本氏の本からの引用ですが、こちらの本の82ページに載っている例は日本の株式収益率、英国の株式収益率、米国の株式収益率でGranger因果性の存在が認められていますね。
経済・ファイナンスデータの計量時系列分析 (統計ライブラリー) | 沖本竜義 |本 | 通販 | Amazon

ジュン、ごめん、そのグレンジャー何とかの前に、こっちの本の145ページの、 SUR モデルっていうのがよくわからなかったんだけど。 $x_t$ と $y_t$ が互いに説明変数になっていなくて、ノイズどうしが相関をもってると SUR モデルっていうのか？でも $x_t$ と $y_t$ の説明変数のセットは共通なのに「見かけ上無相関な回帰モデル」と呼ぶのもよくわかんないし、そう呼んだとしてそれが何っていう…。

調べてみましょう。ウィキペディア（Seemingly unrelated regressions - Wikipedia）によると、SUR モデルというのは時系列の VAR モデルに限らず、一般の線形回帰モデルを複数束ねたものが、被説明変数が互いには説明変数に入っていない状態であって、ノイズが相関をもつかもしれない場合を指しているようですね。例えば以下のようなイメージです。

$\begin{cases} x = a_0 + a_1 x_1 + a_2 x_2 + a_3 x_3 + \varepsilon_x \\ y = b_0 + b_1 y_1 + b_2 y_2 + \varepsilon_y \end{cases}$

この2つの回帰式はぱっと見個別に推定できそうで、「無相関」に見えます。しかし、 $\varepsilon_x$ と $\varepsilon_y$ が無相関である保証はありません。なので、あくまで「見かけ上無相関」であり、それぞれの回帰式を個別に推定するのではなく2つの回帰式を同時にFGLS（実行可能一般化最小二乗法）で推定する必要があります（個別に最小2乗推定すると、パラメータの推定値の有効性が劣ります）。しかし、以下の2つの特殊なケースでは、個別に最小2乗推定してもFGLSと同じ結果が得られるんです。

$\varepsilon_x$ と $\varepsilon_y$ が無相関である場合（この場合は、見かけ上ではなく実際に無相関なので当然ですね）。
全ての回帰式で、右辺に登場する説明変数のセットが共通の場合。なんでこの場合に個別に最小2乗推定できるかは面倒なんで以下の論文の351ページまで見といてください。

Zellner, A. (1962). An Efficient Method of Estimating Seemingly Unrelated Regressions and Tests for Aggregation Bias. Journal of the American Statistical Association, 57(298), 348-368.

そして、ハヤトの言う通り、VARモデルはこの特殊なケースの後者なので、一般的に「見かけ上無相関な回帰モデル」というよりは、「説明変数が共通なのであまり見かけ上無相関にも見えないし、普通に無相関でないが、個別に最小2乗推定しても同時にFGLSしたときと同じ結果が得られる特殊なケース」と言った方がわかりやすいかもしれません。

なっが！

次にGranger因果性ですが、来期の春名さんのドーナツ消費量を、春名さんのドーナツ消費量のみを利用して予測した場合と、春名さんのドーナツ消費量に加えて春名さんの収入も利用して予測した場合で、後者の方が2乗誤差が小さくなるなら、春名さんの収入から春名さんのドーナツ消費量へのGranger因果性が存在するといいます。この検証にはF検定っぽいことをやります。

144ページに書いてあった「普段私たちが使っている因果という言葉の意味とは（後略）」ってのは？

それは146～147ページに説明がありますが、Granger因果性はその定義から「その変数があった方がよい予測ができる」でしかないです。Granger因果性の存在が確かめられても、通常の意味の因果関係、つまり、「春名さんの収入を変化させると、春名さんのドーナツ消費量も変化する」を抑えたわけにはならないんです。因果関係がなくてもGranger因果性の存在が確認されてしまうかもしれませんので。なので、複数時系列間の相互作用を別の切り口で調べる方法として、147ページにインパルス応答が挙げられていますね。

そういうことか。…161ページからは、ARCH・GARCH…また新しいのが出てきた…。

ARCHはノイズの分散を過去のノイズの大きさで回帰するようなモデルですね。ここまでずっとノイズ項の分散は $\sigma^2$ で固定でしたが、もはやこれも時間変化するということです。GARCHはノイズの分散を過去のノイズの大きさ及び過去のノイズの分散でも回帰します。過去のノイズの正負を考慮した拡張モデルも紹介されていますね。

その GJR モデルって実質、過去のノイズの大きさへの回帰係数を、過去のノイズの正負によって別の値に切り替えてるんだよな。それだと 0 で不連続になるし、結構適当な対応だなって思うんだけど、これでも強力ってことなのかな。…ともかくやっと本の前半終わった！！

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

2018-04-01

隼時系列本：ノート5

本読み

以下の本を読みます。

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

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

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

ちょっと調べてみましょうか。まず自己相関の有無の検定から。主に 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があれば）つづく

2018-03-13

隼時系列本：ノート4

本読み

以下の本を読みます。

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

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

時系列のモデルを決めるのに、次数は AIC か何かでゴリゴリ決めるけど、何階差分系列を取るべきかはそれじゃ駄目で、仮説検定をつかって対象データが単位根かどうか判断するって話で…KPSS検定とADF検定は帰無仮説と対立仮説が逆？

有名な単位根検定の中でKPSS検定だけが「単位根がない」の側を帰無仮説とするようですね。なるほど「単位根がある」の方を帰無仮説とするのが自然なように思います。両側検定では $H_0 : \mu = \mu_0 \, \; {\rm vs.} \, \; H_1: \mu \neq \mu_0$ の形式で $H_0$ が正しいかどうかを確かめますが、「単位根がある」とはつまり「その時系列を生成する確率過程を AR(1) で表現したとき $\phi_1 = 1$ である」ですので、 $H_0: \phi_1 = 1$ になりそうです。しかし、「単位根がある」を帰無仮説とするDF検定やADF検定では、帰無仮説を棄却できないケースに困ったようなんです。単位根がないというよほど「強い証拠」がない限り帰無仮説が受容されてしまうと（以下；KPSS検定の原論文のイントロダクションより）。

Testing the null hypothesis of stationarity against the alternative of a unit root

「単位根があることを棄却できない」ってことは「単位根がないと結論づけるには至らない」ってことだよな…それは困るな。棄却されたら単位根過程じゃないって安心できるけど、棄却されなかったらどっちなんだよっていう。

なので、KPSS検定の提唱者らは、帰無仮説と対立仮説を逆にしたんです。彼らは最初からランダムォーク項を仮定し、「ランダムウォークの分散がゼロである」を帰無仮説としました…本の64ページにある通りです。

なるほど、それなら帰無仮説が棄却されれば「単位根がある」ことになる…別の検定では単位根があるのかないのかどっちかわからなかったケースでも、きっぱりと「単位根がある」っていえるようになるかも。…じゃあKPSS検定だけでよくない？何階差分をとるか決めるときには「単位根がある」ってはっきり言ってくれる検定の方がうれしいし…。

実際この本ではKPSS検定を時系列モデル構築手順に据えていますが、65ページにデータに合わせて検定方法を変えるようにとありますね。どの検定を採用すべきかはどのような仮定を置いているかに依存しますし、時系列データの長さや、 $\phi_1$ の大きさにも依存するとも（ $\phi_1 = 0.1$ ならなかなか単位根があるようにはみえないかもしれませんが $\phi_1 = 0.9$ なら単位根があるようにみえてしまいそうですよね）。以下のような論文をWeb上で見つけました。

https://www.czso.cz/documents/10180/32912822/32019716q3047.pdf/09710b90-e1d0-4bb1-816e-5b83faad686b?version=1.0

結局どの検定にするかもケースバイケースってことか…そういえば、65ページのDF検定に「簡単そうに見えますが（略）うまくいきません」ってあったけど、まず簡単そうに見えるって気持ちがわかんない！

DF検定って式の形だけみると単なる単回帰分析にみえます。単回帰分析とは… $(x, y)$ のデータ対がたくさんあったとします。(気温, アイスクリームの売上高) のデータと考えてもいいです。このデータから傾き「1度気温が上がると、アイスの売上高がどれだけ伸びるか」を求めたいとします。単純な仮定の置き方の1つは、 $y$ が確率分布 $N(ax+b, \, \sigma^2)$ から生成されると決めることです。最小2乗法で $a, \, b$ の推定値 $a', b'$ を求めることができるでしょう。全てのデータが $xy$ 平面で $y$ 軸に平行に並んでいるのでもない限り、 $a'$ の値は必ず何かしら求まります。しかし、 $a'$ が求まるからといって先の仮定が正しいとは限りませんよね。本当は直線に当てはめるべきデータではないかもしれません。なので、 $H_0: a=0$ の下でその $a'$ が得られるかどうか確かめるんです。もし $H_0$ が棄却されれば、傾きが $a'$ であることが信頼できるでしょうし、棄却されなければ最初の仮定が信頼しづらいということです。…とまあここまでは簡単な話ですね。ただ、DF検定は単回帰分析のようにみえてこのように簡単にはいかないということなんです。

ごめんその話も別に簡単じゃないんだけど…何がさらに難しいの…。

$y_t= a x_t + \varepsilon_t$ と $y_{t} = \phi_1 y_{t-1} + \varepsilon_t$ は違うんです。後者では、 $\varepsilon_{t=0}$ の影響がいつまで経ってもなくなりません。ので、 $a$ の推定値 $a'$ の分布と、 $\phi_1$ の推定値 $\phi_1'$ の分布の形は違ってきます。さらに、いま帰無仮説は $H_0: \phi_1 = 1$ ですから、この非定常な場合での分布の形を知らないといけないんです。

…よくわかんないけど、単位根過程かどうか調べるのに苦労してるんだな。しかもそれだけ苦労したのにADF検定は帰無仮説が棄却できなかったりするんだろ？

なので、改良されたADF-GLS検定というのもあるそうです。と以下の資料に。

https://repository.kulib.kyoto-u.ac.jp/dspace/bitstream/2433/155556/1/kronso_183_2_35.pdf

へー。まあそれでARIMAモデルが決まったとして、残差をチェックするのか。…残差に自己相関が残ってたら駄目っていうけどさ、計算機が決めてくれたのに ARMA モデルの次数が最適になっていないってあるの？それとも残差が残っちゃったらもう ARMA モデルの限界ってこと？

僕もよくわかりませんが…ARMAのパラメータ最適化も不安定な場合もあるのかもしれませんが…あるいは、階差を取るべきなのに取っていなかったなどということもあるかもしれませんね。

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

2018-03-08

隼時系列本：ノート3

本読み

以下の本を読みます。

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

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

第2部の4章までで色々なモデルを扱って、第2部の5章からは、じゃあどのモデルを選べばいいかって話か…。

ここでモデルという言葉は ARIMA(2,0,0) と ARIMA(1,1,0) は異なるモデルということですね。なので、モデル選択というのは次数の決定のことであると書いてありますね。

あれ、AR モデルか MA モデルか ARIMA モデルか SARIMA モデルか、っていう選択を最初にするのかと思ったんだけど、その選択はないの？

まず、ARIMA モデルは AR モデルと MA モデルを包含しますよ。SARIMA モデルはさらに ARIMA モデルを包含しますが、周期を与えてつかうモデルなので、対象データに周期がある場合には SARIMA をつかってみればいいんじゃないですか。

あ、そうか。じゃあやっぱり最初に取りかかるのは次数の選択なのか。確かにARIMAモデルをつかうとしても次数をどうすればいいのか全然わからないよな…これはどういう手順で決めるのか気になるぞ…えっ、「手当たり次第」？ショックなんだけど…。

ショックを受けることないでしょう。僕たちは計算機が発達した時代に生まれたんです。考えなくても解決できる手順を中途半端に考える必要はありません。下手の考え休むに似たりです。しかし、計算機にまかせるにせよ、モデルのよさの指標については把握しておく必要があるでしょう。

う、うん…まずAICっていう指標が紹介されてるよ…それでその前に、尤度？尤度って何？この61ページのコイン投げの 2/9 っていう尤度は大きいの？小さいの？

コインの表が出る確率が 1/3 ではなく 1/2 だったらその尤度はどうなります？

それは、1/2×1/2＝1/4 になる。さっきの 2/9 と比べると 1/4 の方が大きい？

なので、あるコインについて「2回投げたら1回目が表で2回目が裏が出た」ということがわかっているとき、このコインの表が出る確率は 1/3 よりも 1/2 の方が尤もらしいということです。もっと一般に、コインの表が出る確率を $p$ とすると、「2回投げたら1回目が表で2回目が裏が出た」という結果に対する $p$ の尤度は $f(p)＝p(1-p)$ ですから、これを最大にするのは $p＝1/2$ です。なので結局、2回投げた結果しかわからない範囲では、このコインの表が出る確率として最も尤もらしいのは 1/2 です。

…2回投げて1回目が表で2回目が裏だったんだから、コインの表が出る確率が 1/2 らしいって当たり前なんじゃ？

まあ平均値については、「標本平均の期待値」が「真の平均値」に一致しますからね。でも平均値ではない分布のパラメータについては決して当たり前では…この話は脇道にそれるので置いておきましょう。…ところでハヤト、このコイン投げの結果から考えられる、もっと尤もらしいコインはあるでしょうか。

もっと尤もらしい？「2回投げたら1回目が表で2回目が裏が出た」という結果から考えられる？さっき尤度 $f(p)$ を最大にするのは $p=1/2$ だってジュン自分で言ってたじゃん。だから「表が出る確率が 1/2 のコイン」より尤もらしいコインはないだろ？

$f(p)$ はコインの表が出る確率が1回目も2回目も同じ場合の尤度ですので。2回目に投げるときは1回目に投げるときとはコインの性質が変化しているかもしれません。コインを特徴付けるパラメータを $p_1$ と $p_2$ の2個に増やして、1回目に投げるときは確率1で表が出て、2回目に投げるときは確率1で裏が出るコインを考えるとどうでしょう。こうすれば尤度は 1×1＝1 です。テキストの61ページと同様の例ですが。

えーそれはズルイだろ。

ええ、ずるいんです。ずるいというよりは、このモデリングは3回目に投げた結果を予測するのに全く有用ではないでしょう。なので、なるべくパラメータを増やさずに尤度を最大化したいんです。なるべくパラメータを増やさずに尤度を大きくできたかどうかの指標が62ページの AIC です。

62ページ…これがAICの式か…最大化対数尤度？

61ページに説明がありましたが、あるモデル、例えば ARMA(1,1) でも何でもいいです、ARMA(1,1) だと $c, \, \phi_1, \, \theta_1$ の3つのパラメータが推定対象ですね。対象時系列データの下での尤度が最大になるようにこれらのパラメータを最適化した上で、その尤度に対数をとったものが最大化対数尤度です。最後に対数をとる必要もありませんので、対数尤度が最大になるように3つのパラメータを最適化したときに対数尤度、と言った方がいいでしょうか。次数の異なるモデル間でよさを比較するのだから、パラメータセットはベストにしておかなければならないでしょう。

…時系列データに対するARIMAモデルの尤度ってどう計算するの？時系列データって、時点ごとにデータたくさんあるよな？

さっきのコイン投げの例でも、あるパラメータの下での2回のコイン投げの結果の同時確率がそのパラメータの尤度だったでしょう。2回だろうと100回だろうと一緒ですよ。与えられた100日分の気温に対するあるパラメータセットの尤度は $p($ 1日目～100日目の気温 $|$ パラメータセット $)$ です。もっとも、ARIMAモデルの尤度関数は…インターネットを探せばありますよ。例えばこのリンクなどがそうですね。

http://econ.nsysu.edu.tw/ezfiles/124/1124/img/Chapter17_MaximumLikelihoodEstimation.pdf

…まあいいや。その最大化対数尤度と、パラメータ数で AIC が決まるのか。尤度は大きいほどよくて、パラメータ数は小さいほどよくて、AICは小さいほどよい指標なんだな。…ジュンさっき、なるべくパラメータを増やさずに尤度を大きくしたいって言ってたよな。(2-41) 式をみるとさ、パラメータを1つ増やすんだったら、対数尤度は1よりも大きく増えてくれなきゃ困るってことだよな。なんでパラメータ1個と対数尤度1単位が釣り合うんだろう。あとなんで第1項も第2項も2がかかってるの？

以下のリンクに AIC の導出がありますよ。もっとも、AIC は原論文も探せばありますが。自由度 k のカイ2乗分布の期待値が k なので AIC は (2-41) 式のようになっているんですね。全体に2がかかっているのは、こうすると第1項が漸近的にカイ2乗分布にしたがうからですね。

http://www.is.titech.ac.jp/~shimo/pub/koza2000/ohp-20000703b.pdf

ふーん…まあなんやかんやで計算機が ARMA の次数を決めてくれるんだな。でも階差についてはそういうわけにいかないのか…それが63ページからの検定…帰無仮説って何？

仮説検定ですか。例えば有意水準5%で仮説検定するのだったら、「帰無仮説が正しいとしたら起きる確率が5%以下になるくらい珍しいこと」を「偶然それが起きたと考えるには低い」とみなすんです。与えられたデータが、「帰無仮説が正しいとしたらこんなデータが得られる確率は5%以下だ」というデータだったら、帰無仮説を棄却します。つまり、帰無仮説は正しくないと考えるんです。もし帰無仮説が正しいときに与えられたデータが得られる確率が5%よりも大きいのだったら、帰無仮説を受容します。これは、帰無仮説は正しくないとまではいえない、という状態ですね。

うーん…なんか具体例ない？

今調べたんですけど、高校2年生男子の平均体重は60.4kg、標準偏差は9.93kgらしいです。この平均体重が未知だったとして、僕とハヤトが高校2年生男子の無作為標本だったとき、帰無仮説「高校2年生男子の平均体重は60kg以上」は棄却されるでしょうか。有意水準5%とします。…分散既知の場合の平均の検定では、標本平均が正規分布にしたがうことを利用します。今回は標本平均が下側5%点より小さい領域が棄却域です。正規分布の下側5%点は -1.64 で、サンプル数2なので、棄却域は $60 - 1.64 \times \sqrt{9.93^2 / 2} = 48.48$ 未満ですね。僕のハヤトの体重の平均は 48.5kg なので棄却域になく、帰無仮説は棄却されません。帰無仮説が正しい下で5%も起こりえない珍しいデータではないということですね。

でもすごいギリギリ…。

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