回帰係数の共分散がわからず困っています（確率・統計＞回帰分析）

現在、卒論執筆に際し、以下のような統計の問題にぶつかり、大変困っています。

【問題】
共通の説明変数x1～xnを持つ、N本の回帰モデルを以下の通り考える。

y1 = a1×e + b11×x1 + b12×x2 + … + b1m×xm + ε1
y2 = a2×e + b21×x1 + b22×x2 + … + b2m×xm + ε2
…
yn = an×e + bn1×x1 + bn2×x2 + … + bnm×xm + εn
（x,y,εはベクトル、eは単位ベクトル）

このとき、各回帰式の切片a1～anの分散共分散行列V(a)を求めたいです。

分散共分散行列の対角成分に関しては、Excelの分析ツールで出力される標準誤差の項を基に計算できるのですが、共分散項については、自分で式の展開を試みるもどうにも上手くいきません。

また回帰分析の本は手当たり次第に読みましたが、一本の回帰式の中の係数の共分散行列ならば書かれていましたが、異なる回帰式の共分散行列となると、どの本にも書かれておらず、まったく先の見えない状況です。

そのものずばりの回答でなくても、そのことについて書かれている書籍や論文等がございましたら、お知らせいただければ幸いです。

また、上記の一般化された問題をもっと簡略化した以下の問題でもわかる方がいらっしゃいましたら、どうかお願いします。

共通の説明変数xを持つ2つの回帰モデルを考える。

y = ax + b
z = cx + d

このとき、切片b,dの共分散cov(b,d)を求めよ。

No.2 回答者： quaestio 回答日時： 2009/01/09 21:09

説明し忘れた部分（下記の文中の括弧で括った部分）がありました。

（前略）
ここで、
yijがi行j列目の要素であるn行m列の行列Y
βijがi行j列目の要素であるm行n列の行列β
xijがi行j列目の要素であるk行m列の行列x（ただし、xi1=1）
を定義すると上の式は
（後略）

No.1 回答者： quaestio 回答日時： 2009/01/08 22:54

卒論執筆のためとのことであまり回答しにくい質問ですね。

回答があったとしてももしその回答が間違っていて卒論がうまくかけなかったとしても責任が取れませんので、担当の先生に聞くのが一番だと思います。

とはいえ全く回答がついていないのもあんまりだと思いますので、回答してみましょう。

行列を使って説明しますので、表記を少し変えます。
回帰モデルは

yij = β1j + β2j * xi2 + β3j * xi3 + … + βmj * xim + εij
E(εij * εpq) = σjq (i=p)
E(εij * εpq) = 0 (i≠p)
（表記を少し変更しています。E()は期待値記号です。）

ということでいいでしょうか？　つまり、説明変数がx12～x1mのときy11～y1mという応答が、…説明変数がxk2～xkmのときyk1～ykmという応答が得られたということです。

ここで、
yijがi行j列目の要素であるn行m列の行列Y
βijがi行j列目の要素であるm行n列の行列β
xijがi行j列目の要素であるk行m列の行列x
を定義すると上の式は

Y = Xβ+ε

と簡潔に表記できます。
これの回帰式を求めると係数行列βの推定値βhatは

βhat = (X'X)^(-1)X'Y

となります。ここで、X'は行列Xの転置行列、(X)^(-1)は行列Xの逆行列を意味します。
1行m列の行列J = (1, 0,…, 0)を定義すると、求めたい切片の分散共分散行列Vは

V = E[{Jβhat-E(Jβhat)}'{Jβhat-E(Jβhat)}]

となります。

E(Jβhat) = E{J(X'X)^(-1)X'Y}
= E{J(X'X)^(-1)X'(Xβ+ε)}
= E{J(X'X)^(-1)X'(Xβ+ε)}
= E{J(X'X)^(-1)(X'Xβ+X'ε)}
= E{J(β+(X'X)^(-1)X'ε)}
= Jβ+E{J(X'X)^(-1)X'ε)}
= Jβ

ですので

V = E[{J(X'X)^(-1)X'ε}'{J(X'X)^(-1)X'ε}]
= E[ε'X(X'X)^(-1)J'J(X'X)^(-1)X'ε]

ここで、

α = X(X'X)^(-1)J'J(X'X)^(-1)X'

とおいて、分散共分散行列Vのi行j列の要素をVij、行列αのi行j列の要素をαijとすると、

Vij = E{ΣΣ(εpi * αpq * εqj)} (p, qについての和）

であり、

E(εij * εpq) = σjq (i=p)
E(εij * εpq) = 0 (i≠p)

であることに注意すれば、

Vij = σij * Σαpp (pについての和）
= σij * tr(α)

となる。tr(α)は行列αのトレースを意味する。
行列のトレースの性質より、

Vij = σij * tr(α)
= σij * tr(X(X'X)^(-1)J'J(X'X)^(-1)X')
= σij * tr(J(X'X)^(-1)X'X(X'X)^(-1)J')
= σij * tr(J(X'X)^(-1)J')

となる。したがって、切片β1iとβ1jの分散或いは共分散は、(X'X)^(-1)の1行目の1列目の要素に、i番目の回帰式の誤差項とj番目の回帰式の誤差項の分散或いは共分散をかけたものに等しいことがわかります。

内容に間違いないかはご自分でもご確認ください。

この回答へのお礼

お礼が遅れまして申し訳ございません。

大変丁寧にご解答頂き誠に感謝しております！

自分でも確認しながら再度解いてみます。

本当にありがとうございました。

お礼日時：2009/01/14 02:30