こんにちは、リブセンスでデータサイエンティストをしている北原です。今回はRのmomentfitパッケージを使ってGeneralized Method of Moments(GMM、一般化モーメント法)を実行する方法について紹介します。
GMMはパラメータ推定法の一つで、主に計量経済学で使われています。実務でGMMを直接使うケースはまれだと思いますが、因果推論に関わる機械学習で使われることがあるため、理解しておきたい手法の一つでもあります。今回は線形回帰、操作変数法、ロジスティック回帰を例として実際にGMMでパラメータ推定を行い、GMMをデータ分析や機械学習で利用する際の利点・欠点について考えます。なお、後日公開予定のCovariate Balancing Propensity Scoreの記事では今回紹介するGMMを利用します。
Generalized Method of Moments
GMMは、モーメント条件の距離を最小化することで、未知パラメータの数よりモーメント条件の数が多いときにも使えるようモーメント法を拡張したパラメータ推定法です。モーメント法の拡張であることからわかるように、モーメント条件がわかれば分布を仮定することなく使うことができます。なお、モーメント法は母集団のモーメントを標本モーメントに置き換えることによりパラメータを推定する方法で、計量経済学に関わらず統計学ではよく知られた手法です。確率分布がわからないあるいは計算しにくいがモーメントは計算しやすいときに使われます。モーメントは確率変数のべき乗の期待値ですが、べき乗でない非線形関数の場合も関数を展開すればべき乗になるためか、確率変数がべき乗以外の非線形関数の期待値になっていてもモーメント条件として使われます。
GMMを使うと操作変数法などの計量経済学で使われる推定法を統一的に扱うことができます。そのため計量経済学ではよく知られた手法ですが、それ以外の分野ではあまり知られていないかもしれません。「一般化モーメント法」などのキーワードで検索すれば講義資料や解説記事が見つかるので詳細についてはそちらをご覧ください。なお「GMM」で検索するとガウス混合モデル（Gaussian Mixture Model）がヒットしますが本記事で扱うGMMとは無関係なのでご注意ください。
GMMのポイントは、未知パラメータの数よりモーメント条件の数が多いときに使えるところにあります。モーメント条件の数が多すぎるときでも一部の条件を使わず未知パラメータと同数の条件だけを使えば通常のモーメント法で推定ができます。しかし、使われなかった条件に重要な情報が含まれていた場合この方法は適切ではありません。GMMでは、モーメント条件の距離最小化によって、考慮すべき全てのモーメント条件をなるべく満たすようにパラメータを決めることで、この問題に対処しています。
サンプル数を$N$、未知パラメータ数を$K$、モーメント条件数を$L$とし、未知パラメータ${\mathbf \beta} \in {\mathbb R}^K$とサンプルデータ${\mathbf x} \in {\mathbb R}^{N \times K}$、$L$個のモーメント条件${\mathbf g }({\mathbf \beta}, {\mathbf x})$を使って標本モーメントを
と表記します。ここで$K \le L$です。$L \times L$の正定値行列${\mathbf W}$を使って、GMMでは
として、２次形式で表現される距離を最小にすることで未知パラメータを計算します。線形モデルの場合は解析解が得られますが非線形モデルの場合は数値計算が必要になります。
${\mathbf W}$は重み行列で、どのモーメント条件を重要視するかを表したものです。例えば単位行列を使うと全てのモーメント条件を等しく０に近づけるようにパラメータを決めることになります。どのような${\mathbf W}$を使っても一致性は満たされていますが分散は異なるので分散が小さくなる${\mathbf W}$を使うことが望ましいです。${\mathbf W}=V[{\mathbf g }({\mathbf \beta}, {\mathbf x})]^{-1}$のように分散の逆行列を使うことで推定量の分散が漸近的に最小になることが知られています。$V[{\mathbf g }({\mathbf \beta}, {\mathbf x})]^{-1}$は未知ですが推定量に置き換えて計算することができるので通常はこれを使います。標本分散が大きいモーメント条件の重みを小さくするというのは直感的にも納得しやすいのではないかと思います。
J検定によってモーメント条件が満たされているかを検定することができます。パラメータ数よりモーメント条件数が多い場合、データによってはモーメント条件を満たすパラメータが存在しないこともあります。そのような場合はJ検定を使うことでモーメント条件が適切でないことを検出できます。J検定では
を帰無仮説とします。そのため、p値が小さく帰無仮説が棄却された場合にモーメント条件は満たされないと判断することになります。t検定などの統計的仮説検定ではp値が小さいと意味のある分析結果になることが多いですが、J検定では逆でp値が小さくないときに意味のある結果になります。なお、通常の統計的仮説検定と同じくデータ数が増えると帰無仮説が棄却されやすくなるので注意が必要です。
momentfitパッケージ
momentfitパッケージはGMMを使ってパラメータ推定を行うRパッケージです。GMMを行うパッケージとしてはgmmパッケージがあったのですが、その後継です。モーメント条件とデータを与えるとGMMを実行できるところは同じですが、momentfitでは回帰モデルの扱いやすさ向上などの様々な改良や機能拡張がされています。機能が豊富なため本記事では一部しか紹介できません。使い方の詳細についてはvignettesをご参照ください。また、gmmパッケージのvignettesも参考になると思います。
momentfitパッケージでモデルを指定する方法には大きく分けてformulaを使う方法とモーメント条件を計算する関数を使う方法の２種類があります。
formulaを使う場合、線形モデルではlm()関数などとほぼ同じ表記が使えます。非線形モデルではglm()関数などと同じ表記は使えないものの、回帰係数と変数を使って関数を表した直感的にわかりやすい表記が使えます。そのため通常のRパッケージを利用する感覚で使うことができます。このように通常の回帰モデルの場合はformulaを使えるのですが、他のモーメント条件が必要な場合は使えないので汎用性は低いです。また、formulaを使うと数値計算ではなく解析解を利用した計算が行われます。なお、操作変数法を前提として設計されているのか、基本的には操作変数のモデル式も指定する必要があります。
モーメント条件を計算する関数を使う場合、推定するパラメータとデータを引数とし標本モーメント（より正確には平均する前の各サンプルの計算結果）を計算して返す関数でモデルを表します。明示的にモーメント条件を扱う必要がありますが、標本モーメントを計算できるモデルであればパラメータ推定ができます。
本記事で扱うモデルは全てformulaを使う方法で実行できるのですが、今後様々なモデルを扱えるようにモーメント条件の関数を使う方法を本記事では紹介します。なお、後日公開予定のCovariate Balancing Propensity Scoreの記事ではモーメント条件の関数を使う方法が使われます。
GMMによる推定
本記事では、線形回帰、操作変数法、ロジスティック回帰のパラメータ推定について扱います。線形回帰と操作変数法は線形モデルの例、ロジスティック回帰は非線形モデルの例になっています。また、線形回帰とロジスティック回帰はパラメータ数とモーメント条件数が同数の例、操作変数法はパラメータ数よりモーメント条件数が多いときの例になっています。GMMを使っていると言えるのは操作変数法のときだけですが、モーメント法でもGMMでもモーメント条件を指定して推定するという流れは同じです。
いずれのモデルも既存の関数やパッケージで推定できます。線形回帰はlm()関数、操作変数法は例えばAERパッケージのivreg()関数、ロジスティック回帰はglm()関数で実行できます。GMMの推定結果の確認のために、これらの出力も見ます。
必要なパッケージは事前にインストールし読み込んでおきます。
library(tidyverse) library(momentfit) library(AER)
線形回帰
まず簡単な例として線形回帰のパラメータ推定から見ていきましょう。
サンプルデータは次のようなものを使います。
# データ作成 set.seed(1) n <- 1000 x <- tibble( x1 = rnorm(n, mean = 2, sd = 1), x2 = rnorm(n, mean = 3, sd = 1), x3 = rnorm(n, mean = 4, sd = 1)) %>% mutate(y = rnorm(n, mean = 1 + 2 * x1 + 3 * x2 + 4 * x3, sd = 1))
モデルは
とします。${\mathbf \beta}$は回帰係数、$y$は目的変数、${\mathbf x}$は説明変数、$u$は誤差です。定数項に対応するデータはすべて1になっているものとします。
さらに説明変数は全て外生変数とし誤差項との相関がないとします。説明変数と誤差項の相関が0なので
となります。これが線形回帰モデルのモーメント条件になります。分布を仮定していないところがポイントです。
標本モーメントを計算する関数のRコードは次のようになります。データの１列目は目的変数、２列目以降は説明変数になっていることを仮定しています。
# モーメント条件の関数 g <- function(theta, X_) { y <- X_[, 1] X <- X_[, -1] u <- c(y - X %*% theta) return(u * X) }
momentfitではモデルとデータの指定を行ってモデルのオブジェクトを作成し、そのオブジェクトを使って推定を行います。モデルのオブジェクトはmomentModel()関数を使って次のように作成します。
# モデルオブジェクト作成 model_gmm_lm <- momentModel(g = g, x = x %>% mutate(intercept = 1) %>% # 定数項に対応するデータはすべて1 select(y, intercept, x1, x2, x3) %>% as.matrix(), theta0 = runif(4))
gには標本モーメントを計算する関数、xにはデータ、theta0には数値計算で使う初期値を指定します。線形モデルでは初期値はあまり気にしなくても問題ないです。推定にはgmmFit()を使うのが便利です。
# GMM実行 fit_gmm_lm <- gmmFit(model_gmm_lm) summary(fit_gmm_lm)
実行すると次のような出力が得られます。
Model based on moment conditions
*********************************
Moment type: function
Covariance matrix: iid
Number of regressors: 4
Number of moment conditions: 4
Number of Endogenous Variables: 0
Sample size: 1000
Estimation: One-Step, Just-Identified
Convergence Optim: 0
Sandwich vcov: TRUE
coefficients:
Estimate Std. Error t value Pr(>|t|)
theta1 0.888519 0.176705 5.0282 4.95e-07 ***
theta2 2.023003 0.032981 61.3387 < 2.2e-16 ***
theta3 3.015163 0.031553 95.5586 < 2.2e-16 ***
theta4 4.009270 0.032919 121.7913 < 2.2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
J-Test
Statistics df pvalue
Test E(g)=0: 4.5351e-21 0 NA
coefficientsのところに推定結果が表示されるので確認すると、データ作成時に指定した数値とほぼ一致していることがわかります。パラメータ数とモーメント条件数が同じなのでJ検定は無視して構いません。
念のためlm()の結果も確認します。
model_lm <- lm(y ~ x1 + x2 + x3, x) summary(model_lm)
出力は以下のようになります。推定結果はGMMとほぼ同じですが誤差がわずかに小さいことがわかります。
Call:
lm(formula = y ~ x1 + x2 + x3, data = x)
Residuals:
Min 1Q Median 3Q Max
-3.2427 -0.7139 -0.0122 0.7199 3.0592
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.88852 0.17046 5.213 2.26e-07 ***
x1 2.02300 0.03183 63.564 < 2e-16 ***
x2 3.01516 0.03164 95.295 < 2e-16 ***
x3 4.00927 0.03195 125.482 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 1.04 on 996 degrees of freedom
Multiple R-squared: 0.9682, Adjusted R-squared: 0.9681
F-statistic: 1.011e+04 on 3 and 996 DF, p-value: < 2.2e-16
なお、formulaを使ってGMMを実行するには次のようにします。
# formulaを使ったGMM model_gmm_lm_formula <- momentModel( y ~ x1 + x2 + x3, ~ x1 + x2 + x3, data = x %>% dplyr::select(y, x1, x2, x3) ) fit_gmm_lm_formula <- gmmFit(model_gmm_lm_formula) summary(fit_gmm_lm_formula)
モーメント条件の関数や初期値の指定が不要になるので表記が簡単になっています。ただし、２番目の引数にはすべての外生変数を指定する必要があり、指定しないとエラーになります。
操作変数法
線形回帰モデルでは説明変数は全て外生変数としてモーメント条件を決めていたので、誤差項と相関のある変数である内生変数が含まれていると使えません。例えば目的変数にも説明変数にも影響を与える観測できない交絡変数が存在するとこのような問題が生じます。このようなときに使う方法の一つが操作変数法です。
操作変数法では、説明変数との相関はあるが誤差項とは相関のない変数である操作変数を利用する推定方法です。例として観測できない交絡変数があるときについて考えましょう。目的変数は説明変数と交絡変数の双方の影響を受けて変化するため、説明変数を変化させたときの目的変数の変化を調べても説明変数単独の影響がわかりません。しかし、操作変数を変化させたときの説明変数の変化には交絡変数の影響が含まれていないので、この変化と目的変数の変化の関係を調べれば説明変数単独の影響を取り出すことができます。
サンプルデータは以下のようなものを使います。x1は外生変数ですが、x2は内生変数です。z1とz2は操作変数です。
set.seed(1) n <- 1000 x <- tibble( x1 = rnorm(n, mean = 2, sd = 1), # 外生変数 cv = rnorm(n, mean = 0, sd = 1), # 交絡変数 z1 = rnorm(n, mean = 4, sd = 1), # 操作変数 z2 = rnorm(n, mean = 5, sd = 1) # 操作変数 ) %>% mutate( # 内生変数 x2 = rnorm(n, mean = 50 * cv + 30 * z1 + 30 * z2, sd = 1) ) %>% mutate( y = rnorm(n, mean = 1 + 10 * x1 + 3 * x2 + 50 * cv, sd = 1))
z1とz2はx2との相関はありますが、誤差項との相関はないので操作変数としての条件を満たしています。また、x2は内生変数なので誤差項との相関があります。これらのことは以下で確認できます。
# 内生変数x2と操作変数z1, z2との相関 cor(x$x2 , x) %>% as_tibble() %>% dplyr::select(z1, z2) # 誤差項とx1, x2, z1, z2との相関 cor(x$y - (1 + 10 * x$x1 + 3 * x$x2) , x) %>% as_tibble() %>% dplyr::select(x1, x2, z1, z2)
# A tibble: 1 x 2
z1 z2
<dbl> <dbl>
1 0.468 0.468
# A tibble: 1 x 4
x1 x2 z1 z2
<dbl> <dbl> <dbl> <dbl>
1 0.00604 0.769 0.0230 0.0138
まず通常の線形回帰の結果を見てみましょう。
# 通常の線形回帰 fit_lm <- lm(y ~ x1 + x2, data = x) summary(fit_lm)
出力の回帰係数を確認するとデータ作成時に指定した値とは少し異なっており、バイアスが生じていることがわかります。
Call:
lm(formula = y ~ x1 + x2, data = x)
Residuals:
Min 1Q Median 3Q Max
-98.564 -21.704 -0.561 22.386 98.216
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -154.26208 4.64070 -33.24 <2e-16 ***
x1 8.87271 1.01747 8.72 <2e-16 ***
x2 3.58002 0.01524 234.87 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 33.26 on 997 degrees of freedom
Multiple R-squared: 0.9823, Adjusted R-squared: 0.9823
F-statistic: 2.773e+04 on 2 and 997 DF, p-value: < 2.2e-16
それではGMMを使って推定しましょう。操作変数法のモデルは線形回帰モデル($\ref{eq:lm_model}$)と同じです。モーメント条件は外生変数を${\mathbf z}$としたときに
となります。線形回帰と同じく外生変数の定義をそのまま表したものになっています。この例ではx1とz1、z2が外生変数なので${\mathbf z}$になります。推定するパラメータは３つ、モーメント条件は４つなのでモーメント法ではなくGMMが必要なケースです。
標本モーメントを計算する関数や実行コードは以下のようになります。デフォルトでは分散を最小化した重みを使い２段階GMMが実行されます。
# モーメント条件関数 g <- function(theta, X_) { y <- X_[, 1] X <- X_[, 2:4] # intercept, x1, x2 Z <- X_[, c(2, 3, 5, 6)] # intercept, x1, z1, z2 u <- c(y - X %*% theta) return(u * Z) } # GMM実行 model_gmm_iv <- momentModel(g = g, x %>% mutate(intercept = 1) %>% dplyr::select(y, intercept, x1, x2, z1, z2) %>% as.matrix(), theta0 = runif(3)) fit_gmm_iv <- gmmFit(model_gmm_iv) summary(fit_gmm_iv)
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: 1000
Estimation: Two-Step GMM
Convergence Optim: 0
Sandwich vcov: FALSE
coefficients:
Estimate Std. Error t value Pr(>|t|)
theta1 -8.300325 9.938532 -0.8352 0.4036
theta2 10.231142 1.598819 6.3992 1.562e-10 ***
theta3 3.029649 0.034677 87.3672 < 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.045354 1 0.83135
パラメータの推定結果を見ると線形回帰よりバイアスが小さいことが確認できます。また、J検定の結果を確認するとp値が0.83135で帰無仮説が棄却されていないことがわかります。つまり、モーメント条件が満たされていないわけではないので、操作変数法はうまくいっていると考えてよいでしょう。
念のためにAERパッケージのivreg()関数で操作変数法を行ってみましょう。
library(AER) fit_ivreg <- ivreg(y ~ x1 + x2 | x1 + z1 + z2, data = x) summary(fit_ivreg)
Call:
ivreg(formula = y ~ x1 + x2 | x1 + z1 + z2, data = x)
Residuals:
Min 1Q Median 3Q Max
-156.828 -33.202 -1.226 36.242 175.780
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -8.23095 9.95995 -0.826 0.409
x1 10.23126 1.54740 6.612 6.17e-11 ***
x2 3.02942 0.03521 86.037 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 50.54 on 997 degrees of freedom
Multiple R-Squared: 0.9592, Adjusted R-squared: 0.9592
Wald test: 3767 on 2 and 997 DF, p-value: < 2.2e-16
GMMで推定したときとほぼ同じ結果になっていることが確認できます。また、formulaを使ってGMMを行う場合は次のようにします。
# formulaを使ったGMM model_gmm_iv_formula <- momentModel( y ~ x1 + x2, ~ x1 + z1 + z2, data = x %>% dplyr::select(y, x1, x2, z1, z2)) fit_gmm_iv_formula <- gmmFit(model_gmm_iv_formula) summary(fit_gmm_iv_formula)
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: 1000
Estimation: Two-Stage Least Squares
Sandwich vcov: FALSE
coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -8.230947 9.944995 -0.8276 0.4079
x1 10.231263 1.545081 6.6218 3.548e-11 ***
x2 3.029417 0.035158 86.1664 < 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.043875 1 0.83409
Instrument strength based on the F-Statistics of the first stage OLS
x2 : F( 2 , 996 ) = 379.7742 (P-Vavue = 0 )
ロジスティック回帰
GMMはモーメント条件がわかれば非線形モデルのパラメータ推定もできます。ここではロジスティック回帰のパラメータ推定を行います。
ロジスティック回帰モデルは、シグモイド関数$\sigma(x)=1/(1 + \exp(-x))$を使って
となります。ロジスティック回帰モデルのモーメント条件には
を使います。なお、標本モーメントは、ロジスティック回帰モデルの対数尤度を$\beta$で微分して0と置いたものと等しくなります。
サンプルデータには次のものを使います。
set.seed(1) n <- 1000 x <- tibble( x1 = rnorm(n, mean = 2, sd = 1), x2 = rnorm(n, mean = 1, sd = 1)) %>% mutate(y = rnorm(n, mean = 4 + 2 * x1 - 5 * x2, sd = 2)) %>% mutate(y = ifelse(1 / (1 + exp(-y)) >= 0.5, 1, 0))
モーメント条件がわかれば計算方法は今までと同じです。標本モーメントを計算する関数や実行コードは以下のようになります。
# モーメント条件関数 g <- function(theta, X_) { y <- X_[, 1] X <- X_[, -1] p <- 1 / (1 + exp(-X %*% theta)) u <- c(y - p) return(u * X) } # GMM実行 model_gmm_logistic <- momentModel( g = g, x = x %>% mutate(intercetp = 1) %>% select(y, intercetp, x1, x2) %>% as.matrix(), theta0 = c(1, 1, 1)) fit_gmm_logistic <- gmmFit(model_gmm_logistic) summary(fit_gmm_logistic)
Model based on moment conditions
*********************************
Moment type: function
Covariance matrix: iid
Number of regressors: 3
Number of moment conditions: 3
Number of Endogenous Variables: 0
Sample size: 1000
Estimation: One-Step, Just-Identified
Convergence Optim: 0
Sandwich vcov: TRUE
coefficients:
Estimate Std. Error t value Pr(>|t|)
theta1 3.67658 0.38262 9.609 < 2.2e-16 ***
theta2 1.84935 0.17974 10.289 < 2.2e-16 ***
theta3 -4.51761 0.32962 -13.706 < 2.2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
J-Test
Statistics df pvalue
Test E(g)=0: 8.5131e-13 0 NA
データ作成で指定した値と近い推定結果が得られていることがわかります。ここでは初期値にすべて1を指定していますが、他の初期値を使うと異なる結果が得られることがあります。モデルやデータに依存するのかもしれませんが初期値依存性が強いので、非線形モデルでGMMを使う場合は注意が必要です。最尤法などの推定値を初期値として使うことが多いようです。
念のためにglm()関数でロジスティック回帰を行うと
fit_glm <- glm(y ~ x1 + x2, x, family = "binomial") summary(fit_glm)
Call:
glm(formula = y ~ x1 + x2, family = "binomial", data = x)
Deviance Residuals:
Min 1Q Median 3Q Max
-2.75957 -0.05933 0.04549 0.24290 2.47470
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.6766 0.3855 9.536 <2e-16 ***
x1 1.8493 0.1868 9.901 <2e-16 ***
x2 -4.5176 0.3406 -13.263 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 1206.09 on 999 degrees of freedom
Residual deviance: 403.93 on 997 degrees of freedom
AIC: 409.93
Number of Fisher Scoring iterations: 7
となり、GMMとほぼ同じ結果が得られていることがわかります。formulaを使ってGMMを行う場合は次のようにします。非線形モデルの場合、回帰係数も含めてモデル式を指定し、初期値theta0の名前をモデル式の回帰係数名と同じにする必要があります。
# formulaを使ったGMM theta0 <- c(theta1=1, theta2=1, theta3=1) model_gmm_logistic_formula <- momentModel( y ~ 1 / (1 + exp(-(theta1 + theta2 * x1 + theta3 * x2))), ~ x1 + x2, data = x %>% dplyr::select(y, x1, x2), theta0 = theta0) fit_gmm_logistic_formula <- gmmFit(model_gmm_logistic_formula) summary(fit_gmm_logistic_formula)
Model based on moment conditions
*********************************
Moment type: nonlinear
Covariance matrix: iid
Number of regressors: 3
Number of moment conditions: 3
Number of Endogenous Variables: 1
Sample size: 1000
Estimation: One-Step, Just-Identified
Convergence Optim: 0
Sandwich vcov: TRUE
coefficients:
Estimate Std. Error t value Pr(>|t|)
theta1 3.67658 0.90931 4.0433 5.271e-05 ***
theta2 1.84935 0.51377 3.5996 0.0003187 ***
theta3 -4.51761 1.14864 -3.9330 8.388e-05 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
J-Test
Statistics df pvalue
Test E(g)=0: 8.3205e-13 0 NA
まとめ
今回はGMMを使って線形回帰、操作変数法、ロジスティック回帰のパラメータ推定を行いました。
GMMの利点はモーメント条件さえわかればどのようなモデルも統一的に扱うことができるところです。実務上の制約を考慮したモデルを推定したいことはよくあるので、それらの制約をモーメント条件として表現できれば役立つことがありそうです。また、適切な操作変数を見つけるのが大変なので実務で使える場面は限られますが、既存の回帰モデルに操作変数法を組み合わせて使うのも容易です。
欠点は非線形モデルのGMMの初期値依存性が強いところです。実データで使う場合は最尤推定など他の推定方法と比較しないと使い物にならないぐらい推定が不安定でした。データやモデルに依存するとは思いますが、モーメント条件はなるべく線形にしたほうがよさそうです。他にも推定値の分散が最尤法より大きいなどといった欠点もありますがサンプルデータで確認した範囲ではよほどデータ数が少なくない限り問題にはならないようでした。
