共分散行列$Q$は2x2で、右下に$\sigma_e^2$がはいり、残りが0。こいつをカルマンフィルタで推定する。Rコードはこんな感じ。遷移行列$T$に未知パラメータが入っているので、ちょっと面倒くさい。Mplusでベイズ推定。実はN=1の時系列分析についても便利な機能を持っているのである。Mplusのコードはこんな感じ。y@0;MODELコマンドがわかりにくいが、1行目は「vは潜在変数ですが指標を持っていません。モデルのなかでラグ1を使わせて下さい」。2行目は「その残差分散(つまり$\sigma_e$)を自由推定したいです」。3行目で$v_t$の自己回帰を定義し(v&1とは$v_{t-1}$を表す)、4行目で$y_t$のモデルを定義している($v_t$の回帰係数は1に固定している)。ほっとくと$y_t$の切片と残差分散を推定してしまうので、最後の2行でそれを抑止している。Stanでベイズ推定。Stanファイルはこんな感じ。$\sigma_e$より大きくなるからである。(1-\rho^2)$である。このモデルとともに、parametersブロックでrhoの範囲制約をなくし、かつ下から3行目を消したモデルも推定してみた。まとめると、出場選手は以下の10名である。ベイズ推定(Stan)推定の様子この10人の選手に、11水準の$\rho$の各1000個のデータセットについて、$\beta$を推定させた。Stanは時間がかかるので、1000個のうち100個だけについて推定するだけで勘弁してやった。1$と仮定)による推定結果を図にしてみた。縦軸は$\rho$の推定値である。真値が赤で表現してある。してみると、真の$\rho$が1に近づくにつれ、$\rho$の推定値は0に向かって歪むみたいですね。横軸と縦軸の間にはあまり相関がなさそうだ。$\beta$の推定値$\hat{\beta}$とのずれである。その大きさを次の2つの指標で評価しよう。$\hat{\beta}$の平均から1を引いた値。RMSE、すなわち$(\hat{\beta}-1)^2$の平均の平方根。0に近いほうがありがたい。F5のみ試行回数が100である点に注意。シンプルなチャートだが、なかなか情報量が多いので、順にみていこう。0$と勝手に想定し、OLS回帰をやってしまった場合。真の$\rho$がそれほど大きくなければバイアスはそれほど大きくない。RMSEは$\rho$とともに増大する。$\rho$が1に近づくと、バイアス・RMSEともに急増する。撹乱項が非定常に近づいているわけだから、パラメータがうまく推定できないのも道理である。あれですね、時系列を分析する際には、きちんと自己相関を考慮することが大事ですね。次に、$|\rho|<1$と決め打ちした場合(図のパネル2)と、$\rho=1$と決め打ちした場合(パネル3)について。バイアスはそれほど大きくないようだ。RMSEのみ拡大して示す。選手F2はあとで観察することにして、ここでは省略する。選手B3とA3はぴったり重なってしまっている。$|\rho|<1$と決め打ちする場合、最尤推定するかFGLS推定するかにはたいしたちがいがないわけだ。そりゃそうか、きっと推定値は同じだ。$\rho$が小さいとき、$\rho=1$と決め打ちするとRMSEが大きくなる。誤ったモデルを指定したわけだから、これは当然である。これをみて大変に意外だった点がふたつある。第1に、$\rho$が大きくなるにつれ、つまり撹乱項の自己相関が大きくなるにつれ、$\beta$のRMSEが小さくなっているという点である。???もうひとつ、ふへええと言葉にならないため息をついたのは、$\rho$が1に近づくほど、$|\rho|<1$と仮定して定常AR(1)撹乱項を持つ回帰モデルを推定しても、$\rho=1$と仮定して差分時系列について回帰モデルを推定しても、違いがなくなってしまう点である。さきにみたように、単位根検定がしくじりやすいのは$\rho$が1に近いときである。いうときほど、実は間違えたところでたいした実害はないわけだ。このように、単位根検定でモデルを切り替えようが、頭から$|\rho|<1$と信じこんで定常AR(1)撹乱項を持つ回帰モデルを使い続けようが、ほとんどちがいがないようだ。$\rho=1$だったとしても、である。単位根検定をやったのに、意味なかったわけ...?最後に、$\rho$の範囲について制約しないモデル。元の図にもどると、選手E5がとんでもないバイアスを持っていることがわかる。冗談はともかく、ひょっとすると私がコードを間違えているのかもしれない。残念ながらMplusくんは休場とし、他の選手についてRMSEを拡大してみよう。