概要
前回の記事では, 回帰式モデルの説明変数に未観測な変数があった場合に発生する欠落変数バイアスを対処する手法として, 操作変数法(Instrument variables; IV法)と２段階最小二乗法(two stage least squares; 2SLS ) を紹介した.
https://zenn.dev/akitek/articles/72acbf42ac804a
今回は, それらの一般形であるGMM(generalized moment method; 一般化積率法, 一般化モーメント法） を考えていく.
（一般化ではない）モーメント法とは
まず, モーメント法とはなにかをとりあげよう. （私も初めて聞いた名前であるため, 各種リソースから説明を引用する）
モーメント法は, モーメント条件（直交条件とも呼ばれる）を用いることにより推定を行う方法. モーメント条件とはモデルより導かれる, 確率変数とパラメーターのある関数の期待値が0となるという条件のことである.
あるいは,
モーメント法とは, モデルが満たすべきモーメント条件を, 標本モーメント条件で置き換えることでパラメータを推定する方法である.
...とのことだが, いまいちピンとこない.もう少し具体的な説明が欲しいので, 以下の記事を参考に, （まずは一般化されていない）モーメント法のやりたいことを考える.
https://zenn.dev/nigimitama/articles/2024-04-23-mme-and-mle
いま, 確率変数Xと, k次元のパラメータ\theta = (\theta_1, \theta_2, ..., \theta_k)があるとすると, 次のモーメント条件を満たすようなk個のスコア関数(m_k(\cdot, \cdot))を考えて, 標本モーメント(例：モーメントの平均）に置き換えて, パラメータ\thetaを解く, という方法らしい.
E[m(X, \theta)] =E
\begin{pmatrix}
m_1(X, \theta) \\
\vdots \\
m_k(X, \theta) \\
\end{pmatrix} = 0
例：正規分布の場合
具体例で, モーメント法の流れを確認する. 例えば, データX_iが, 未知の平均\muと分散\sigma^2をパラメータにもつ正規分布N(\mu, \sigma^2)から得られたと仮定する.
この正規分布のパラメータをモーメント法で推定したい.ここで, \theta = [\mu,\sigma^2]^Tとする.
ここで, 例えば二つのスコア関数として,
m_1(X,\theta) = X_i - \mu, m_2(X,\theta) = (X_i - \mu)^2 - \sigma^2
を考える. そして, 二つのスコア関数の期待値の計算を求めると
\begin{align} E[m_1(X,\theta)] &= E[X_i - \mu] = E[X_i] - \mu = 0 \\
E[m_2(X,\theta)] &= E[(X_i - \mu)^2 - \sigma^2] = E[(X_i - \mu)^2] - \sigma^2 = 0 \end{align}
となる. このE[m_i(X,\theta)] = 0がモーメント条件である. そして, この期待値計算を, 標本によって類似させた標本モーメント条件を解くことで, パラメータの推定量を得ることができる.
よって,
\begin{align} E[m_1(X,\theta)] &\simeq \frac{1}{N}\sum_{k=1}^{N}m_1(X_k, \theta) = \frac{1}{N}\sum_{k=1}^{N} X_i - \mu = 0 \\
E[m_2(X,\theta)] &\simeq \frac{1}{N}\sum_{k=1}^{N}m_2(X_k, \theta) = \frac{1}{N}\sum_{k=1}^{N}(X_i - \mu)^2 - \sigma^2 = 0 \end{align}
この等式を解けば, 正規分布の平均\muと分散\sigma^2の推定量は,
\begin{align}
\hat{\mu} &= \frac{1}{N}\sum_{k=1}^{N} X_i \\
\hat{\sigma^2} &= \frac{1}{N}\sum_{k=1}^{N} (X_i - \mu)^2
\end{align}
最小二乗法や最尤法などでよく見慣れた推定量が, モーメント法からも算出できている
例：線形回帰モデル
次に, 回帰モデルでも同様に考える. 今, 観測データY_iが, k個の説明変数X = (X_1, X_2, ..., X_kと, その回帰係数\beta = (\beta_1, \beta_2, ..., \beta_k)および誤差項\varepsilon_i ~ N(0, \sigma^2)の和で求められる回帰式を考える.
Y_i = X_i^{\top}\beta + \varepsilon_i
ここで, 説明変数X_iと, 誤差項\varepsilon_iが互いに無相関であるとする E[X_i \varepsilon_i] = 0
正規分布の例と同様に, パラメータ\betaをモーメント法で推定する. 回帰モデルでは, 説明変数と誤差項とが無相関であることをそのままモーメント条件として扱うと,
E[X_i\varepsilon_i] = E[X_i(Y_i - X_i^{\top}\beta = 0
という等式がk個作ることができる.(回帰係数\betaの数だけ)
そして, 標本モーメント条件（標本平均）で, この等式を置き換えた式,
\frac{1}{n}\sum_{i=1}^{n}X_i\varepsilon_i = \frac{1}{n}\sum_{i=1}^{n})X_i(Y_i - X_i^{\top}\beta) = 0
を, 行列式 で表すと
\frac{1}{n}X^{\top}\varepsilon = \frac{1}{n} X^{\top}(Y - X\beta) = 0
後ろの等式を解くと,
\begin{align}
\frac{1}{n} X^{\top}Y &= \frac{1}{n} X^{\top}X\beta) \\
\beta &= (X^{\top}X)^{-1}X^{\top}Y
\end{align}
と, よく知られた最小二乗推定量が得られる.
モーメント法の限界
正規分布モデルや線形回帰モデルに対し, モーメント法をつかっても, しっかりとパラメータの推定量が計算できることがわかった.
しかし, モーメント法にはある限界点が存在する. それは, モーメント条件の数J（つまり、スコア関数の数）が, 未知のパラメータの数Kと一致していないといけないのだ.
例えば, 正規分布に対しては, 他にも, 4次のモーメントに対応する m_3(X,\theta)= (X_i - \mu^4) - 3\sigma^3 もひとつのスコア関数, つまりは, この期待値を0とおいたモーメント条件がつくれてしまう. モーメント条件を三つ用意して, パラメータ\mu, \sigma^2を推定しようとすると, J(=3) > K(=2)の関係から, 連立不等式を解くことができない.
この一例にあるように, モーメント条件の数Jと, 未知のパラメータの数Kの大小関係はとても重要である.
-
K > Jのとき：未知パラメータの数に対し, 条件式が足りないため, パラメータを一意に求められない(過少識別; under-identified)
-
K = Jのとき：未知パラメータの数と, 条件式の数が一致する(丁度識別; just or exact identified)場合に限り, モーメント法が有効
-
K < Jのとき:未知パラメータの数に対し, 条件式が多いため(過剰識別; over identified), 条件式を同時に満たすパラメータはそもそも存在しないため求まらない
特に,(3)の場合, 通常のモーメント法では, 条件式を, 未知パラメータの数に合うように選別することで対応している
でも, できるならば,すべてのモーメント条件を使用して推定の効率をあげたい という思いに応えたのが, GMMというものらしい.
一般化モーメント法
話を一般化するために, モーメント条件を一般的な関数に置き換える
h(y_i, X_i, \beta) = 0_{[J\times1]}
ここで, h(y_i, X_i, \beta)は, J\times1のベクトル値関数である(つまり, モーメント条件がJ個あるとしている）.
そして, 標本平均をとった標本モーメント条件は以下の通りとなる.
h_s(y,X,\beta)=\frac{1}{n}\sum_{i}^{n}h_s(y_i, X_i, \beta)
ここまでは, いままでの話を一般化しただけであるが, GMMでは次の関数を最小化するパラメータを求めることを考える.これによって得られるパラメータの推定量がGMM推定量
\hat{\beta}_{gmm} = \min_{\beta}{h_s(y,X,\beta)^{\top}V h_s(y,X,\beta)}
上記式のなかにあるVは, 各モーメント条件に対する重みづけであり, Hansen(1982)によれば, 以下の式をVの値として採用することで, GMM推定量の分散が最小化されることを示している.
V=[\frac{1}{n}\sum_{i=1}^{n}h(y_i, X_i, \beta)h(y_i, X_i, \beta)^{\top}]^{-1}
このVを具体的に求める場合, Vの式に, 推定したいパラメータ\beta自身が含まれているので, 計算には工夫が必要である. 『空間統計学 - 自然科学から人文・社会科学まで』(朝倉書店）で紹介されている方法として,
1.適当な初期値となる重み（単位行列など）をVの値として, 先にGMM推定りょう\hat{\beta}_{gmm}^{(0)}を計算(第一段階目）
2. この暫定的な推定値\hat{\beta}_{gmm}^{(0)}を用いて, 今度はVの推定量V = S^{-1}を計算
3. そして最後に,改めて, このS^{^-1}をウェイト行列にもちいて, GMM推定量をもとめる
この推定方法は, optimal GMM, efficiency GMM two-step GMM などと呼ばれる.
GMM推定量は一致性・漸近正規性を持ち, そして, GMM推定量は, モーメント条件のみを使って得られる漸近正規性を満たす推定量の中で分散が最小である
GMMの計算例： 丁度識別 のとき
それでは,GMMの流れを見ていく. 最初は, K = Jのとき, つまり, 未知パラメータの数と, 条件式の数が一致する(丁度識別; just or exact identified)場合を考える.
この場合は, ウェイト行列Vが解析的に求められる.
例として, 標準的な回帰モデルを再び考える
\begin{align}
y_t &= \beta_0 + \beta_1 x_{1t} + ... + \beta_K x_{Kt} + \varepsilon_t \\
&= x_t^{\top}\beta + \varepsilon_t
\end{align}
ここで, x_t = [1, x_{1t}, ..., x_{Kt}]^{\top}, \beta = [\beta_0, \beta_1, ..., \beta_K]^{\top} であり, 未知数はK+1(\betaの数だけ）である.
スコア関数を, h(y_t, x_t, \beta) = x_t \varepsilon_t = [\varepsilon_t, x_{1t}\varepsilon_t,...,x_{Kt}\varepsilon_{t}]^{\top} とすると, モーメント条件は,
E[h(y_t, x_t, \beta]E[x_t \varepsilon_t] = [E(\varepsilon_t),E(x_{1t}, \varepsilon_t), ..., E(x_{Kt}\varepsilon_t)]^{\top} = 0
であり, ちょうどJ = K+1個ある. つまりは, 未知パラメータの数と, 条件式の数が一致する.
そして, 期待値を標本平均で置き換えると,
\begin{align}
h(y_t,x_t,\beta) &= \frac{1}{T} \sum_{t=1}^{T}x_t(y_t - x_t^{\top} \beta) \\
&= \frac{1}{T} \sum_{t=1}^{T}x_t y_t - (\frac{1}{T}\sum_{t=1}^{T} x_t x_t^{\top} )\beta \\
&= s_{xy} - S_{xx}\beta \\
\end{align}
とする. ここで, s_{xy} = \frac{1}{T}\sum_{t=1}^{T}x_t y_t ,S_{xx} = \frac{1}{T}\sum_{t=1}^{T}x_t x_t^{\top} とする.
GMMの目的関数において, ウェイト行列Vを, 単位行列I(つまり, 全ての条件を均等に重みづけ）とおくと,
\begin{align}
h(y_t, x_t, \beta)^{\top} I_{K} h(y_t, x_t, \beta) &= [s_{xy} - S_{xx}\beta]^{\top}[s_{xy} - S_{xx}\beta] \\
&= [s_{xy}^{\top} - \beta^{\top} S_{xx}^{\top}][s_{xy} - S_{xx}\beta] \\
&= s_{xy}^{\top}s_{xy} - s_{xy}^{\top}S_{xx}\beta - \beta^{\top}S_{xx}^{\top}s_{xy} + \beta^{\top}S_{xx}^{\top}S_{xx}\beta \\
\end{align}
最後の式を両辺\betaで微分すると,
-2S_{xx}^{\top}s_{xy} + 2S_{xx}^{\top}S_{xx}\beta = 0 \\
\Leftrightarrow S_{xx}^{\top}S_{xx}\beta = S_{xx}^{\top}s_{xy} \\
\begin{align}
\Leftrightarrow \hat{\beta}_{OLS} &= (S_{xx}^{\top}S_{xx})^{-1}S_{xx}^{\top}s_{xy} \\
&= S_{xx}^{-1}(S_{xx}^{\top})^{-1}S_{xx}^{\top}s_{xy} \\
&= S_{xx}^{-1}s_{xy} \\
&= [\frac{1}{T}\sum_{t=1}^{T}x_t x_t^{\top}]^{-1}\frac{1}{T}\sum_{t=1}^{T}x_t y_t
\end{align}
となり, よくある最小二乗推定量と一致した.
GMMの計算例： 過剰識別 のとき
本題で, モデルが過剰識別の場合, つまり, 未知パラメータの数に対し, モーメント条件式が多い場合を考える.
先ほどと同様に, 線形モデル
\begin{align}
y_t &= \beta_0 + \beta_1 x_{1t} + ... + \beta_K x_{Kt} + \varepsilon_t \\
&= x_t^{\top}\beta + \varepsilon_t
\end{align}
ここで, x_t = [1, x_{1t}, ..., x_{Kt}]^{\top}, \beta = [\beta_0, \beta_1, ..., \beta_K]^{\top} であり, 未知数はK+1(\betaの数だけ）である.
しかし, モデルの内生性の問題により, 操作変数z_tを使って,パラメータ\betaを推定するとする.ただし操作変数の数(z_tの次元)のJ は J > K+1であるとする.
さきほどと同じく, モーメント条件と期待平均を考えると,
\begin{align}
E[z_t \varepsilon_t] &= E[z_t(y_t - x_t^{\top}\beta)] = 0 \\
h(y_t, x_t, z_t, \beta) &= \frac{1}{T}\sum_{t=1}^{T}z_t(y_t - x_t^{\top}\beta) = 0 \\
&= \frac{1}{T}\sum_{t=1}^{T}z_t y_t - \frac{1}{T}\sum_{t=1}^{T}z_tx_t^{\top}\beta) \\
&= s_{xy} - S_{zx}\beta \\
\end{align}
ここで, s_{zy} = \frac{1}{T}\sum_{t=1}^{T}z_t y_t, S_{zx} =\frac{1}{T}\sum_{t=1}^{T}z_tx_t^{\top} とおいた.
一般には上記の方程式を満たす\betaは必ずしも存在しない.
この式を, 目的関数に代入して,
\begin{align}
h(y_t, x_t, z_t, \beta)^{\top} V h(y_t, x_t, z_t, \beta) &= [s_{xy} - S_{zx}\beta]^{\top} V [s_{xy} - S_{zx}\beta] \\
&= [s_{xy}^{\top} - \beta^{\top} S_{zx}^{\top}] V [s_{xy} - S_{zx}\beta] \\
&= s_{xy}^{\top}Vs_{xy} - s_{xy}^{\top}VS_{zx}\beta - \beta^{\top} S_{zx}^{\top}V s_{xy} + \beta^{\top} S_{zx}^{\top} V S_{zx}\beta \\
\end{align}
よって,両辺をパラメータ\betaで微分すると,
-S_{zx}^{\top}V s_{xy} - S_{zx}^{\top}V s_{xy} + 2S_{zx}^{\top} V S_{zx}\beta = 0 \\
\Leftrightarrow S_{zx}^{\top} V S_{zx}\beta = S_{zx}^{\top}V s_{xy} \\
\Leftrightarrow \hat{\beta}_{gmm} = (S_{zx}^{\top} V S_{zx})^{-1}S_{zx}^{\top}V s_{xy}
この推定量の式で, V = Iかつ, z_t = x_t(つまり, 全ての変数が外生変数（操作変数）とみなせる場合には, 最小二乗推定量(OLS推定量）と一致することも確認できる.
２段階最小二乗法との比較
冒頭, GMMは２段階最小二乗法の一般化と述べた. そこで, 本記事の最後として, ２段階最小二乗法(2SLS)とGMMの関連性について確認する.
まず２段階最小二乗法およびGMMの適用例として, 以下の回帰モデルを考える.
y = X\beta + \dot{X}\dot{\beta} + \varepsilon
ここで, Xは, 定数項と外生変数からなるn \times kの説明変数行列(つまり, 外生変数の数はk個）, \dot{X}は, 内生変数からなるn \times lの説明変数行列(つまり, 内生変数の数はl個）. \beta, \dot{\beta}はそれぞれ, 外生変数X,内生変数\dot{X}に対応する回帰係数である.
ここで, 外生変数, 内生変数の違いについてだが, 外生変数は, 内生変数に影響を与える因果関係にあるもので, 自身は他の変数からの影響を受けない変数のことであり, 内生変数は, 外生変数や他の内生変数からの因果的な影響を受けている変数である.
つまり, 内生変数の場合は, 誤差項との相関をもってしまう(C(\dot{X}, \varepsilon)).
さて, ２段階最小二乗法, そして操作変数法は, 操作変数を導入することで, この問題を対処しようとするものだった. ここで操作変数Zは, 内生変数\dot{X}とは相関をもつが, 誤差項varepsilonとは相関を持たないものを選ぶ.
２段階最小二乗法では, ２段階の推定処理を行う.
１段階では, 内生変数と外生変数の集まりR = [X ; \dot{X}]を, 誤差項と相関のない変数群, つまり, 外生変数と操作変数の集まりS = [X ; Z]で推定する.
2段階目では, 従属変数yを推定した\hat{R}に回帰する.
処理の役割として, １段階目で, Rと誤差項\varepsilonの相関を取り除き, 誤差項と無相関な成分\hat{R}を用いて分析を行うというアイデアである.
推定の結果については, \hat{\ddot{\beta}}_{2sls} = (\hat{R}^{\top}\hat{R})^{-1}\hat{R}^{\top}y となり, 最小二乗推定量に類似した表現の推定量を得られる. ここで, \hat{R} = S (S^{\top}S)^{-1}S^{\top}Rである.
そして, 今度はGMMで考える(問題設定は同じである）
モーメント条件は, 誤差項\varepsilonと無相関である, n\times (k+p)の行列S = [X ; Z]をもちいて, 以下の式となる.
E[S^{\top}\varepsilon] = 0_{[(k+p) \times 1]}
ここで, 変数行列R = [X ; \dot{X}], 誤差項行列\ddot{\beta} = [\beta ; \dot{\beta}]とおき,標本モーメント条件（標本平均）をつくると
\frac{S^{\top}(y - R\ddot{\beta})}{n}
である. これを, スコア関数h(y,X,\dot{X}, Z, \ddot{\beta})とおくと, 解くべき目的関すは以下のようにかける.
h(y,X,\dot{X}, Z, \ddot{\beta})^{\top}Vh(y,X,\dot{X}, Z, \ddot{\beta}) \\
(\frac{S^{\top}(y - R\ddot{\beta})}{n})^{\top}V\frac{S^{\top}(y - R\ddot{\beta})}{n}
ここで, ウェイト行列Vを, V=(\frac{\sigma_{\varepsilon}^2 S^{\top}S}{n})として, 先ほどの目的関数を\ddot{\beta}で微分したものを 0とおくとGMM推定量が得られる.
\hat{\ddot{\beta}}_{gmm} = [R^{\top}S(S^{\top}S)^{-1}S^{\top}R]^{-1}R^{\top}S(S^{\top}S)^{-1}S^{\top}y
ここで, 便宜的にP_z = S(S^{\top}S)^{-1}S^{\top}と置き換える. P_zは, べき等行列(P_z P_z = P_z)であるため, P_z^{\top}P_z = P_zが成り立つことを覚えておく.
すると, 上式は,
\hat{\ddot{\beta}}_{gmm} = [R^{\top}P_z R]^{-1}R^{\top}P_z y
２段階最小二乗法においては, \hat{R} = P_z Rとかける. これを, \hat{\ddot{\beta}}_{2sls} = (\hat{R}^{\top}\hat{R})^{-1}\hat{R}^{\top}yに代入すると,
\begin{align}
\hat{\ddot{\beta}}_{2sls} &= (\hat{R}^{\top}\hat{R})^{-1}\hat{R}^{\top}y \\
&= ((P_z R)^{\top}P_z R)^{-1}(P_z R)^{\top} y \\
&= (R^{\top}P_z ^{\top} P_z R)^{^1}R^{\top}P_z y \\
&= (R^{\top} P_z R)^{^1}R^{\top}P_z y \\
&= \hat{\ddot{\beta}}_{gmm}
\end{align}
と, GMM推定量と2SLS推定量が一致することがわかった.
より詳細に2SLS推定量との違いについては, 以下の記事を参考にしてほしい.
https://ill-identified.hatenablog.com/entry/2015/02/22/203650
残りの話題
最後に, 実用的な観点で, 関連する内容をいくつか紹介しておく
Efficient GMM推定量 と J統計量
ウェイト行列Vは, GMM推定量の推定効率（漸近分散） に大きな影響をもたらす. そして, 推定効率をもっともよくするウェイト行列Vは, 以下の式で求まるとわかっている.
V = S^{-1} = (E[h(y,x,\beta][h(y,x,\beta]^{\top})^{-1}
残念ながら, この式を直接解くのは難しいため, 何らかの方法で推定値\hat{S}_{T}を用いるのが一般である.
この推定量\hat{S}_{T}を用いたGMM目的関数を最小にするパラメータを, Efficient GMM 推定量と呼ぶ.そして, この推定量を用いたGMM目的関数の値に, 標本数Tをかけたものを, J-統計量とよぶ
カイ二乗分布とJ検定
J 統計量はモデルの全ての仮定が正しいとき, 自由度 K–Lのカイ二乗分布に従う。
つまり, この統計量が, このカイ二乗分布からの標本の値として大きすぎるときには,(モーメント条件などの)モデルの仮定のどれかが満たされていないと考えられるので、改めてモデルを見直す必要が出てくる。
このように J 統計量を用いて、モデルの定式化の誤りがあるかどうかを検定することを J 検定という
モーメント（条件）選択基準
実証分析ではしばしば J 統計量によってモデルの仮定が棄却される.
この理由として, 使用したモーメント条件のうち正しくないものが含まれている可能性がある. そのような場合, どのモーメント条件が正しくてどのモーメント条件が正しくないかを推定するような方法があれば有用であろう. ここではそのような方法のうち Andrews (1999) と Hall,Inoue, Jana, and Shin (2007) によるモーメント(条件)選択基準（Moment Selection Criteria; MSC) と呼ばれる方法がある, らしい.
弱い操作変数の問題
別の問題として, 操作変数zを用いる場合, 元のモデルで使用されている内生変数xとの相関が非常に弱い操作変数を採用してしまうと, 非常に大きな有限標本においても、β の推定量に非常に大きいバイアスがかかり、かつ標準誤差(の推定量)が非常に小さくなる, という問題が生じる.
これを, 弱い操作変数（Weak Instrument) と呼ばれる.
弱い操作変数の存在を検定する仕組みには, Cragg-Donald推定量がある.
まとめ
今回は, 操作変数を導入して欠落変数バイアスに対象する操作変数法,２段階最小二乗法の一般形である一般化モーメント法(GMM)を紹介した. モーメント条件に基づき, その標本平均をもとめ, パラメータを推定しようというものである. 特に, 操作変数の数が, 未知のパラメータ数よりも多い場合にも有効な手法である.
より詳細な手法の説明は専門書にお任せしよう
Hayashi, Fumio (2000) Econometrics: Princeton University Press.
雑記
かなり数式が多くなってしまった. その理由として, 省略されがちな数式展開などを（個人的な趣味で）やっていたためであるが, そのぶん理解が深まったと思う.
ようやく「空間統計手法」の準備ができた...💦
参考サイト
https://zenn.dev/nigimitama/articles/2024-04-23-mme-and-mle
コメント：モーメント法の説明としてまずイメージがしやすい記事だと思う
https://blog.neko-ni-naritai.com/entry/idea-of-gmm#f-2f38e239
GMMのイメージをつかむために多いに活用させていただいた
https://ill-identified.hatenablog.com/entry/2015/02/22/203650#f-07376afa
GMM法の詳細な説明と, ２段階最小二乗法との関連性について詳細に説明している
また, 非線形モデルへのGMMへの適用記事も書かれているのおすすめだ
https://user.keio.ac.jp/~nagakura/zemi/ts5_slide_2017A.pdf
記事内の式の導出で大変参考になった. また, 発展的な内容にも踏み込んでいる
https://www.anarchive-beta.com/entry/2024/08/02/210000
行列の式に関してのテクニックはこの記事を参考にした. べき等行列などの定理に.
Discussion