一般化モーメント法（GMM）は、計量経済学（Econometrics）で開発され発展した手法で、非線形モデルが推定でき、不均一分散、内生性に対処することができ、さらに誤差項が従う分布を特定しなくても推定ができる場合があるため、経済学の計量分析でよく使われています1。
社会科学分野の統計を見る事が多い人は、目にする機会がそこそこあるはずの手法なので、概要を把握しておくのが良いでしょう。推定手順を説明し、簡単なものを実践してみるための説明を用意しました。
1 GMMの構成
モーメント法を多数のモーメント条件を使って推定できるように一般化したのがGMMです。学部の統計学の教科書にも紹介があるモーメント法ですが、だいたいすっかり忘れ去られていると思うので、モーメント条件から確認していきましょう。数理的に詳しい説明が必要な場合は、「GMM — 計量経済学大学院講義ノート」に詳しい導出があったので、参考にしてください。
1.1 モーメント条件
モーメント条件は、期待値をとったらゼロになる観測値とパラメーターの関係です2。代表的なものを幾つか見ておきましょう。
1.1.1 よくある例（正規分布）
GMMの紹介では平均と分散と正規分布のモーメント特性を使った3つの式を例として挙げることが多いです。
\[\begin{align} E[x - \mu] &= 0 \\ E[(x - \mu)^2 - \sigma] &= 0 \\ E[(x - \mu)^4 - 3 \sigma^4] &= 0 \end{align}\]
\(x\)が観測値、\(\mu\)は母平均、\(\sigma^2\)は母分散です。左辺の期待値をモーメントと呼びます。
モーメント条件がパラメーターよりも多い場合があることを示すために3つになっています。正規分布は2変数なのでモーメント条件2つで推定でき、この2つで推定するとモーメント法になります。3つ使うためには、GMMが必要です。
1.1.2 実践的な例（線形モデル）
教科書のGMMの節では実践的に線形回帰モデルのモーメント条件を
\[ E[{}^t\!Z(y - X\beta)] = 0 \]
として説明していることが多いです。
\(Z\)は外生変数の行列、\(y\)は従属変数ベクトル、\(X\)は内生変数を含む説明変数の行列、\(\beta\)は係数ベクトルです。観測値の数を\(n\)、外生変数の数を\(l\)、説明変数の数を\(k\)とすると、\(Z\)は\(n\)行\(l\)列、\(y\)の長さは\(n\)、\(X\)は\(n\)行\(k\)列、\(\beta\)の長さは\(k\)です。\(y - X\beta\)は回帰モデルの残差なので誤差項になり、\(E[{}^t\!Z(y - X\beta)]\)がモーメントになります。
この式は、外生変数と誤差項が直交していることを表していますが、和がゼロになると言う直交条件は、期待値がゼロになると言う事なので、モーメント条件でもあります。また、式は1つですが、行列表記になっているので、モーメント条件は\(l\)あります。\(l \ge k\)であれば、分布に依存するモーメント条件なくして推定できます。
1.1.3 非線形の例（ロジットモデル）
非線形モデルでもモーメント条件が導出できれば推定できます。ロジットモデルだと以下のモーメント条件が使えます。
\[ E\bigg[{}^t\!Z\bigg( y - \frac{1}{1 + e^{-X\beta}}\bigg)\bigg] = 0 \]
分布を特定しているので有り難みが薄いですが、操作変数が使えるのは利点かも知れません。
1.2 ウェイト行列
丁度識別と言いますが、\(l = k\)であれば、すべてのモーメント条件を満たすパラメーターが存在します。しかし、\(l > k\)である場合は、すべての条件を満たせない事がほとんどです。
素朴な解決策として、すべてのモーメントの二乗和を足した値が一番小さくなるパラメーターを推定量とすることを考えます。非線形モデルでも同様に議論できますが、後で試すのが線形モデルなので、線形モデルのモーメント関数を以下のように\(g(y, X, Z, \beta)\)と置きます。
\[ g(y, X, Z, \beta) = E[ {}^t\!Z(y - X\beta) ] \]
教科書的な表記は\(g(\beta)\)ですが、後のRのコードにあわせています。
真のパラメーター\(\beta\)を任意の\(b\)に置き換え、期待値ではなく算術平均をとった標本モーメント関数を定義します3。
\[ \bar{g}(y, X, Z, b) = \frac{1}{n}\sum^n_{i=1}{}^t\!Z_i(y_i - X_ib) \]
これを使い係数の推定量を定義すると以下になります。
\[ \newcommand{\argmin}{\mathop{\rm arg~min}\limits} \hat{\beta}^{identity} = \argmin_{b} \ {}^t\! \bar{g}(y, X, Z, b) \bar{g}(y, X, Z, b) \]
この\(\hat{\beta}^{identity}\)は一致推定量になるのでそんなには悪く無いのですが、効率的ではないことが分かっています。\(l\)あるモーメントの分散は異なり、不均一分散があるからです。
GMMでは不均一分散を是正するために、\(l\)行\(l\)列のウェイト行列\(W\)を導入します。
\[ \hat{\beta} = \argmin_{b} \ {}^t\! \bar{g}(y, X, Z, b) W \bar{g}(y, X, Z, b) \]
この式は、GMM推定量の説明でよく見かけるもので、ここでおおよそ完成です。しかし、\(W\)をどう求めれば良いのか課題が残ります。恒等行列\(I\)を入れても推定できるのですが、それでは\(\hat{\beta}^{identity}\)が得られてしまうので改善になりません。
1.2.1 Two-Step GMM / Efficient GMM
効率的なウェイト行列の計算方法は分かっていて、外生変数の分散共分散行列の
\[ S = E[ \bar{g}(y_i, X_i, Z_i, \beta)\ {}^t\!\bar{g}(y_i, X_i, Z_i, \beta)] \]
の逆行列を用いるのが効率的であることが分かっています。\(\bar{g}(y_i, X_i, Z_i, \beta)\)は、\(i\)番目（行目）の観測値だけを使った\(\beta\)を中心とする標本モーメントで、\(1\)行\(l\)列になります。\(1\)行\(l\)列の行列に\(l\)行\(1\)列の行列を乗じているので、\(l\)行\(l\)列の行列になるのに注意してください。
しかしながら、\(\beta\)は未知の値なので\(S\)は計算できません。\(S\)の推定量を得る様々な手法が提案されていますが、応用ではTwo StepもしくはEfficient GMMと呼ばれる手法が一般的です。
1.2.1.1 非線形モデルにも使える方法
ウェイト行列\(W^0\)に恒等行列\(I\)を用いて、推定を行いone step推定量\(\hat{\beta}^{one}\)を得ます。次に、\[ \hat{S} = \frac{1}{n} \sum^n_{i=1} \bar{g}(y_i, X_i, Z_i, \hat{\beta}^{one})\ W^0 \ {}^t\!\bar{g}(y_i, X_i, Z_i, \hat{\beta}^{one}) \]
と\(\hat{S}\)を得た後、\(W^0=\hat{S}^{-1}\)として推定をしなおし、Efficient GMM推定量\(\hat{\beta}^{two}\)を得ます。
\[ \hat{\beta}^{two} = \argmin_{b} \ {}^t\! \bar{g}(y, X, Z, b) \hat{S}^{-1} \bar{g}(y, X, Z, b) \]
この手順を\(\hat{S}\)が収束したと見なせる値まで繰り返すこともでき、Iterated GMMと呼ばれます。
1.2.1.2 線形モデルに使える効率の良い方法
線形モデルの場合は\(W^0\)に\((Z{}^t\!Z)^{-1}\)の逆行列を使います。TSLS推定量に一致するので、TSLS-GMM推定量と呼ばれます。この方法は2つの意味で効率が良いです。
モーメントの標本分散共分散行列\(\hat{S}\)は\(S\)に収束するのですが、モーメント\(E[{}^t\!Z(y - Xb)]\)の\(y - Xb\)は誤差項です。\(Z{}^t\!Z\)の収束先と（定数である）誤差項の母分散を乗じたものが、\(S\)となります。誤差項の母分散は定数ですが、係数の点推定量を得るためのウェイト\(W^0\)は、どのモーメント条件を相対的にどの程度重視するか定めるものです。定数倍変化しても計算結果は同じになります。\(\hat{S}\)ではなく\(Z{}^t\!Z\)を用いても大差は生じません。むしろ、誤差項を乗じない\(Z\)の分散共分散行列の方が収束が早いです。つまり、\(Z\)の分散共分散行列を使う方が、統計学の意味で効率がよくなります。
計算量の意味でも効率がよいです。two stepと言いつつ最適化は一回で済みます。\(\hat{\beta}\)によって、\(Z{}^t\!Z\)は変化しないからです。係数の分散共分散行列を計算するには\(\hat{S}\)が必要になりますが、\(Z{}^t\!Z\)に誤差項の標本分散を乗じれば\(\hat{S}\)になります。
1.2.2 CU-GMM
CU-GMMは、ウェイト行列を係数の推定量の関数\(W(b)\)と見做して、係数と同時に推定する方法です。
1.2.3 Newey-West一致推定量
自己回帰項がある時系列データ向けのウェイト行列の作り方です。
1.3 過剰識別検定（J検定）
操作変数を積み上げるなどしてモーメント条件を大量に入れると、モーメント条件とモーメント条件が直交しないモデル特定化の失敗が生じ、推定結果が不安定になるなどの問題が生じるときがあります。
この問題が生じていないかを調べる過剰識別検定は、J統計量を
\[ J(b) = n \bar{g}(y, X, Z, b) \hat{S}^{-1} \bar{g}(y, X, Z, b) \sim \chi^2(l - k) \]
を用いて行います[^keiikegami]。ウェイト付き標本モーメントに標本サイズを乗じていますが、モーメントが十分にゼロに近いかを見る検定です。モデル定式化は適切であると言うのが帰無仮説なので、χ二乗検定を棄却してしまったときはモデル（操作変数）を見直す必要があります。
なお、\(l=k\)の丁度識別のときはモーメントはすべて\(0\)になるため、J統計量も\(0\)になりモデル定式化が不適切になることはありません。
また、GMM推定量はこのJ統計量を用いて定義する事もできます。
\[ \hat{\beta}^{GMM} = \argmin_{b}\ J(b) \]
1.4 分散共分散行列
統計学的仮説検定や信頼区間の計算に必要な推定されたパラメーターの分散共分散行列は、最小化する目的関数の一階微分がゼロになることと、\(\sqrt{n} \bar{g} \xrightarrow[d]{} \mathcal{N}(0, S)\)になることから求められます。
\[ V_\beta = ({}^t\!G W G)^{-1}({}^t\!G W S^{-1} W G) ({}^t\!G W G)^{-1} \]
ここで、
\[ G = E\bigg[ \frac{\partial}{\partial \beta} g(y, X, Z, \beta) \bigg] \]
です。推定するときは、\(S\)は\(\hat{S}\)を用い、\(G\)は標本平均\(\bar{G}\)を代わりに用います。
\[ \bar{G} = \frac{1}{n} \sum^n_{i=1} \frac{\partial}{\partial \beta} \bar{g}(y_i, X_i, Z_i, \beta) \]
線形モデルの推定においては、
\[ G = E\big[ {}^t\! Z X \big],\quad S = E\big[ Z {}^t\! Z e^2 \big] \]
を用いて計算します。Iterated GMMだと\(W=S^{-1}\)になるので、\(({}^t\!G W G)^{-1}\)まで計算が簡潔になります。
2 推定
教科書の説明と対応が分かるように簡単なモデルを推定した後、パッケージを使って同じ推定をします。
2.1 データセット
内生変数1つ、操作変数2つの線形モデルのデータセットをつくります。
2.2 推定
線形モデルとしてTSLS-GMM推定量を、非線形モデルとしてTwo Step GMM推定量を計算します。
まず、任意のパラメーター\(b\)からモーメントを計算する関数を用意します。
# データフレームの変数を展開
X <- model.matrix(~ x + z, df01)
Z <- model.matrix(~ z + d1 + d2, df01)
y <- df01$y
n <- nrow(df01)
# モーメント関数を定義
g <- function(y, X, Z, b) (t(Z) %*% (y - X %*% b))/length(y)
# 試しに計算してみる
g(y, X, Z, c(1, 1, -1)) [,1]
(Intercept) 0.5473707
z -0.6451392
d1 -0.0934314
d2 0.56377332.2.1 線形モデルとして推定
ニュートン=ラフソン法などの最適化アルゴリズムで、y,
X,
Zを所与として、gを最小化するbを求めます。
# 最小化する目的関数
objf <- function(b, W) t(g(y, X, Z, b)) %*% W %*% g(y, X, Z, b)
# Step 1
W0_tsls <- solve(t(Z) %*% Z / n)
r_tsls <- nlm(objf, c(0, 0, 0), W = W0_tsls, gradtol = 1e-10, steptol = 1e-10) # momentfitパッケージの出力との比較のために制度をあげている
(b_tsls <- r_tsls$estimate)[1] 1.181857 1.947318 -2.946895# 誤差項の影響をウェイトに加える
include_variance <- function(W, b){
e <- y - X %*% b
W / (drop(t(e) %*% e) / n)
}
# Step 2はウェイト行列の補正だけ
W_tsls <- include_variance(W0_tsls, b_tsls)
k <- ncol(X)
l <- ncol(Z)
# J検定
J <- function(b, W) n * t(g(y, X, Z, b)) %*% W %*% g(y, X, Z, b)
sprintf("J-Statistics (One Step): %5f", (J_tsls <- J(b_tsls, W_tsls)))[1] "J-Statistics (One Step): 0.222662"[1] "p-value of J-test: 0.63702"TSLS-GMM推定量が計算できました。係数の分散共分散行列も出しておきましょう。
G <- t(Z)%*%X/n
inv_GWG <- solve(t(G) %*% W0_tsls %*% G)
S_hat_tsls <- solve(W_tsls) / n
(V_tsls <- inv_GWG %*% (t(G) %*% W0_tsls %*% S_hat_tsls %*% W0_tsls %*% G) %*% inv_GWG) (Intercept) x z
(Intercept) 0.059537083 -0.03871773 0.008202819
x -0.038717727 0.03490428 -0.015916071
z 0.008202819 -0.01591607 0.028152453概ねDGPの係数\(\beta = {}^t\!(1, 2, -3)\)に近い値が推定できています。また、モーメント条件の妥当性は\(J\)検定で棄却されず、信頼がおけると言うことになります。
2.2.2 非線形モデルとして推定
\(W^0\)を恒等行列にして推定します。
W0_nl <- diag(ncol(Z))
r_nlm_step1 <- nlm(objf, c(0, 0, 0), W = W0_nl, gradtol = 1e-10, steptol = 1e-10)
(b_step1 <- r_nlm_step1$estimate)[1] 1.152770 1.978351 -2.973269続けてモーメントの分散共分散行列の\(\hat{S}\)を求めます。
estimate_S <- function(b){
g_i <- matrix(NA, n, l)
# 1行ごとのモーメントを計算
for(i in 1:n){
g_i[i, ] <- g(y[i], X[i, ,drop = FALSE], Z[i, ,drop = FALSE], b)
}
# 列（外生変数）ごとの平均を出す
g_mean <- apply(g_i, 2, mean)
# モーメントの分散共分散行列を計算
S_hat <- matrix(0.0, l, l)
for(i in 1:n){
# Centeredしているが、無くてもよい
S_hat <- S_hat + (g_i[i, ] - g_mean) %*% t(g_i[i, ] - g_mean) / n
}
S_hat
}
S_hat <- estimate_S(b_step1)\(\hat{S}\)の逆行列をウェイトとします。
推定したウェイトを用いて、2段階目の推定を行います。
r_nlm_step2 <- nlm(objf, b_step1, W = W_step2, gradtol = 1e-10, steptol = 1e-10) # b_step1は初期値に使っているだけ
(b_step2 <- r_nlm_step2$estimate)[1] 1.152502 1.963083 -2.951443J推定量を出します。
[1] "J-Statistics (Two Step): 0.194305"[1] "p-value of J-test: 0.65936"推定結果自体は概ね良好ですが、TSLS-GMM推定量とは異なる値になりました。
分散共分散行列は、非線形であることを前提にするので、数値微分した値を用いて計算します。
# Zの要素ごとのモーメントをXの要素で数値微分
derivX <- function(y, X, Z, b){
G <- matrix(NA, ncol(Z), ncol(X),
dimnames = list(colnames(Z), colnames(X)))
for(j in 1:ncol(X)){
d <- numeric(ncol(X))
delta <- 1e-6
d[j] <- delta
g0 <- g(y, X, Z, b - d)
g2 <- g(y, X, Z, b + d)
G[,j] <- (g2 - g0)/2/delta
}
G
}
G <- derivX(y, X, Z, b_step1)
inv_GWG <- solve(t(G) %*% W_step2 %*% G)
S_hat_step2 <- estimate_S(b_step2) / n
(V_2step <- inv_GWG %*% (t(G) %*% W_step2 %*% S_hat_step2 %*% W_step2 %*% G) %*% inv_GWG) (Intercept) x z
(Intercept) 0.05199936 -0.03823610 0.01282015
x -0.03823610 0.03866967 -0.01931077
z 0.01282015 -0.01931077 0.022537762.3
momentfitを使った推定
最尤法と比べてコードの分量は多くはないですが、やはりパッケージを使って推定する方が無難です。推定だけならばまだしも、制約の検定などもあります。momentfitパッケージ4を使ってみましょう。
2.3.1 TSLS推定量
線形モデルはsem::tslsのように推定できます。
要求されたパッケージ sandwich をロード中ですmodel_lm <- momentModel(y ~ x + z, ~ z + d1 + d2, data = df01)
r_model_lm <- gmmFit(model_lm, init = "tsls")
summary(r_model_lm)Model based on moment conditions
*********************************
Moment type: linear
Covariance matrix: iid
Number of regressors: 3
Number of moment conditions: 4
Number of Endogenous Variables: 1
Sample size: 100
Estimation: Two-Stage Least Squares
Sandwich vcov: FALSE
coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.18186 0.24400 4.8436 1.275e-06 ***
x 1.94732 0.18683 10.4231 < 2.2e-16 ***
z -2.94689 0.16779 -17.5633 < 2.2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
J-Test
Statistics df pvalue
Test E(g)=0: 0.22266 1 0.63702
Instrument strength based on the F-Statistics of the first stage OLS
x : F( 2 , 96 ) = 55.76589 (P-Vavue = 1.110223e-16 ) (Intercept) x z
(Intercept) 0.059537083 -0.03871773 0.008202819
x -0.038717727 0.03490428 -0.015916071
z 0.008202819 -0.01591607 0.028152453
attr(,"type")
attr(,"type")$sandwich
[1] FALSE
attr(,"type")$df.adj
[1] FALSE
attr(,"type")$breadOnly
[1] FALSE簡単に計算できました。パッケージを使わない場合と同じ値ですね。
2.3.2 TWO-STEP GMM推定量
非線形モデルの推定の場合、行ごとのモーメントを計算する関数を用意します。
moment_function <- function(b, M) {
y <- M[, 1]
X <- M[, 2:(1+ncol(X))]
Z <- M[, (2+ncol(X)):(1+ncol(X)+ncol(Z))]
e <- drop(y - X %*% b)
Z * e # Zの各列にeを乗じる
}
model_nl <- momentModel(
g = moment_function,
x = cbind(y, X, Z),
theta0 = c(1, 1, 1))
r_model_nl <- gmmFit(model_nl)
summary(r_model_nl)Model based on moment conditions
*********************************
Moment type: function
Covariance matrix: iid
Number of regressors: 3
Number of moment conditions: 4
Number of Endogenous Variables: 0
Sample size: 100
Estimation: Two-Step GMM
Convergence code: 0 (see help(optim))
Sandwich vcov: FALSE
coefficients:
Estimate Std. Error t value Pr(>|t|)
theta1 1.15250 0.22803 5.0541 4.325e-07 ***
theta2 1.96308 0.19665 9.9828 < 2.2e-16 ***
theta3 -2.95144 0.15013 -19.6598 < 2.2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
J-Test
Statistics df pvalue
Test E(g)=0: 0.1943 1 0.65936 theta1 theta2 theta3
theta1 0.05199935 -0.03823610 0.01282015
theta2 -0.03823610 0.03866966 -0.01931076
theta3 0.01282015 -0.01931076 0.02253775
attr(,"type")
attr(,"type")$sandwich
[1] FALSE
attr(,"type")$df.adj
[1] FALSE
attr(,"type")$breadOnly
[1] FALSE簡単に計算できました。パッケージを使わない場合と同じ値ですね。
2.3.3 制約付きモデル
momentfitでは、パラメーターを\(\theta\)として、\(R\theta =
q\)という形式の制約をつけることができます。
\(x = 0\)と言う制約を置いた線形モデルを推定してみましょう。
Model based on moment conditions
*********************************
Moment type: rlinear
Covariance matrix: iid
Number of regressors: 2
Number of moment conditions: 4
Number of Endogenous Variables: 1
Sample size: 100
Constraints:
x = 0
Restricted regression:
y = (Intercept)+z Model based on moment conditions
*********************************
Moment type: rlinear
Covariance matrix: iid
Number of regressors: 2
Number of moment conditions: 4
Number of Endogenous Variables: 1
Sample size: 100
Constraints:
x = 0
Restricted regression:
y = (Intercept)+z
Estimation: Two-Stage Least Squares
coefficients:
(Intercept) z
3.341928 -2.058934 \(x = 1 - z\)と言う制約を置いた非線形モデルを推定してみましょう。
Model based on moment conditions
*********************************
Moment type: rfunction
Covariance matrix: iid
Number of regressors: 2
Number of moment conditions: 4
Number of Endogenous Variables: 0
Sample size: 100
Constraints:
theta2 ~ 1 - theta3Model based on moment conditions
*********************************
Moment type: rfunction
Covariance matrix: iid
Number of regressors: 2
Number of moment conditions: 4
Number of Endogenous Variables: 0
Sample size: 100
Constraints:
theta2 ~ 1 - theta3
Estimation: Two-Step GMM
Convergence code: 0 (see help(optim))
coefficients:
theta1 theta3
-1.376305 -3.021013 変数名がtheta1, theta2,
theta3となってしまうこと、~を挟んで左辺と右辺に値が必要なので注意してください。
2.3.4 ウェイトを指定した推定
momentfitを使う上でウェイトを強く意識する必要はないと思いますが、ウェイトを指定した推定もできます。
wObj_tsls <- evalWeights(model_lm, w = solve(t(Z) %*% Z))
(theta_tsls <- solveGmm(model_lm, wObj_tsls)$theta)(Intercept) x z
1.181857 1.947318 -2.946895 (Intercept) x z
1.152720 1.978390 -2.973282 wObj2 <- evalWeights(model_lm, theta = r_solveGmm1$theta, w = "optimal")
(solveGmm(model_lm, wObj2)$theta)(Intercept) x z
1.181857 1.947318 -2.946895 3 まとめ
GMMは、最尤法と比較すると、分布を特定しなくても済む場合がある一方、効率性で劣る事が知られています。また、同じモデルでもウェイト行列の置き方によって推定にバリエーションができるため、とっつきにくさがあります。しかし、資産価格や投資のEuler方程式の推定や、動学的パネルデータモデルの推定ではよく使われています。
経済学以外で見かける事はほとんどありません。社会学で明らかに経済学の手法を利用した場合ぐらいです。歴史的な理由と、最尤法などに対して圧倒的に優位と言う事ではないためだと思います。↩︎
物理学用語のモーメントと概念がややずれるので、このように把握する方が混乱が少ないと思います。↩︎
\(\bar{g}(b)\)と書くのが通例ですが、後の説明とRのコードにあわせるために冗長に書いています。↩︎
momentfitは同じ作者のgmmパッケージの後継で、vignetteに計量経済学の有名テキストの例を実行する方法の紹介があるので、GMMを研究に用いた人はbrowseVignettes("momentfit")で添付のドキュメントに目を通しておくと良いと思います。↩︎