餡子付゛録゛

ソフトウェア開発ツールの便利な使い方を紹介。

完全情報最尤推定法(FIML)で操作変数のある連立方程式モデルを推定

計量経済学では同時方程式を考えることが多く、そのために一般最小二乗法(OLS)の発展である操作変数法(IVM)とそのバリエーションである二段階最小二乗法(2SLS)*1や三段階最小二乗法(3SLS)、モーメント法の発展である一般化モーメント法(GMM)がよく使われているのですが、これらは尤度函数を用いないのでベイズ推定にできません。頻度主義者がベイジアンになるのは難しいです。
しかし、完全情報最尤推定法(FIML)と言うのがあり、実際のところは簡単に最尤法で同時方程式を推定できます。FIMLはRでは欠損値処理アルゴリズムのように扱われている気がするのですが、教科書的な手法なので確認しておきましょう*2。最尤法だけに不均一分散や欠損値処理の対処も工夫できますし、ベイズ推定しない場合でも、便利に使える局面は多そうです。

理論モデルとOLS推定量に入る同時性バイアス

推定の例をつくるわけですが、需要曲線と供給曲線を考えましょう*3


需要関数: d = \alpha_{10} + \alpha_{11} p + \xi_1
供給関数: s = \alpha_{20} + \alpha_{21} p  + \alpha_{22} z + \xi_2
均衡条件: d=s
 dは需要、 sは供給、 pは価格、 zは天候などの操作変数、 \alphaは係数、 \xiは誤差項です。ここで均衡価格を計算すると、

均衡価格: p^*=\frac{\alpha_{20} - \alpha_{10}}{\alpha_{11} - \alpha_{21}} + \frac{ \alpha_{22}}{\alpha_{11} - \alpha_{21}} z  + \frac{ \xi_2 - \xi_1 }{\alpha_{11} - \alpha_{21}}
となります。観測データは均衡価格 p^*になるわけですが、例えば需要関数の推定で p^*を説明変数に用いると、 p^*の第3項に \xi_1があるため、 COV(p^*, \xi_1) \neq 0となりOLSの前提を満たしません。需要関数の価格の係数の推定量 \hat{\alpha_{01}}_{OLS}=\alpha_{01}+COV(p^*, \xi_1) とバイアスが入ります。

完全情報最尤推定

需要関数と均衡価格式を同時推定することを考えましょう。推定する係数を \beta、誤差項を \muと置きなおし、方程式一本ごとに誤差項があると考え、二本の方程式を行列にまとめて書くと、次のようになります。 nは観測数です。

 \begin{pmatrix}
d_1 & p_1 \\
d_2 & p_2 \\
\vdots & \vdots \\
d_n & p_n \end{pmatrix}  = \begin{pmatrix}
1 & z_1 & p_1\\
1 & z_2 & p_2\\
\vdots & \vdots & \vdots \\
1 & z_n & p_n
\end{pmatrix}  \begin{pmatrix}
\beta_{11}  & \beta_{12}  \\
0  & \beta_{22} \\
\beta_{31}  & 0 \end{pmatrix} +  \begin{pmatrix}
\mu_{11}  & \mu_{12}\\
\mu_{21}  & \mu_{22}\\
\vdots & \vdots \\
\mu_{n1} & \mu_{n2} \end{pmatrix}

前節の係数とは、 \beta_{11}=\alpha_{10} \beta_{12}=(\alpha_{20} - \alpha_{10})/(\alpha_{11} - \alpha_{21})と言った対応関係があります。
右辺の内生変数と外生変数を分離します。

 \begin{pmatrix}
d_1 & p_1 \\
d_2 & p_2 \\
\vdots & \vdots \\
d_n & p_n \end{pmatrix}  = \begin{pmatrix}
1 & z_1 \\
1 & z_2 \\
\vdots & \vdots \\
1 & z_n
\end{pmatrix}  \begin{pmatrix}
\beta_{11}  & \beta_{12}  \\
0  & \beta_{22} \end{pmatrix} + \begin{pmatrix}
p_1\\
p_2\\
\vdots \\
p_n
\end{pmatrix}  \begin{pmatrix}
\beta_{31}  & 0 \end{pmatrix} +  \begin{pmatrix}
\mu_{11}  & \mu_{12}\\
\mu_{21}  & \mu_{22}\\
\vdots & \vdots \\
\mu_{n1} & \mu_{n2} \end{pmatrix}

内生変数を左辺にまとめます。

 \begin{pmatrix}
d_1 & p_1 \\
d_2 & p_2 \\
\vdots & \vdots \\
d_n & p_n \end{pmatrix} - \begin{pmatrix}
p_1\\
p_2\\
\vdots \\
p_n
\end{pmatrix}  \begin{pmatrix}
\beta_{31}  & 0 \end{pmatrix} = \begin{pmatrix}
1 & z_1 \\
1 & z_2 \\
\vdots & \vdots \\
1 & z_n
\end{pmatrix}  \begin{pmatrix}
\beta_{11}  & \beta_{12}  \\
0  & \beta_{22} \end{pmatrix} +  \begin{pmatrix}
\mu_{11}  & \mu_{12}\\
\mu_{21}  & \mu_{22}\\
\vdots & \vdots \\
\mu_{n1} & \mu_{n2} \end{pmatrix}

左辺第1項と第2項をまとめましょう。

 \begin{pmatrix}
d_1 & p_1 \\
d_2 & p_2 \\
\vdots & \vdots \\
d_n & p_n \end{pmatrix} \begin{pmatrix}
1 & 0 \\
 - \beta_{31}   & 1 \end{pmatrix} = \begin{pmatrix}
1 & z_1 \\
1 & z_2 \\
\vdots & \vdots \\
1 & z_n
\end{pmatrix}  \begin{pmatrix}
\beta_{11}  & \beta_{12}  \\
0  & \beta_{22} \end{pmatrix} +  \begin{pmatrix}
\mu_{11}  & \mu_{12}\\
\mu_{21}  & \mu_{22}\\
\vdots & \vdots \\
\mu_{n1} & \mu_{n2} \end{pmatrix}

この行列をいちいち書いていられないので、書き直します。

 Y\Gamma = W B + U

FIMLでは、この構造方程式を推定します。こう書くとOLSで入るバイアスが入らない理由が分かりづらいですが、この構造方程式の尤度函数は説明変数に内生変数が入らない誘導形 Y = W B \Gamma^{-1} + U \Gamma^{-1} = W\Pi + Vの尤度函数から導かれます*4。また、数理的な特性を使って整理すると以下のようなconcentrated loglikelihood function*5にできることが知られています*6


 \Sigma = \frac{1}{n} (Y\Gamma - WB)^{t} (Y\Gamma - WB)
 \log L(\beta) = \frac{-gn}{2}(\log 2\pi + 1)  + n \log | \Gamma | - \frac{n}{2} | \Sigma |
 gは内生変数(もしくは左辺にまとめられた被説明変)の数で、ここでは2です。 \beta = ( \beta_{11},  \beta_{12},  \beta_{22}, -\beta_{31} ) の4つのパラメーターを上の対数尤度函数で推定すればよいです。

データ

理論モデルに沿った式で、乱数を使ってデータを生成します。

set.seed(1231)
n <- 300
a10 <- 100
a11 <- -1
a20 <- 2
a21 <- 3
a22 <- 4
z <- runif(n, min=0, max=3)
mu <- rnorm(n, mean=0, sd=2)
nu <- rnorm(n, mean=0, sd=1)
p <- (a20 - a10 + a22*z + nu - mu)/(a11 - a21)
d <- a10 + a11*p + mu
s <- a20 + a21*p + a22*z + nu

推定

concentrated loglikelihood functionをoptim関数で推定します。

Y <- cbind(d, p)
W <- cbind(rep(1, n), z)
g <- 2

llf <- function(p){
    B <- matrix(c(p[1], 0, p[2], p[3]), 2, 2)
    G <- matrix(c(1, p[4], 0, 1), 2, 2) # Γ
    E <- Y%*%G - W%*%B
    S <- t(E)%*%E/n # ∑
    -g*n/2*(log(2*pi) + 1) + n*log(det(G)) - n/2*log(det(S))
}

r_optim <- optim(c(1, 0, 0, 0), function(p){ -llf(p) }, method="BFGS")

r_optim$parが推定された \hat{\beta}になりますが、内生変数の係数となる4番目のパラメーターの符号が逆になることに注意してください。なお、optim関数の引数にhessian=TRUEをつけるとヘッセ行列がとれるので、標準誤差を計算することができます*7

OLSとIVとの比較

推定結果の比較をしておきましょう。FIMLはOLSよりはかなりマシ、IVよりは僅かに劣るといった感じでしょうか。

係数 真の値 OLS IV FIML
 a_{10} 100 75.95501 100.11595 100.1228135
 a_{11} -1 0.03949 -1.01101 -1.0113085

なお、最尤法は漸近的に有効だけに、サンプルサイズが小さいとIVとの差が広がる傾向があります。IV推定量は以下のように計算できるので、nを変えて試してみてください。

zm <- matrix(c(rep.int(1, n), c(z)), n, 2)
xm <- matrix(c(rep.int(1, n), c(p)), n, 2)
solve(t(zm) %*% xm) %*% (t(zm) %*% d)

*1:操作変数が1つの場合はIVMと同値なのですが、一応、別なので。

*2:昔読んだGreene (2003)やWooldridge (2002)では最尤法で連立方程式モデルの推定の仕方は書いていなかったのですが、Davidson and MacKinnon (2003)では第12章で詳しく説明されています。

*3:Hayashi (2000)にある例を踏襲しています。

*4:Davidson and MacKinnon (2003)もしくはオンラインで配布している2021年版のpp.503–505とpp.544–545を参照。誘導形はSURと同じ多変量正規分布の尤度函数になるので、 \Gammaによる変形の影響を整理すれば構造形の尤度函数が導出できます。

*5:定訳不明

*6:一階条件の節の一階条件に直接関係のないExercise 12.24を \Sigma \Sigma^{-1}=Iに注意して \Sigma^{-1}の成分を考えて解くことで得ることができる・・・ハズ。

*7:Rと手作業で覚える最尤法 - 餡子付゛録゛