経済学では複数の条件式を同時に満たす均衡点を考えるため、計量経済学では同時方程式を取り扱うことが多いのですが、含まれる内生変数は一般最小二乗法（OLS）では説明変数として扱えません。そのため、操作変数法（IVM）を扱うことが多いです。

1 同時方程式

内生変数が出てくる状況は多様ですが、需給の分析が典型的です。供給関数と需要関数を定義してみましょう。

\[\begin{align} \mbox{需要関数: }d &= \alpha_{10} + \alpha_{11} p + \xi_1 \\ \mbox{供給関数: }s &= \alpha_{20} + \alpha_{21} p + \alpha_{22} z + \xi_2 \end{align}\]

\(d\)は需要、\(s\)は供給、\(p\)は価格、\(z\)は天候などの外生変数、\(\alpha\)は係数、\(\xi\)は誤差項です。需給が一致したときに取引は行われるはずなので、観測データは

\[ \mbox{均衡条件: }d = s \]

を満たします。

1.1 価格の内生性と同時性バイアス

この均衡条件を満たす均衡価格\(p^*\)を計算すると、

\[ \mbox{均衡価格: }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\)があるため、\(\mathrm{COV}(p^*, \xi_1) \neq 0\)となり一般最小二乗法（OLS）の前提を満たしません。需要関数の価格の係数の推定量は

\[ E[\hat{\alpha_{11}}_{OLS}]=\alpha_{11}+\frac{\mathrm{COV}(p^*, \xi_1)}{\mathrm{VAR}(p^*)} \]

とバイアスが入ります。

この価格のように連立方程式の中で定まる変数のことを内生変数と呼びます。

2 二段階最小二乗法（TSLS）

価格が需要に影響すると同時に、需要が価格に影響するために、価格と誤差項の共分散が非ゼロとなりました。OLSが一致推定量になる前提が満たされません。

この解決策として、2段階最小二乗法（2SLS/TSLS）と呼ばれる手法が広く知られています。具体的には

内生変数である価格と相関し、かつ従属変数である需要で説明されないデータ（操作変数）を探す
内生変数以外の説明変数（外生変数）と操作変数で、内生変数にOLSをかける
OLSの結果から内生変数の予測値をつくる
内生変数の予測値を、内生変数の代わりに使って、従属変数を回帰

します。内生変数の予測値は、従属変数で説明されない変数で構成されるので、内生性の問題は生じません。理屈上はTSLSで解決です。

ただし、実用しやすいかと言うとそうでもないです。

操作変数が見つかりません。需要と供給の分析では、よく天候データが使われていますが、農作物以外はどうしたらいいかわかりませんし、農作物でも厳密には生産地と消費地の天候が異なるものである必要があります¹。また、操作変数が内生変数と十分に相関しない（弱相関）とダメですし、複数の操作変数を使うと内生性が出てくる過剰識別の問題が生じます。今まで使われてこなかった有用な操作変数を見つけると、経済学の博士論文が書けると言われているぐらいです。
多重共線性の問題が生じやすいです。内生変数の予測値は、通常、すべての外生変数を用いて作成します。操作変数以外の外生変数を省いた場合、予測値の誤差が外生変数と相関していると、二段階目のOLSで外生変数の係数の推定値に同時性バイアスが入りこむことになるからです。一方、すべての外生変数を用いて予測値をつくると、予測値と外生変数の相関係数があがり、分散拡大係数（VIF）が上昇し、推定係数の標準誤差が跳ね上がっていきます。解決のためには、大きなサンプルサイズが必要です。

出版されている学術論文でも、TSLS（と後述するIVM）を使った分析は、苦しいものが多いです。

3 操作変数法（IVM）

TSLSの操作変数と内生変数の数は一致しなくても操作変数の方が大きければ推定可能ですが、一致している場合を丁度識別と呼びます。また、口頭でTSLSとIVMを厳密に言い分けているわけではないですが、丁度識別のときの推定方法をIVMと呼びます。まず、IVM推定量の導出を確認しましょう。

3.1 係数

上の需給の連立方程式モデルよりも、より一般的な形で考えます。

従属変数を\(y\)とおき、内生性のある説明変数を含む行列を\(X\)とします。ここで、\(y\)とは独立で、かつ内生変数と強く相関する操作変数を見つけてきます。\(X\)の内生変数の列を操作変数の列に置き換えた行列を\(Z\)とします。

\(Z\)を用いた\(X\)のOLSからの予測値は\(\hat{X}\)は、推定量に\(X\)を左から乗じたものになるので

\[ \hat{X} = Z({}^t\!ZZ)^{-1}{}^t\!ZX \]

となります。プロジェクション行列を\(P = Z({}^t\!ZZ)^{-1}{}^t\!Z\)と置いて

\[ \hat{X} = PX \]

\(\hat{X}\)で\(y\)を回帰して係数\(\beta\)を求めると、

\[\begin{align} \hat{\beta}_{IV} &= ({}^t\!\hat{X}\hat{X})^{-1}{}^t\!\hat{X}y \\ &= ({}^t\!(PX)PX)^{-1}{}^t\!(PX)y \\ &= ({}^t\!ZX)^{-1} {}^t\!Z y \end{align}\]

となります。

3.1.1 一致性の確認

線形モデルが正しければ観測値\(y\)は説明変数\(X\)とその係数\(\beta\)と誤差項\(\epsilon\)で、

\[ y = X\beta + \epsilon \]

と書けます。IVMによる推定量は

\[\begin{align} \hat{\beta}_{IV} &= ({}^t\!ZX)^{-1} {}^t\!Z y \\ &= ({}^t\!ZX)^{-1} {}^t\!Z (X\beta + \epsilon) \\ &= \beta + ({}^t\!ZX)^{-1} {}^t\!Z \epsilon \\ &= \beta + \bigg(\frac{1}{n}{}^t\!ZX\bigg)^{-1} \frac{1}{n}{}^t\!Z \epsilon \end{align}\]

と展開できます。ここで\(n\)はサンプルサイズです。

\(Z\)と\(X\)に相関があれば、

\[ \lim_{n \rightarrow \infty} \bigg(\frac{1}{n}{}^t\!ZX\bigg)^{-1} \ne 0 \]

です。\(Z\)と\(\epsilon\)に相関が無ければ、

\[ \lim_{n \rightarrow \infty} \frac{1}{n}{}^t\!Z \epsilon = 0 \]

です。よって、操作変数が誤差項と相関せず、内生変数と相関すると言う操作変数の条件を満たしていれば、

\[ \lim_{n \rightarrow \infty} \hat{\beta}_{IV} = \beta \]

となるため、IVMによる推定量は一致推定量になります。

3.2 分散共分散行列

OLSと同様に導出できます。

\[\begin{align} \text{VAR}[ \hat{\beta}_{IV} \mid X ] &= E[^t\!(\hat{\beta}_{IV} - \beta)(\hat{\beta}_{IV} - \beta) \mid X]\\ &= E[(^t\!ZX)^{-1}{^t\!Z} \epsilon {^t\!\epsilon} Z (^t\!ZX)^{-1} \mid X] \\ &= (^t\!ZX)^{-1}{^t\!Z} E[\epsilon {^t\!\epsilon} \mid X] Z (^t\!ZX)^{-1}\\ &= (^t\!ZX)^{-1}{^t\!Z} (\sigma^2 I) Z (^t\!ZX)^{-1}\\ &= \sigma^2 (^t\!ZX)^{-1}{^t\!Z} Z (^t\!ZX)^{-1} \end{align}\]

この値は\(\sigma^2 ({}^t\!\hat{X}\hat{X})^{-1}\)と一致します。ただし、\(\sigma^2\)は、\(\hat{X}\)では無く、\(X\)の条件付き分散であることに注意してください。つまり、\(n\)をサンプルサイズ、\(k\)を説明変数の数として、

\[ \hat{\sigma}^2 = \frac{(y - X \hat{\beta}_{IV})^2}{n - k} \]

となります。

パッケージを用いずTSLSで推定したときに、\(X\)の代わりに\(\hat{X}\)を入れてしまうと、分散を小さく推定してしまいます。\(\hat{X}\)と\(X\)の差になる誤差を無視することになるからです。

第一段階の推定の相関係数が高いか、第二段階の推定における内生変数の予測値の寄与が（操作変数以外の外生変数のおかげで）小さいか、もしくはその両方であれば、誤って\(\hat{X}\)を使っても\(\hat{\sigma}^2\)は何割か程度しか変わらない事になりますが、パッケージを用いて推定する方が無難です。

3.3 IVMの例

IVMを試してみましょう。上手く推定できることを確認するために、上の需給モデルに沿って乱数から生成したデータセットを用います。

3.3.1 データセット

需要関数は\(d = 100 - p\), 供給関数は\(s = 2 + 3p\)を仮定します。

set.seed(1231)
df01 <- with(new.env(), {
    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
    data.frame(s = s, d = d, p = p, z = z)
})

3.4 推定

行列計算で需要関数の係数と分散共分散行列を出します。

X <- model.matrix(~p, df01)
Z <- model.matrix(~z, df01)
(beta <- solve(t(Z) %*% X) %*% t(Z) %*% df01$d)

                 [,1]
(Intercept) 100.11595
p            -1.01101

s2 <- sum((df01$d - X %*% beta)^2) / (nrow(df01) - ncol(Z))
(vcov <- s2 * solve(t(Z) %*% X) %*% t(Z) %*% Z %*% t(solve(t(Z) %*% X)))

            (Intercept)           p
(Intercept)  10.8557453 -0.47138460
p            -0.4713846  0.02049553

概ね、DGPの係数通りになりました。

3.4.1 誤った分散共分散行列

上で説明した間違った分散共分散行列の計算をして、誤った係数の標準誤差を見てみましょう。

r_lm_1st <- lm(p ~ z, df01)
(r_lm_2nd <- lm(d ~ predict(r_lm_1st), df01))


Call:
lm(formula = d ~ predict(r_lm_1st), data = df01)

Coefficients:
      (Intercept)  predict(r_lm_1st)  
          100.116             -1.011

# 誤った係数の標準誤差
sqrt(diag(vcov(r_lm_2nd)))

      (Intercept) predict(r_lm_1st) 
        2.4767573         0.1076176

# 正しい係数の標準誤差
sqrt(diag(vcov))

(Intercept)           p 
  3.2948058   0.1431626

# 誤った場合と正しい場合の比
(sqrt(diag(vcov(r_lm_2nd)))/sqrt(diag(vcov)))[1]

(Intercept) 
  0.7517157

操作変数しかないわけですが、標準誤差は25%ぐらいの過小評価にとどまっています。

3.5 不均一分散

ウェイトをつけることで不均一分散に対処するようにTSLS/IVMを発展させることもできますが、不均一分散がどのように入るかの仮定によって手法が多様になるので本稿では割愛します。

4 TSLSの実践

だいぶ前に、semパッケージのtsls関数を使って、農林水産省の統計情報と気象庁の気象統計情報から、だいこん需要関数の推定を行ったときのRのコードを紹介します。

# semパッケージを読み出す
library(sem)
# テキスト・ファイルからデータを読み込む
df02 <- read.table("https://wh.anlyznews.com/R/dataset/for_tsls.txt", header=TRUE, sep="\t", fileEncoding="UTF-8")
# 年次ダミーを作成
month <- factor(df02$月)
# 年次ダミーを作成（i年10月から(i+1)年9月までをi年とする）
df02["年度"] <- with(df02, ifelse(月 >= 10, 年 + 1, 年))
year <- factor(df02$年度)
# 被説明変数と内生変数と外生変数の推定式
# 係数が需要の価格弾力性をそのまま表すように、変数は対数をとっている
frml <- log(出荷量) ~ log(価格) + month + year
# 操作変数と外生変数を指定
zm <- ~ log(降水量.0) + log(日照時間.3) + month + year
# 二段階最小二乗法で推定
r_tsls <- tsls(frml, zm, data = df02)
# 推定結果を表示
summary(r_tsls)


 2SLS Estimates

Model Formula: log(出荷量) ~ log(価格) + month + year

Instruments: ~log(降水量.0) + log(日照時間.3) + month + year

Residuals:
     Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
-0.136924 -0.029823 -0.002417  0.000000  0.028451  0.102018 

               Estimate  Std. Error  t value   Pr(>|t|)    
(Intercept) 12.54087534  0.44125411 28.42098 < 2.22e-16 ***
log(価格)   -0.33224585  0.10338916 -3.21355 0.00205453 ** 
month2       0.02824776  0.03002783  0.94072 0.35038603    
month3       0.10241160  0.03407292  3.00566 0.00378196 ** 
month4       0.04585596  0.04128979  1.11059 0.27090071    
month5      -0.02160695  0.03289950 -0.65676 0.51369235    
month6      -0.05684150  0.02872577 -1.97876 0.05214928 .  
month7      -0.03827124  0.03008154 -1.27225 0.20788747    
month8       0.05719311  0.03754859  1.52318 0.13264093    
month9       0.28282606  0.03871574  7.30520 5.4093e-10 ***
month10      0.39872012  0.03151535 12.65162 < 2.22e-16 ***
month11      0.13308131  0.02740148  4.85672 8.0400e-06 ***
month12      0.10968422  0.03154317  3.47727 0.00091584 ***
year2005    -0.03049104  0.02241152 -1.36051 0.17844172    
year2006     0.01249078  0.02284169  0.54684 0.58638949    
year2007    -0.13721037  0.03190916 -4.30003 5.9475e-05 ***
year2008    -0.07893549  0.02425310 -3.25466 0.00181582 ** 
year2009    -0.08253713  0.02208109 -3.73791 0.00039809 ***
year2010    -0.10265711  0.02228360 -4.60685 1.9996e-05 ***
year2011    -0.05214745  0.03801860 -1.37163 0.17496834    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.0494751 on 64 degrees of freedom

最近はsem::tslsではなく、AER::ivregを使うのが一般的のようです。

library(AER)
r_ivreg <- ivreg(log(出荷量) ~ log(価格) + month + year | log(降水量.0) + log(日照時間.3) + month + year, data = df02)
# 推定量と後述する検定に対応する検定の結果を表示する
summary(r_ivreg, diagnostics = TRUE)

4.1 行列による計算

操作変数ではない外生変数を一段階目の推定に入れるのを忘れたり、誤差項の分散の推定に予測値を使ってしまったり、ありがちな間違いがあるので推奨されませんが、行列で計算することもできます。

X <- model.matrix(frml, df02)
Z <- model.matrix(zm, df02)
P <- Z %*% solve(t(Z) %*% Z) %*% t(Z)
y <- log(df02$出荷量)
beta_tsls <- solve(t(X) %*% P %*% X) %*% t(X) %*% P %*% y
s2_tsls <- sum((y - X %*% beta_tsls)^2)/(nrow(X) - ncol(X))
vcov_tsls <- s2_tsls * solve(t(P %*% X) %*% (P %*% X))
se_tsls <- sqrt(diag(vcov_tsls))

なお、\(X\)と\(Z\)のランクが異なるため、IVMの導出途中のような式で計算しますが、手順は同様です。

4.2 弱相関検定

前提条件のひとつ、内生変数と操作変数がしっかり相関しているかテストします。具体的には、操作変数の変わりに定数\(1\)を用いて推定し、分散分析で有意に差があるかを確認します。

#  操作変数の変わりに定数 1 を用いて推定し、F検定で有意に差があるかを確認
r_weak_iv <- anova(
    lm(log(価格) ~ 1 + month + year, data = df02),
    lm(log(価格) ~ log(降水量.0) + log(日照時間.3) + month + year, data = df02)
)
print(paste("弱相関テスト", r_weak_iv["Pr(>F)"]))

[1] "弱相関テスト c(NA, 0.0305046662216007)"

4.3 過剰識別検定

TSLSで丁度識別ではなく、内生変数の数よりも操作変数の数が多い場合、誤差項と操作変数が相関してしまうときがあります。過剰識別です。LM統計量を用いて過剰識別ではないかχ二乗検定を行い確認します。

# 操作変数の数が多い場合は過剰識別になり、操作変数うちの一つと誤差項の相関が疑われる
r_oi_test <- lm(r_tsls$residuals ~ log(降水量.0) + log(日照時間.3) + month + year, data=df02)
p_oit <- 1 - pchisq(summary(r_oi_test)$r.squared*length(r_oi_test$residuals), 1)
print(sprintf("過剰識別検定 %0.5f", p_oit))

[1] "過剰識別検定 0.83026"

4.4 DWH検定

OLSでもTSLSでも推定結果が同じと言うことがあります。内生性がないか、操作変数も実は内生変数である場合です。DWH検定を用いて、推定結果に相違があるか確認します。

# DWH検定
#  まずOLSで推定し、次にOLSとTSLSで有意に差があるかを確認
#  棄却するとOLS推定量が一致推定量でないと言え、TSLS推定量を採用すべきとなる
r_ols <- lm(frml, data=df02) 
h_value <- t(r_tsls$coefficients - r_ols$coefficients) %*% solve(vcov(r_tsls) - vcov(r_ols)) %*% (r_tsls$coefficients - r_ols$coefficients)
p_dwh <- 1 - pchisq(h_value, 1)
print(sprintf("DWH検定 %0.5f", p_dwh))

[1] "DWH検定 0.07121"

5 制御関数法（CFM）

二段階残差算入（2SRI）、ヘックマン補正とも呼ばれる手法です。

TSLSが二段階目の推定で内生変数の予測値を内生変数の代わりに用いる一方で、2SRIは内生変数はそのままに一段階目の推定で得られた残差を説明変数に加えます。線形回帰においては、両者の推定される係数はほぼ同じになりますが、2SRIは非線形モデルにも使える利点があります。ロジットモデルと組み合わせる事が多いようです。

IVMの例のデータを用いて推定してみましょう。

cfm_1st <- lm(p ~ z, df01)
summary(cfm_2nd <- lm(d ~ p + residuals(cfm_1st), df01))


Call:
lm(formula = d ~ p + residuals(cfm_1st), data = df01)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.74052 -0.54686  0.00292  0.58100  2.20684 

Coefficients:
                    Estimate Std. Error t value Pr(>|t|)    
(Intercept)        100.11595    1.38950   72.05   <2e-16 ***
p                   -1.01101    0.06038  -16.75   <2e-16 ***
residuals(cfm_1st)   3.22561    0.10580   30.49   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.87 on 297 degrees of freedom
Multiple R-squared:  0.758, Adjusted R-squared:  0.7564 
F-statistic: 465.1 on 2 and 297 DF,  p-value: < 2.2e-16

標準誤差（分散共分散行列）は適切に推定されないので補正をするか、ブートストラッピングを用いる事が推奨されています（Palmer et al. (2017)）。ここではブートストラッピングを用います。

library(boot)
set.seed(603)
boot_results <- boot(data = df01, statistic = function(data, i){
  s <- data[i, ]
  cfm_1st <- lm(p ~ z, s)
  cfm_2nd <- lm(d ~ p + residuals(cfm_1st), s)
  coef(cfm_2nd)[1:2]
}, R = 1000)
(cfm_bootstrapping_stderr <- apply(boot_results$t, 2, sd))

[1] 3.1203543 0.1354724

過小ではありますが、IVMに近い、かなり改善した標準誤差が得られました。

説明変数	IVM	CFM-raw	CFM-Bootstrapping
Const.	3.2948058	1.3895014	3.1203543
p	0.1431626	0.0603752	0.1354724

6 完全情報最尤法（FIML）

尤度関数を用いることもできます。需要関数と均衡価格式を同時推定することを考えましょう。推定する係数を\(\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} \]

IVMの例のDGPの係数とは、\(\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\)の尤度函数から導かれます²。また、数理的な特性を使って整理すると以下のようなconcentrated loglikelihood function³にできることが知られています⁴。

\[\begin{align} \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 | \end{align}\]

\(g\)は内生変数（もしくは左辺にまとめられた被説明変）の数で、ここでは\(2\)です。\(\beta = ( \beta_{11}, \beta_{12}, \beta_{22}, -\beta_{31} )\)の4つのパラメーターを上の対数尤度函数で推定すればよいです。

n <- nrow(df01)
Y <- with(df01, cbind(d, p))
W <- with(df01, 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", hessian = TRUE)
(fiml_beta <- r_optim$par[c(1, 4)])

[1] 100.128351   1.011543

(fiml_se <- sqrt(diag(solve(r_optim$hessian)))[c(1, 4)])

[1] 3.2771152 0.1424937

推定量も標準誤差も概ねIVMと同じものが得られました。\(\beta_{31}\)は符号が逆に推定されることに注意してください。

最尤法は非線形モデルにも適用できるため、応用範囲は広いです。

天候は生産に大きく響くが、消費には大した影響はないと仮定してもよいと思います。↩︎
Davidson and MacKinnon (2003)もしくはオンラインで配布している2021年版のpp.503–505とpp.544–545を参照。誘導形はSURと同じ多変量正規分布の尤度函数になるので、\(\Gamma\)による変形の影響を整理すれば構造形の尤度函数が導出できます。↩︎
定訳不明。↩︎
一階条件の節の一階条件に直接関係のないExercise 12.24を\(\Sigma \Sigma^{-1}=I\)に注意して\(\Sigma^{-1}\)の成分を考えて解くことで得ることができるはずです。↩︎

操作変数法を学んで、計量経済学に踏み出そう

1 同時方程式

1.1 価格の内生性と同時性バイアス

2 二段階最小二乗法（TSLS）

3 操作変数法（IVM）

3.1 係数

3.1.1 一致性の確認

3.2 分散共分散行列

3.3 IVMの例

3.3.1 データセット

3.4 推定

3.4.1 誤った分散共分散行列

3.5 不均一分散

4 TSLSの実践

4.1 行列による計算

4.2 弱相関検定

4.3 過剰識別検定

4.4 DWH検定

5 制御関数法（CFM）

6 完全情報最尤法（FIML）