using Random
using Statistics
using LinearAlgebra
using DataFrames
using CSV
using ProjectRoot
using Optimization
using ForwardDiff
import OptimizationNLopt
import NLopt
using Distributions
using FastGaussQuadrature
using Roots
using PrettyTables
using Plots
using LaTeXStrings4 GMM
Generalized Method of Moments
構造推定 (structural estimation) では, 経済モデルのパラメータ \(\boldsymbol{\theta}\) を, モデルがデータをよく再現するように定めます. その代表的な方法が, モデルの最適化条件 (一階条件) から導かれるモーメント条件を使う一般化モーメント法 (GMM) です. この章では, 前半でモーメント条件を個票の上で直接評価できる場合の GMM と, モーメントがパラメータを定めるかという識別の問題, そして推定量の標準誤差を扱います. 後半では, モデル側の対応物が数値的にしか計算できない場合を扱います. これは最小距離推定 (minimum distance estimation), モデル側をシミュレーションで評価するなら SMM と呼ばれますが, モーメント条件の書き方が違うだけで, やはり GMM の一種です. 家族経済学では, 家計内の配分や交渉力といった直接観察できない対象を推定するためにこれらの手法が広く使われます.
4.1 一階条件による GMM
データ \(\{\boldsymbol{x}_i\}_{i=1}^N\) と, 真のパラメータ \(\boldsymbol{\theta}_0\) のもとで期待値がゼロになるモーメント条件 \(\mathbb{E}[\boldsymbol{g}(\boldsymbol{x}_i, \boldsymbol{\theta}_0)] = \boldsymbol{0}\) が与えられたとします. 典型的には \(\boldsymbol{g}\) は最適化の一階条件で, 複数の条件を束ねてベクトルにできます. その標本版 \(\bar{\boldsymbol{g}}_N(\boldsymbol{\theta}) = \frac{1}{N} \sum_{i} \boldsymbol{g}(\boldsymbol{x}_i, \boldsymbol{\theta})\) を可能な限りゼロに近づける \(\boldsymbol{\theta}\) を選ぶのが GMM です.
定義 4.1 (GMM 推定量) GMM 推定量 \(\hat{\boldsymbol{\theta}}^{\text{GMM}}\) は, 正定値の重み行列 \(W\) に対して, 次のように定義される: \[ \hat{\boldsymbol{\theta}}^{\text{GMM}} := \arg\min_{\boldsymbol{\theta}}\ \bar{\boldsymbol{g}}_N\left(\boldsymbol{\theta}\right)^\top W\, \bar{\boldsymbol{g}}_N\left(\boldsymbol{\theta}\right). \] パラメータ数と同じだけモーメント条件があれば (丁度識別, just-identified), \(W\) に依らず \(\bar{\boldsymbol{g}}_N(\hat{\boldsymbol{\theta}}) = \boldsymbol{0}\) となる. モーメント条件の方が多い場合 (過剰識別, over-identified) は \(W\) が推定効率に影響し, 最適な \(W\) はモーメントの分散共分散行列の逆行列である.
4.1.1 推定量の性質
まず漸近分布です. GMM 推定量は, 観測が定常かつエルゴード的であれば一致性と漸近正規性を持ちます (Hansen 1982). その漸近分散は, ヤコビアン \(G\) とモーメントの分散共分散行列 \(V\) から計算できます.
命題 4.1 (GMM 推定量の漸近分布) モーメント条件 \(\mathbb{E}\left[\boldsymbol{g}\left(\boldsymbol{x}_i, \boldsymbol{\theta}_0\right)\right] = \boldsymbol{0}\) が正しく, \(\boldsymbol{\theta}_0\) が識別され, 標準的な正則条件が成り立つとする. \(G = \partial\, \mathbb{E}\left[\boldsymbol{g}\right] / \partial \boldsymbol{\theta}^\top\) (\(p \times k\)), \(V = \lim_{N \to \infty} \text{Var}\left(\sqrt{N}\, \bar{\boldsymbol{g}}_N\left(\boldsymbol{\theta}_0\right)\right)\) とおくと, \[ \sqrt{N}\left(\hat{\boldsymbol{\theta}} - \boldsymbol{\theta}_0\right) \xrightarrow{\ d\ } \mathcal{N}\left(\boldsymbol{0},\ \left(G^\top W G\right)^{-1} G^\top W V W G \left(G^\top W G\right)^{-1}\right). \tag{4.1}\] 特に \(W = V^{-1}\) のとき, 分散は \(\left(G^\top V^{-1} G\right)^{-1}\) に簡約され, これが \(W\) の選び方の中で最小になる (効率的 GMM).
証明は Newey and McFadden (1994) などを参考にしてください. 式 4.1 は, 重み行列 \(W\) の選び方が推定量の効率を左右することを言っています. 最適なのは \(W = V^{-1}\) ですが, ここに循環があります. \(V = \lim_N \operatorname{Var}\left(\sqrt{N}\, \bar{\boldsymbol{g}}_N\left(\boldsymbol{\theta}_0\right)\right)\) は真値のまわりで定義される量なので, \(V\) を推定するには \(\boldsymbol{\theta}\) の推定値が必要です. 一方, 効率的に \(\boldsymbol{\theta}\) を推定するには \(V\) が必要です. この鶏と卵の順番を決めるのが2段階 GMM です.
アルゴリズム 4.1 (2段階 GMM)
- 第1段階. 適当な正定値行列 (基本は \(W = I\)) を使って最小化し, \(\hat{\boldsymbol{\theta}}_1\) を得る. 効率的ではないが一致推定量である.
- 重み行列の推定. モーメント条件が独立同時分布なら, \(\hat{S} = \frac{1}{N} \sum_i \boldsymbol{g}\left(\boldsymbol{x}_i, \hat{\boldsymbol{\theta}}_1\right) \boldsymbol{g}\left(\boldsymbol{x}_i, \hat{\boldsymbol{\theta}}_1\right)^\top\) が \(V\) の一致推定量になる.
- 第2段階. \(W = \hat{S}^{-1}\) で最小化し直す. これが効率的 GMM 推定量で, その漸近分散は \(\frac{1}{N}\left(G^\top V^{-1} G\right)^{-1}\) である.
実務上の注意が2つあります. 第一に, 第1段階の \(W = I\) はモーメントの単位に依存する恣意的な選択です. スケールの大きく違うモーメントを混ぜるときは, 各モーメントを標本標準偏差で割った対角行列を第1段階に使うほうが安定します. 第二に, \(\hat{S}\) は \(p \times p\) の行列を \(N\) 個の観測から推定したものです. モーメントの本数が多く標本が小さいと \(\hat{S}^{-1}\) の推定誤差が大きくなり, 漸近的には効率的なはずの第2段階が, 有限標本ではかえって悪い推定量になることが知られています. 共分散構造の推定でこの小標本バイアスを詳しく調べたのが Altonji and Segal (1996) です. その場合は非対角成分を捨てて対角の重みで済ませる, という保守的な選択がよく採られます. この章の後半 (4.2節) で使うのがまさにそれです.
最後に, 効率的な重み行列にはおまけがついてきます. 最小化された目的関数そのものが検定統計量になるのです. モーメントがパラメータより多いとき, すべてのモーメントを同時にゼロにはできません. その「残り」がモデルと整合的な大きさに収まっているかを見るのが過剰識別の検定です.
命題 4.2 (過剰識別の検定 (Hansen の \(J\))) \(W = \hat{V}^{-1}\) のとき, 最小化された目的関数に \(N\) を掛けたものは, \(p\) 本のモーメント条件がすべて正しいという帰無仮説の下で自由度 \(p - k\) の \(\chi^2\) 分布に従う. \[ J = N\, \bar{\boldsymbol{g}}_N\left(\hat{\boldsymbol{\theta}}\right)^\top \hat{V}^{-1} \bar{\boldsymbol{g}}_N\left(\hat{\boldsymbol{\theta}}\right) \ \xrightarrow{\ d\ }\ \chi^2\left(p - k\right) \]
\(J\) が棄却したときに分かるのは, 「\(p\) 本のモーメント条件を1つの \(\boldsymbol{\theta}\) で同時に満たせない」ということだけで, どれが誤っているのかは教えてくれません. 原因の候補は, モデルの定式化の誤り, モーメント条件そのもの (たとえば操作変数の外生性) の誤り, 標本の偶然, そして標本の大きさの4つです. 最後は見落とされがちですが重要で, 現実のモデルはどれも近似なので, \(N\) が大きくなるほど \(J\) はごく小さな定式化の誤りにも反応し, たいていのモデルは棄却されます.
そこで実務では, 推定値における各モーメントのずれ \(\bar{g}_j\left(\hat{\boldsymbol{\theta}}\right)\) を標準誤差で割って並べ, どの条件が外れているのかを見ます. 外れたモーメントが特定できれば, 次に疑うのはその式を支えている経済的な仮定です. そのうえで, モデルを一般化するか, 問題のモーメントを推定から落とすか, 棄却を報告したうえで結果を示すかを選ぶことになります. 2番目は手軽に見えますが, 落とすと識別が壊れることがあるので, どのモーメントがどのパラメータを支えているかを知らないままでは選べません. 逆に, 検定を通ったからモデルが正しいということにもなりません. 識別が弱ければ検出力もないからです. Andrews et al. (2017) が推定値の感応度を報告することを勧めるのも, 検定の合否だけでモデルを評価する慣行への処方箋でした. 検定は診断の入口であって結論ではない, というくらいの距離感が実務的です.
4.1.2 モデル
GMM の真価は, 閉じた形の回帰式に落ちない複数の非線形な最適化条件をまとめてモーメント条件にできる点にあります. Lise and Yamada (2019) は, 夫と妻の個人別消費・時間配分のパネルデータを使い, 家計内の資源配分を特徴づける一階条件の体系を非線形 GMM で推定しています. ここでは, そのモデルの最小限のトイモデルを作ります.
家計は妻 \(W\) と夫 \(H\) からなります. メンバー \(j \in \{W, H\}\) は, 自分の私的消費 \(c_j\) と自分の余暇 \(\ell_j\) から効用を得ます. 効用関数は, \(c_j\) と \(\ell_j\) の CES 集計を CRRA で包んだ形とします. \[ u_j\left(c_j, \ell_j\right) = \frac{1}{1 - \sigma_j}\left[\alpha_j c_j^{\phi_j} + \left(1 - \alpha_j\right) \ell_j^{\phi_j}\right]^{\frac{1 - \sigma_j}{\phi_j}}. \tag{4.2}\] \(\alpha_j\) は消費と余暇の相対的な選好, \(\phi_j\) は両者の代替の弾力性 (\(=1/(1-\phi_j)\)) を決めます. ここでは, 妻と夫の選好に異質性があることが許容されています (\(\alpha_W \neq \alpha_H\), \(\phi_W \neq \phi_H\)).
家計は効率的に資源を配分する (collective model) とし, 妻への Pareto ウェイトを \(\mu \in (0, 1)\) とします. 資産 \(a\) を状態変数として, 家計の動的な問題は次の価値関数で書けます. \[ V(a) = \max_{c_W, \ell_W, c_H, \ell_H, a'} \mu_W u_W\left(c_W, \ell_W\right) + \mu_H u_H\left(c_H, \ell_H\right) + \delta\, \mathbb{E} V(a'), \tag{4.3}\] subject to \[ c_W + c_H + w_W \ell_W + w_H \ell_H + a' = w_W + w_H + \left(1 + r\right) a =: y. \] 各メンバーの時間の賦存量を1に基準化し, 余暇を自分の賃金 \(w_j\) で評価しています. 右辺の \(y\) は full income です. また, \(\mu_W, \mu_H \in (0, 1)\) は, 家計内の交渉力を表す Pareto ウェイトで, \(\mu_W + \mu_H = 1\) です.
モーメント条件
予算制約のラグランジュ乗数を \(\lambda\) とすると, 内点解での一階条件は次のようになります. \[ \mu \frac{\partial u_j}{\partial c_j} = \lambda, \quad \mu \frac{\partial u_j}{\partial \ell_j} = \lambda w_j, \quad j \in \{W, H\}, \tag{4.4}\]
\[ \lambda = \delta\, \mathbb{E} V'(a') \qquad (\text{Euler equation}). \] ここで \(\lambda\) は 富の限界価値 (marginal value of wealth) で, 家計の効率的な配分を特徴づける共通の尺度です. すべての一階条件は \(\lambda\) を含むので, 式 4.4 の比を取ると \(\lambda\) が消え, 家計内の配分の条件が残ります.
(1) 各メンバーの期内 MRS 条件. \(j \in \{W, H\}\) について, 消費と余暇の比を取ると \(\lambda\) も \(\mu\) も消えます. \[ \frac{\alpha_j}{1 - \alpha_j}\left(\frac{c_j}{\ell_j}\right)^{\phi_j - 1} = \frac{1}{w_j}. \tag{4.5}\]
(2) 家計内の消費配分条件. 妻と夫の消費の一階条件の比を取ると, \(\lambda\) が消えて \(\mu\) が残ります. 余暇についての同様の条件は 式 4.5 と 式 4.6 から導かれるので, 独立な条件はこの3本です.
Noteなぜオイラー方程式を使わないのか
式 4.3 は動的な問題なので, オイラー方程式 \(\lambda = \delta \mathbb{E}V'(a')\) も一階条件の一部です. しかし Lise and Yamada (2019) は, これを推定にあえて使いません. 理由は2つあります.
第一に, 頑健性です. 流動性制約に直面する家計ではオイラー方程式は成立しませんが, 期内条件 式 4.5・式 4.6 は制約の有無に関わらず成立します. 効率的な配分でありさえすればよいからです.
第二に, 識別です. 異時点間の代替の弾力性 (すなわち \(\sigma_j\)) を信頼できる形で識別するには, 家計が直面する金利の変動を使う必要があります. しかし家計パネルには通常そうした情報がありません. そこで Lise and Yamada (2019) は先行研究にならい \(\sigma_W = \sigma_H = 1.5\) に固定します. 以下でも同じく \(\sigma_j = 1.5\) を既知として扱い, 残りのパラメータ \(\left(\phi_W, \phi_H, \alpha_W, \alpha_H, \mu\right)\) を推定します.
式 4.5 と 式 4.6 の対数を取り, \(A_j := \alpha_j c_j^{\phi_j} + (1 - \alpha_j)\ell_j^{\phi_j}\) とおくと, \[ \begin{aligned} e_{ji} &= \log \frac{\alpha_j}{1 - \alpha_j} + \left(\phi_j - 1\right) \log \frac{c_{ji}}{\ell_{ji}} + \log w_{ji} & j \in \{W, H\},\\ e_{3i} &= \log \frac{\mu}{1 - \mu} + \log \frac{\alpha_W}{\alpha_H} + \frac{1 - \sigma_W - \phi_W}{\phi_W} \log A_{W_i} - \frac{1 - \sigma_H - \phi_H}{\phi_H} \log A_{H_i} \\ &+ \left(\phi_W - 1\right) \log c_{W_i} - \left(\phi_H - 1\right) \log c_{H_i}. \end{aligned} \] モデルが正しければ, これらの残差は平均ゼロになります. つまり, \(\boldsymbol{e}_i := (e_{Wi}, e_{Hi}, e_{3i})^\top\) に対して, \(\mathbb{E}\left[\boldsymbol{e}_i\right] = \boldsymbol{0}\) です. ただし, これだけではパラメータ5個に対して3本のモーメントしかなく, 識別できません. そこで, 残差と賃金の直交条件を加えます. 賃金は外生なので, モデルが正しければ残差と相関しないはずです.
ただし, どの賃金でもよいわけではありません. 賃金にも測定誤差が乗るからです. 観測される賃金を \(w_{ji} = w_{ji}^{*} e^{\varepsilon_{ji}^{w}}\) とすると, 式 4.5 の残差には \(\log w_{ji}\) がそのまま入っているので, \(e_{ji}\) は測定誤差 \(\varepsilon_{ji}^{w}\) を含みます. 自分の賃金を掛けた条件 \(\mathbb{E}\left[e_{ji} \log w_{ji}\right]\) は, 残差と操作変数の両方に同じ誤差が入るぶんだけゼロから外れます. 一方, 配偶者の賃金なら測定誤差は独立なので直交条件は成り立ちます. 残差 \(e_{3i}\) には賃金が入らないので, こちらは両方の賃金を使えます. 使えるモーメントは7本です. \[ \begin{aligned} \mathbb{E}\left[e_{Wi}\left(1, \log w_{Hi}\right)^{\top}\right] &= \boldsymbol{0},\\ \mathbb{E}\left[e_{Hi}\left(1, \log w_{Wi}\right)^{\top}\right] &= \boldsymbol{0},\\ \mathbb{E}\left[e_{3i}\left(1, \log w_{Wi}, \log w_{Hi}\right)^{\top}\right] &= \boldsymbol{0}. \end{aligned} \tag{4.7}\] 自分の賃金を捨てるぶん, 曲率 \(\phi_j\) の情報は減ります. パネルデータがあれば, 測定誤差が時点間で独立であることを使って1年前の賃金を操作変数にでき, 情報を取り戻せます. Lise and Yamada (2019) が実際に採っているのがこの戦略です.
データ
真値を Lise and Yamada (2019) の推定値に近い水準に取り, \(N = 4000\) 世帯の人工的なクロスセクションデータを作ります. \[ \left(\phi_W, \phi_H, \alpha_W, \alpha_H, \mu\right) = \left(0.2,\ 0.6,\ 0.45,\ 0.55,\ 0.45\right). \]
データ生成過程 (DGP) は次の通りです. まず, 各家計 \(i\) の賃金を対数正規分布から独立に引き, フルインカムを定めます. \[ \begin{aligned} \log w_{Wi} &\sim \mathcal{N}\left(-0.15,\ 0.15^2\right),\\ \log w_{Hi} &\sim \mathcal{N}\left(0,\ 0.15^2\right), \\ y_i &= 1.05 \left(w_{Wi} + w_{Hi}\right). \end{aligned} \] 次に, 真の賃金とフルインカムを所与として, 真の配分 \(\left(c_{ji}^*, \ell_{ji}^*\right)\) を一階条件 式 4.5・式 4.6 と予算制約から閉じた形で解きます. 最後に, 消費・余暇・賃金に乗法的な測定誤差を乗せたものを「観測値」とします. \[ c_{ji} = c_{ji}^* e^{\varepsilon_{ji}^c}, \quad \ell_{ji} = \ell_{ji}^* e^{\varepsilon_{ji}^\ell}, \quad w_{ji} = w_{ji}^* e^{\varepsilon_{ji}^w}, \quad \varepsilon_{ji}^c,\ \varepsilon_{ji}^\ell \overset{\text{iid}}{\sim} \mathcal{N}\left(0,\ 0.03^2\right), \quad \varepsilon_{ji}^w \overset{\text{iid}}{\sim} \mathcal{N}\left(0,\ 0.01^2\right). \] 賃金の測定誤差は消費や余暇より小さく取っていますが, 大きさは問題ではありません. どれだけ小さくても, 自分の賃金は操作変数として無効になります.
Code
σ = 1.5 # fixed: unidentified without the Euler equation
φW0, φH0, αW0, αH0, μ0 = 0.2, 0.6, 0.45, 0.55, 0.45
function alloc(wW, wH, y, φW, φH, αW, αH, μ)
# c/ℓ from the own-MRS condition
c_over_ℓ_W = ((1 - αW) / (αW * wW))^(1 / (φW - 1))
c_over_ℓ_H = ((1 - αH) / (αH * wH))^(1 / (φH - 1))
# A_j / ℓ_j^φ_j, i.e. A_j at ℓ_j = 1
A_unit_W = αW * c_over_ℓ_W^φW + (1 - αW)
A_unit_H = αH * c_over_ℓ_H^φH + (1 - αH)
# leisure up to a common scale, from the sharing condition
ℓrel_W = (μ * A_unit_W^((1 - σ - φW) / φW) * αW *
c_over_ℓ_W^(φW - 1))^(1 / σ)
ℓrel_H = ((1 - μ) * A_unit_H^((1 - σ - φH) / φH) * αH *
c_over_ℓ_H^(φH - 1))^(1 / σ)
# full-income cost of that bundle; the budget pins down the scale
cost = (c_over_ℓ_W + wW) * ℓrel_W + (c_over_ℓ_H + wH) * ℓrel_H
ℓW = ℓrel_W * y / cost
ℓH = ℓrel_H * y / cost
return (c_over_ℓ_W * ℓW, c_over_ℓ_H * ℓH, ℓW, ℓH) # (cW, cH, ℓW, ℓH)
end
Random.seed!(9)
N = 4000
ws_W★ = exp.(-0.15 .+ 0.15 .* randn(N)) # true wages
ws_H★ = exp.(0.15 .* randn(N))
ys = (ws_W★ .+ ws_H★) .* 1.05 # full income
allocations = [alloc(ws_W★[i], ws_H★[i], ys[i], φW0, φH0, αW0, αH0, μ0)
for i in 1:N]
σ_ε_cℓ, σ_ε_w = 0.03, 0.01 # measurement error
cs_W = getindex.(allocations, 1) .* exp.(σ_ε_cℓ .* randn(N))
ℓs_W = getindex.(allocations, 3) .* exp.(σ_ε_cℓ .* randn(N))
cs_H = getindex.(allocations, 2) .* exp.(σ_ε_cℓ .* randn(N))
ℓs_H = getindex.(allocations, 4) .* exp.(σ_ε_cℓ .* randn(N))
ws_W = ws_W★ .* exp.(σ_ε_w .* randn(N)) # observed wages
ws_H = ws_H★ .* exp.(σ_ε_w .* randn(N))
households = DataFrame(; id=1:N, wW=ws_W, wH=ws_H, cW=cs_W, cH=cs_H,
lW=ℓs_W, lH=ℓs_H)
println("First 5 of the $(N) simulated households:")
first(households, 5)First 5 of the 4000 simulated households:5×7 DataFrame
| Row | id | wW | wH | cW | cH | lW | lH |
|---|---|---|---|---|---|---|---|
| Int64 | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | |
| 1 | 1 | 0.915595 | 1.08042 | 0.425136 | 0.783264 | 0.596384 | 0.362711 |
| 2 | 2 | 0.99009 | 0.975259 | 0.43097 | 0.694002 | 0.542136 | 0.438244 |
| 3 | 3 | 0.701922 | 0.905238 | 0.332111 | 0.532203 | 0.617971 | 0.421421 |
| 4 | 4 | 1.10593 | 1.22219 | 0.492344 | 0.91595 | 0.587731 | 0.334796 |
| 5 | 5 | 0.87725 | 1.17662 | 0.413801 | 0.769628 | 0.630747 | 0.309514 |
Code
scatter(log.(ws_W), log.(cs_W ./ ℓs_W), ms=1.5, ma=0.25, label="Wife", msw=0)
scatter!(log.(ws_H), log.(cs_H ./ ℓs_H), ms=1.5, ma=0.25, label="Husband",
msw=0)
plot!(xlabel=L"\log w_j", ylabel=L"\log (c_j / \ell_j)")注: The intra-temporal MRS condition (式 4.5) in logs. The slope identifies \(\phi_j\) and is free of the Pareto weight.
図 4.1 から識別の直観が得られます. \(\log \left(c_j / \ell_j\right)\) を \(\log w_j\) に対して描くと直線に乗り, その傾き \(1/\left(1 - \phi_j\right)\) が夫と妻で異なることが, 選好の異質性を直接示しています.
4.1.3 推定
式 4.7 の7本のモーメントでパラメータ5個を推定します. 過剰識別 (\(p = 7 > k = 5\)) なので, 4.1.1節の手順どおり2段階で最小化します.
点推定 残差 \(\boldsymbol{e}_i\) は \(\boldsymbol{\theta}\) の対数とべき乗だけでできた閉形式なので, 目的関数 \(\bar{\boldsymbol{g}}\left(\boldsymbol{\theta}\right)^\top W \bar{\boldsymbol{g}}\left(\boldsymbol{\theta}\right)\) の勾配は解析的に計算できます. て計算で勾配を導出するのは面倒なので, 自動微分を (2.2節) を用いた、 準ニュートン法 (L-BFGS, 1.3節) を使って最小化します. 最適化は Optimization.jl の共通インターフェース (第1章) で行い, OptimizationFunction に AutoForwardDiff() を渡すと, ForwardDiff.jl が目的関数のコードから勾配を計算します. パラメータの制約 \(\phi_j < 1\), \(\alpha_j, \mu \in \left(0, 1\right)\) は, 変数変換ではなく箱制約 lb / ub としてソルバーに渡します.
instruments_W = hcat(ones(N), log.(ws_H)) # wife's equation: husband's wage
instruments_H = hcat(ones(N), log.(ws_W)) # husband's equation: wife's wage
instruments_3 = hcat(ones(N), log.(ws_W), log.(ws_H)) # sharing: both wages
function resid(θs)
φW, φH, αW, αH, μ = θs
As_W = αW .* cs_W .^ φW .+ (1 - αW) .* ℓs_W .^ φW
As_H = αH .* cs_H .^ φH .+ (1 - αH) .* ℓs_H .^ φH
es_W = log(αW / (1 - αW)) .+ (φW - 1) .* log.(cs_W ./ ℓs_W) .+ log.(ws_W)
es_H = log(αH / (1 - αH)) .+ (φH - 1) .* log.(cs_H ./ ℓs_H) .+ log.(ws_H)
es_3 = log(μ / (1 - μ)) + log(αW / αH) .+
((1 - σ - φW) / φW) .* log.(As_W) .-
((1 - σ - φH) / φH) .* log.(As_H) .+
(φW - 1) .* log.(cs_W) .- (φH - 1) .* log.(cs_H)
return (es_W, es_H, es_3)
end
# stack the seven moments into an N×7 matrix
function moments(θs)
es_W, es_H, es_3 = resid(θs)
return hcat(es_W .* instruments_W, es_H .* instruments_H,
es_3 .* instruments_3)
end
gbar(θs) = vec(mean(moments(θs), dims=1))
lb = [-5.0, -5.0, 0.01, 0.01, 0.01] # φj < 1, αj, μ ∈ (0, 1)
ub = [0.99, 0.99, 0.99, 0.99, 0.99]
# quasi-Newton with box constraints; the gradient comes from ForwardDiff
minimize(f, θs0) = solve(
OptimizationProblem(
OptimizationFunction((θs, _) -> f(θs), Optimization.AutoForwardDiff()),
θs0; lb, ub),
NLopt.LD_LBFGS(), maxiters=20_000)
θs0 = [0.3, 0.3, 0.5, 0.5, 0.5]
# step 1: W = I
solution_step1 = minimize(θs -> (gs=gbar(θs); gs' * I(7) * gs), θs0)
g_obs_1 = moments(solution_step1.u) # row i: g(x_i, θ̂_1)'
Ŝ_1 = (g_obs_1' * g_obs_1) ./ N # moment covariance at θ̂_1
# step 2: W = Ŝ_1⁻¹
solution_step2 = minimize(θs -> (gs=gbar(θs); gs' * inv(Ŝ_1) * gs),
solution_step1.u)
θ̂s = solution_step2.u
param_labels_gmm = [raw"$\phi_W$", raw"$\phi_H$", raw"$\alpha_W$",
raw"$\alpha_H$", raw"$\mu$"]
θs_true = [φW0, φH0, αW0, αH0, μ0]標準誤差 二段階 GMM を用いているため, 分散は \(\left(\hat{G}^\top \hat{S}^{-1} \hat{G}\right)^{-1} / N\) に簡約化されます. ヤコビアン \(\hat{G} = \partial \bar{\boldsymbol{g}} / \partial \boldsymbol{\theta}^\top\) は, \(\bar{\boldsymbol{g}}\) を計算する関数 gbar がすでにあるので, 自動微分で求められます. 分散共分散行列は \(\hat{S} = \frac{1}{N}\sum_i \boldsymbol{g}\left(\boldsymbol{x}_i, \hat{\boldsymbol{\theta}}\right) \boldsymbol{g}\left(\boldsymbol{x}_i, \hat{\boldsymbol{\theta}}\right)^\top\) です.
g_obs = moments(θ̂s) # row i: g(x_i, θ̂)'
ses_m = vec(std(g_obs, dims=1)) ./ sqrt(N) # standard error of each moment
Ĝ = ForwardDiff.jacobian(gbar, θ̂s) # 7×5, ∂ḡ/∂θ' at θ̂s
Ŝ = (g_obs' * g_obs) ./ N # moment covariance at θ̂s
ses_gmm = sqrt.(diag(inv(Ĝ' * inv(Ŝ) * Ĝ) ./ N))Code
pretty_table(hcat(param_labels_gmm, round.(θ̂s, digits=4),
round.(ses_gmm, sigdigits=2), θs_true);
column_labels=["パラメータ", "推定値", "標準誤差", "真値"],
backend=:markdown, allow_markdown_in_cells=true)| パラメータ | 推定値 | 標準誤差 | 真値 |
|---|---|---|---|
| \(\phi_W\) | 0.195 | 0.0071 | 0.2 |
| \(\phi_H\) | 0.5983 | 0.0052 | 0.6 |
| \(\alpha_W\) | 0.4494 | 0.00078 | 0.45 |
| \(\alpha_H\) | 0.5502 | 0.00065 | 0.55 |
| \(\mu\) | 0.45 | 0.00018 | 0.45 |
| 注: データは真値のもとでモデルから生成したもの. | |||
Code
param_plot_labels_gmm = [L"\phi_W", L"\phi_H", L"\alpha_W", L"\alpha_H",
L"\mu"]
panels = []
for k in 1:5
pad = 4 * ses_gmm[k]
x_min = min(θ̂s[k], θs_true[k]) - pad
x_max = max(θ̂s[k], θs_true[k]) + pad
p = scatter([θ̂s[k]], [1.0], xerror=[1.96 * ses_gmm[k]], ms=5, msw=1.2,
c=1, label="", xlims=(x_min, x_max), ylims=(0.5, 1.5), yticks=false,
title=param_plot_labels_gmm[k])
plot!(p, [θs_true[k], θs_true[k]], [0.5, 1.5], ls=:dash, c=:black, label="")
push!(panels, p)
end
plot(panels..., layout=(5, 1), size=(520, 620), left_margin=4Plots.mm)注: Dot: point estimate. Horizontal bar: 95% confidence interval. Dashed line: true value. Each panel has its own scale.
5つのパラメータはどれも真値のごく近くに戻っています. 選好の曲率 \(\left(\phi_W, \phi_H\right)\) と選好の水準 \(\left(\alpha_W, \alpha_H\right)\) は夫婦で異なる値が, そして家計内の交渉力 \(\mu\) も, それぞれ正しく推定されています. 複数の非線形な期内条件をまとめてモーメント条件にすることで, 選好と交渉力が同時に識別されるわけです. これが Lise and Yamada (2019) が実際の家計パネルで行っている構造推定の核となる仕組みです.1
図 4.2 で目を引くのは, 信頼区間の幅がパラメータによって大きく違うことです (パネルごとに横軸の目盛りが違うことに注意してください). \(\phi_W\) の標準誤差だけが飛び抜けて大きく, 推定値に対する比で見ると \(\phi_H\) の4倍, \(\alpha_H\) の30倍あります. なぜ妻の曲率だけが決まりにくいのか. その答えは次項の識別の分析にあります.
4.1.4 識別
ここまでの推定では, 最小化がうまくいけばパラメータが求まる, ということを暗黙の前提にしてきました. しかしその前提には, そもそもモーメントがパラメータを一意に定めるのかという問いが隠れています. これが識別 (identification) の問題です. どれだけ標本を増やしても区別できないパラメータは, どんな推定法でも推定できません.
推定の入力は「パラメータからモーメントへの写像」\(\boldsymbol{m}\left(\boldsymbol{\theta}\right) = \mathbb{E}\left[\boldsymbol{g}\left(\boldsymbol{x}_i, \boldsymbol{\theta}\right)\right]\) でした. 識別はこの写像が反転できるかという問いとして書けます.
定義 4.2 (識別) モーメント写像 \(\boldsymbol{m} \colon \Theta \to \mathbb{R}^p\) (\(p\) 本のモーメント, \(k\) 個のパラメータ) に対して,
- 大域識別: \(\boldsymbol{m}\left(\boldsymbol{\theta}\right) = \boldsymbol{m}\left(\boldsymbol{\theta}_0\right)\) ならば \(\boldsymbol{\theta} = \boldsymbol{\theta}_0\).
- 局所識別: ヤコビアン \(G = \dfrac{\partial \boldsymbol{m}\left(\boldsymbol{\theta}_0\right)}{\partial \boldsymbol{\theta}^\top}\) (\(p \times k\) 行列) のランクが \(k\) (ランク条件) であれば, \(\boldsymbol{\theta}_0\) の近傍で識別される.
ランク条件が成り立つには少なくとも \(p \geq k\) (次数条件) が必要である.
ランク条件の解釈
真値 \(\boldsymbol{\theta}_0\) の近くでは, モーメント写像は1次近似できます. \[ \boldsymbol{m}\left(\boldsymbol{\theta}\right) \approx \boldsymbol{m}\left(\boldsymbol{\theta}_0\right) + G\left(\boldsymbol{\theta} - \boldsymbol{\theta}_0\right). \] 推定は, データから得た \(\hat{\boldsymbol{m}}\) に \(\boldsymbol{m}\left(\boldsymbol{\theta}\right)\) を合わせる \(\boldsymbol{\theta}\) を探すことなので, 真値の近くでは連立1次方程式 \[ G\left(\boldsymbol{\theta} - \boldsymbol{\theta}_0\right) = \hat{\boldsymbol{m}} - \boldsymbol{m}\left(\boldsymbol{\theta}_0\right) \] を \(\boldsymbol{\theta} - \boldsymbol{\theta}_0\) について解く問題に帰着します. この解 (\(p > k\) なら両辺の距離を最小にする解) が一意に決まるのは, \(G\) の \(k\) 本の列が1次独立なとき, すなわち \(\operatorname{rank} G = k\) のときです. これがランク条件の意味です.
パラメータが2つ \(\boldsymbol{\theta} = \left(\theta_1, \theta_2\right)^\top\), モーメントが2本の場合を考えます. 真値からのずれを \(\Delta\boldsymbol{\theta} = \boldsymbol{\theta} - \boldsymbol{\theta}_0\), \(\Delta\boldsymbol{m} = \hat{\boldsymbol{m}} - \boldsymbol{m}\left(\boldsymbol{\theta}_0\right)\) と書くと, 解くのは \(G \Delta\boldsymbol{\theta} = \Delta\boldsymbol{m}\) です.
ヤコビアンのランクが2の場合. \[ \begin{pmatrix} \Delta m_1 \\ \Delta m_2 \end{pmatrix} = \begin{pmatrix} 1 & 1 \\ 1 & -1 \end{pmatrix} \begin{pmatrix} \Delta\theta_1 \\ \Delta\theta_2 \end{pmatrix} \quad \Longrightarrow \quad \Delta\theta_1 = \frac{\Delta m_1 + \Delta m_2}{2}, \qquad \Delta\theta_2 = \frac{\Delta m_1 - \Delta m_2}{2}. \] 2本の式は独立な情報を持っているので, \(\Delta \boldsymbol{\theta}\) は一意に決まります.
また, ヤコビアンのランクが1の場合, \[ \begin{pmatrix} \Delta m_1 \\ \Delta m_2 \end{pmatrix} = \begin{pmatrix} 1 & 1 \\ 2 & 2 \end{pmatrix} \begin{pmatrix} \Delta\theta_1 \\ \Delta\theta_2 \end{pmatrix}. \] これは自明に一意に定まりません.
識別の強さと特異値
ここまでは 「解ける・解けない」という2値の話でした. しかし, 実際に数値計算の場面で遭遇するのは, その中間である 「理屈上解けはするが, 現実的な時間で収束しない」という 「弱識別」の状態になります. 極端な話, 上のランク条件は, 2本目の式の係数を \(\left(2, 2\right)\) から \(\left(2, 2.001\right)\) にズラせばランク条件を満たしたことになります. しかし, 2本目の式は1本目の式とほとんど同じなので, 識別という観点からするとほとんど意味がありません. そのため, ランクという離散的な判定ではなく, 「いちばん効かない方向でも, モーメントはどれだけ動くのか」という連続的な量が要ります. それを与えるのが特異値分解です.
命題 4.3 (特異値分解) \(p \times k\) 行列 \(G\) (\(p \geq k\)) は, 列が正規直交な \(U\) (\(p \times k\)) と直交行列 \(V\) (\(k \times k\)) を使って \[ G = U \Sigma V^\top, \qquad \Sigma = \operatorname{diag}\left(s_1, \ldots, s_k\right), \quad s_1 \geq \cdots \geq s_k \geq 0 \] と分解できる. \(s_j\) を特異値, \(V\) の第 \(j\) 列 \(\boldsymbol{v}_j\) を右特異ベクトルという. \(G \boldsymbol{v}_j = s_j \boldsymbol{u}_j\) であり, \(\boldsymbol{u}_j\) は単位ベクトルだから, \(s_j\) は「パラメータを \(\boldsymbol{v}_j\) 方向に1単位動かしたときにモーメントが動く大きさ」である. 一般に \[ \max_{\left\|\boldsymbol{v}\right\| = 1} \left\|G \boldsymbol{v}\right\| = s_1, \qquad \min_{\left\|\boldsymbol{v}\right\| = 1} \left\|G \boldsymbol{v}\right\| = s_k \] で, 最小は最終列 \(\boldsymbol{v}_k\) で達成される. したがって
- \(s_k > 0 \iff \operatorname{rank} G = k\): どの方向に動かしてもモーメントが動く (局所識別),
- \(s_k = 0 \iff G \boldsymbol{v}_k = \boldsymbol{0}\): \(\boldsymbol{v}_k\) 方向にはモーメントがまったく動かない (ランク落ち),
- \(s_k\) が \(s_1\) に比べて極端に小さい: \(\boldsymbol{v}_k\) 方向にはほとんど動かない (弱識別).
さきほどの2つの \(G\) に, 2行目を \(\left(2, 2.001\right)\) にずらしたものを加えた3つで計算します.
Code
Gs = ([1.0 1.0; 1.0 -1.0], [1.0 1.0; 2.0 2.001], [1.0 1.0; 2.0 2.0])
svd_rows = Any[]
for G in Gs
svd_G = svd(G)
push!(svd_rows, Any["($(G[2, 1]), $(G[2, 2]))", rank(G),
round(svd_G.S[1], digits=3), round(svd_G.S[2], sigdigits=3),
round(svd_G.S[1] / svd_G.S[2], sigdigits=3),
"("*join(round.(svd_G.V[:, end], digits=3), ", ")*")"])
end
pretty_table(permutedims(hcat(svd_rows...));
column_labels=[raw"$G$ の2行目", "ランク", raw"$s_1$", raw"$s_2$",
raw"$s_1 / s_2$", raw"$\boldsymbol{v}_2$"],
backend=:markdown, allow_markdown_in_cells=true)| \(G\) の2行目 | ランク | \(s_1\) | \(s_2\) | \(s_1 / s_2\) | \(\boldsymbol{v}_2\) |
|---|---|---|---|---|---|
| (1.0, -1.0) | 2 | 1.414 | 1.41 | 1.0 | (-0.0, -1.0) |
| (2.0, 2.001) | 2 | 3.163 | 0.000316 | 10000.0 | (-0.707, 0.707) |
| (2.0, 2.0) | 1 | 3.162 | 8.23e-17 | 3.84e16 | (-0.707, 0.707) |
| 注: どの \(G\) も1行目は \((1, 1)\). | |||||
1行目が局所識別であり, 特異値が十分大きい値であることがわかります. 3行目がランク落ちのケースであり, \(s_2\) がほとんどゼロであり, 識別されないことが表現されています (解析的には \(s_2 = 0\) ですが, 浮動小数点の計算ではごく小さい値になります). そして, 2行目が弱識別のケースです. ランク条件は満たしていますが, \(v_2\) の方向にほとんどモーメントが動かないことが \(s_2\) が非常に小さいことからわかります.
目的関数 \(\left\|G\left(\boldsymbol{\theta} - \boldsymbol{\theta}_0\right)\right\|^2\) の等高線を描くと, このことがそのまま図形になります. 等高線は \(\boldsymbol{v}_1, \ldots, \boldsymbol{v}_k\) を軸とする楕円で, 軸の長さは \(1/s_j\) に比例します. 特異値が小さい方向ほど楕円は細長く伸び, \(s_k = 0\) では無限に伸びて平行な溝になります.
Code
d_grid = range(-1, 1, length=201)
objective_lin(G, Δθ1, Δθ2) = sum(abs2, G * [Δθ1, Δθ2])
titles = ["rank 2", "rank 2, weak"]
ps = [contour(d_grid, d_grid, (Δθ1, Δθ2) -> objective_lin(G, Δθ1, Δθ2),
levels=12, colorbar=false, title=panel_title, aspect_ratio=1,
xlabel=L"\theta_1 - \theta_{10}",
ylabel=(k == 1 ? L"\theta_2 - \theta_{20}" : ""))
for (k, (G, panel_title)) in enumerate(zip(Gs[1:2], titles))]
scatter!(ps[1], [0], [0], ms=4, c=:black, label="")
plot!(ps[2], [-1, 1], [1, -1], lw=2, ls=:dash, c=:black, label="")
Plots.annotate!(ps[2], 0.42, -0.2, Plots.text(L"\mathbf{v}_2", 11))
plot(ps..., layout=(1, 2), size=(540, 300), bottom_margin=3Plots.mm,
left_margin=3Plots.mm)注: Contours of \(|G(\boldsymbol{\theta} - \boldsymbol{\theta}_0)|^2\) around the truth, for the first two \(G\) of 表 4.2. Both have full rank, but at this scale the right one looks flat along \(\boldsymbol{v}_2\).
左は同心円で, 底が1点に定まっています. 右は \(\boldsymbol{v}_2\) 方向に伸びた溝で, 目で見るかぎり, 谷底に沿ってどこまで動いても目的関数の値は変わりません. 実際には \(\boldsymbol{v}_2\) 方向にも \(\left(3.2 \times 10^{-4}\right)^2\) の傾きがあり, 谷は \(\left|\boldsymbol{\theta} - \boldsymbol{\theta}_0\right| \sim 3000\) のはるか遠方で閉じているのですが, この縮尺ではランク落ち (\(s_2 = 0\)) の平らな溝とほぼ同じものです.
最適化アルゴリズムは谷底のどこかで「収束」しますが, どこに止まるかを決めているのは初期値と丸め誤差です. 識別ができているかどうかの判定は, 収束した点の近くでのヤコビアンの特異値を見ればできます. また, 特異値の計算は候補となるパラメータ値でモデルを微分するだけなので, 推定の前にモデルとモーメントの設計を点検する道具になります.
設計を点検する
標準誤差を計算するときに, 7本のモーメントの5つのパラメータに関するヤコビアン \(\hat{G}\) を自動微分で求めてありました. 同じものを識別の点検に使い回します. 目的関数と同じ距離で測るため, 行を各モーメントの標準誤差で割り, 列をパラメータの値でスケールしたうえで, 特異値を見ます.
G_scaled = (Ĝ ./ ses_m) .* θ̂s'
svals_gmm = svdvals(G_scaled)
(singular_values=round.(svals_gmm, sigdigits=3),)(singular_values = [8760.0, 4070.0, 2230.0, 108.0, 15.6],)特異値は5つともゼロから明確に離れており, ランク条件は満たされています. 前節の推定がうまくいったのは幸運ではなかったわけです.
ヤコビアンの可視化
ヤコビアンには, ランク条件の成否以上の情報が詰まっています. 「パラメータを1%動かしたとき, 各モーメントが標準誤差の何倍動くか」を行列として可視化すると, どのモーメントがどのパラメータの情報源かが一目で分かります.
Code
G_percent = (Ĝ .* θ̂s') ./ ses_m ./ 100
moment_plot_labels_gmm = [L"e_W \times 1", L"e_W \times \ln w_H",
L"e_H \times 1", L"e_H \times \ln w_W",
L"e_3 \times 1", L"e_3 \times \ln w_W", L"e_3 \times \ln w_H"]
color_limit = maximum(abs, G_percent)
heatmap(1:5, 1:7, G_percent, c=:RdBu, clims=(-color_limit, color_limit),
colorbar=false, xticks=(1:5, param_plot_labels_gmm),
yticks=(1:7, moment_plot_labels_gmm), yflip=true, size=(520, 340))
# PrettyTables also exports annotate!, so name the module
Plots.annotate!(vec([(l, k, text(round(G_percent[k, l], digits=1), 7,
abs(G_percent[k, l]) > 0.6 * color_limit ? :white : :black))
for k in 1:7, l in 1:5]))
plot!(xlabel="Parameter", ylabel="Moment")注: Elasticity of each moment with respect to each parameter.
図 4.4 は見事なブロック構造をしています. 妻の MRS 条件 (\(e_W\) の行) は \(\left(\phi_W, \alpha_W\right)\) だけを, 夫の行は \(\left(\phi_H, \alpha_H\right)\) だけを動かし, \(\mu\) の列は \(e_3\) の行にしか値を持ちません. 各自の MRS 条件は一階条件の比をとって \(\lambda\) と一緒に \(\mu\) も消した式 (式 4.5) なので, 交渉力の情報は夫婦の限界効用を突き合わせる 式 4.6 だけが持っているのです. 各自の選好は各自の MRS 条件から, 交渉力は夫婦間の条件から, という識別の分業が, 行列の形としてそのまま見えています. \(e_W \times \ln w_H\) の行がまるごとゼロなのは, 2つの賃金が独立で, しかも夫の対数賃金の平均がゼロだからです. この条件は有効 (残差と直交する) ではありますが, パラメータをまったく動かさないので情報を持っていません. 操作変数が有効であることと, 有益であることは別だということです. 感応度を定量化して見せるこの種の分析は, Andrews et al. (2017) が構造推定の透明性のために提案したものと同じ発想であり, モーメントの取捨選択の指針になります.
過剰識別の検定
4.1.1節の \(J\) 統計量を計算します. 効率的な重み行列 \(\hat{S}^{-1}\) は前項で作ったものをそのまま使います. モーメントが7本, パラメータが5個なので, 自由度は \(7 - 5 = 2\) です.
J = N * gbar(θ̂s)' * inv(Ŝ) * gbar(θ̂s)
println("J = ", round(J, digits=2),
" χ²(2) の5%臨界値 = ", round(quantile(Chisq(2), 0.95), digits=2),
" p 値 = ", round(1 - cdf(Chisq(2), J), digits=3))J = 1.39 χ²(2) の5%臨界値 = 5.99 p 値 = 0.498\(J\) は臨界値を大きく下回り, 過剰識別の制約はデータと整合的です. データをモデルから作ったのだから当然ではありますが, 実装が意図どおり動いていることの確認にはなります. 前項で述べたとおり, 検定が通ったからといってモデルが正しいと結論できるわけではありません.
4.2 最小距離推定と SMM
前節の GMM が機能したのは, モーメント関数 \(\boldsymbol{g}\left(\boldsymbol{x}_i, \boldsymbol{\theta}\right)\) (残差 \(\times\) 操作変数) が観測データとパラメータだけで計算できたからです. 一階条件を「残差 = 0」の形に書き直し, 世帯ごとの数字を代入すれば評価できました. 構造推定では, これが書けることの方がむしろ例外です.
書けなくなる理由は2つあります. 第一に, モデルの含意が数値的にしか出せない場合です. 動学モデルの政策関数, 均衡モデルの価格, 定常状態から計算される結婚率や失業率は, どれもパラメータを与えるたびにモデルを解き直さないと分かりません. データ側の結婚率は調査から数えれば済むのに, それに対応するモデル側の量は閉じた形で書けない, という非対称が生じます. しかも結婚率は誰か1人について定義される量ではないので, 個人 \(i\) ごとの残差という形にそもそもなりません. 第二に, 世帯 \(i\) について観測できない変数がモデルに入る場合です. 未観測の異質性や潜在的な状態変数があると, 観測できるのはそれらについて平均を取った後の分布だけになります. どちらも, モーメント関数 \(\boldsymbol{g}\left(\boldsymbol{x}_i, \boldsymbol{\theta}\right)\) が個人ごとに書けないという点で同じであり, モデルが含意する統計量を元に推定する必要があります.
定義 4.3 (最小距離推定) データから計算した統計量 \(\hat{\boldsymbol{m}} = \frac{1}{N}\sum_i \boldsymbol{h}\left(\boldsymbol{x}_i\right)\) と, モデルが含意するその対応物 \(\boldsymbol{m}\left(\boldsymbol{\theta}\right)\) の距離を最小化する推定量を最小距離推定量 (minimum distance estimator) という. \[ \hat{\boldsymbol{\theta}} = \arg\min_{\boldsymbol{\theta}}\ \left(\hat{\boldsymbol{m}} - \boldsymbol{m}\left(\boldsymbol{\theta}\right)\right)^\top W \left(\hat{\boldsymbol{m}} - \boldsymbol{m}\left(\boldsymbol{\theta}\right)\right). \] これは \(\boldsymbol{g}\left(\boldsymbol{x}_i, \boldsymbol{\theta}\right) = \boldsymbol{h}\left(\boldsymbol{x}_i\right) - \boldsymbol{m}\left(\boldsymbol{\theta}\right)\) と置いた GMM である.
命題 4.4 (最小距離推定量の漸近分布と最適な重み) 命題 4.1 の条件が成り立つとし, \(G = \partial \boldsymbol{m}\left(\boldsymbol{\theta}_0\right) / \partial \boldsymbol{\theta}^\top\) (\(p \times k\)), \(V = \lim_{N \to \infty} \operatorname{Var}\left(\sqrt{N}\left(\hat{\boldsymbol{m}} - \boldsymbol{m}\left(\boldsymbol{\theta}_0\right)\right)\right)\) とおくと, \[ \sqrt{N}\left(\hat{\boldsymbol{\theta}} - \boldsymbol{\theta}_0\right) \xrightarrow{\ d\ } \mathcal{N}\left(\boldsymbol{0},\ \Lambda V \Lambda^\top\right), \qquad \Lambda = \left(G^\top W G\right)^{-1} G^\top W. \tag{4.8}\] 特に \(W = V^{-1}\) のとき, 分散は \(\left(G^\top V^{-1} G\right)^{-1}\) に簡約され, これが \(W\) の選び方の中で最小になる. すなわち, 最適な重み行列はデータのモーメントの分散共分散行列の逆行列である.
これは 命題 4.1 を最小距離推定に当てはめただけです. \(\boldsymbol{g}\left(\boldsymbol{x}_i, \boldsymbol{\theta}\right) = \boldsymbol{h}\left(\boldsymbol{x}_i\right) - \boldsymbol{m}\left(\boldsymbol{\theta}\right)\) なので \(\bar{\boldsymbol{g}}_N\left(\boldsymbol{\theta}\right) = \hat{\boldsymbol{m}} - \boldsymbol{m}\left(\boldsymbol{\theta}\right)\) です. したがって 命題 4.1 のヤコビアンは \(-G\) で, モーメントの分散共分散行列は \[ V = \lim_{N \to \infty} \operatorname{Var}\left(\sqrt{N}\left(\hat{\boldsymbol{m}} - \boldsymbol{m}\left(\boldsymbol{\theta}_0\right)\right)\right) = \lim_{N \to \infty} \operatorname{Var}\left(\sqrt{N}\, \hat{\boldsymbol{m}}\right) \tag{4.9}\] です. \(\boldsymbol{m}\left(\boldsymbol{\theta}_0\right)\) は確率変数ではなくただの定数なので, 分散を取ると消えます. つまり \(V\) はデータのモーメント \(\hat{\boldsymbol{m}}\) の分散共分散行列そのもので, モデルにもパラメータにも依存しません. 式 4.1 に代入すると \(-G\) の符号は打ち消し合って \(\left(G^\top W G\right)^{-1} G^\top W V W G \left(G^\top W G\right)^{-1} = \Lambda V \Lambda^\top\) となり, 最適な重みも 命題 4.1 の \(W = V^{-1}\) そのままです. 最適な重みがパラメータに依存せずデータから計算できるので, 最小距離推定の場合は2段階推定をする必要はありません.
\(\Lambda\) は「データのモーメントが1単位動いたとき推定値がどれだけ動くか」を並べた \(k \times p\) 行列で, 推定量は漸近的に \(\hat{\boldsymbol{\theta}} - \boldsymbol{\theta}_0 \simeq \Lambda\left(\hat{\boldsymbol{m}} - \boldsymbol{m}\left(\boldsymbol{\theta}_0\right)\right)\) と書けます. モーメントの誤差の分散 \(V\) が \(\Lambda\) を通って推定値の分散になる, ということです. この \(\Lambda\) は, 標準誤差 (4.2.2節) にも感応度 (4.2.4.2節) にも出てきます. 最適な重み \(V^{-1}\) は, 分散の小さい, つまり精度よく測れているモーメントのずれほど重く数え, 互いに相関するモーメントは重なっている情報を割り引いて数える重み付けです.
モデル側 \(\boldsymbol{m}\left(\boldsymbol{\theta}\right)\) を, シミュレーションによって計算する場合を, とくに SMM (simulated method of moments) と呼びます. 推定量が変わるわけではなく, \(\boldsymbol{m}\left(\boldsymbol{\theta}\right)\) の計算方法の名前です. 以下の例ではモデル側の量が2次元の積分になるので, シミュレーション (モンテカルロ法) は使わず求積で計算します.
4.2.1 モデル
静学的な労働供給を考えます. 個人は消費 \(c\) と労働時間 \(h\) を選びます. 労働時間は週40時間を1とする単位で測ります. \[ \max_{c,\, h} U(c, h; \psi_i) = \frac{c^{1-\gamma}}{1-\gamma} - \psi_i \frac{h^{1+1/\eta}}{1+1/\eta} \qquad \text{s.t.} \quad c = w h \tag{4.10}\] \(\psi_i\) は働くことの不効用の重み, \(\gamma\) は消費の曲率, \(\eta\) は労働時間の Frisch 弾力性です. 一階条件は \[ w \left(w h\right)^{-\gamma} = \psi_i\, h^{1/\eta} \tag{4.11}\] で, これは閉じた形で解けます. \[ h^*\left(w, \psi\right) = \left(\frac{w^{1-\gamma}}{\psi}\right)^{\frac{\eta}{1 + \gamma}}. \tag{4.12}\] この解を希望労働時間と呼ぶことにします. \(\gamma > 1\) なら \(w\) の指数が負で, 賃金が上がると希望労働時間は減ります. 所得効果が代替効果に勝つ領域です.
ここで, 個人の異質性を考えます. 賃金 \(w_i\) と労働の不効用の重み \(\psi_i\) を独立な確率変数とします. \[ \begin{aligned} \log w_i \sim \mathcal{N}\left(\mu_w, \sigma_w^2\right), \qquad \log \psi_i \sim \mathcal{N}\left(\mu_\psi, \sigma_\psi^2\right) \end{aligned} \tag{4.13}\] 式 4.12 の対数を取れば \(h_i^* := \log h^*(w_i, \psi_i)\) は \(\log w_i\) と \(\log \psi_i\) の線形結合なので, 希望労働時間の分布も対数正規です.
ここに標準労働時間を入れます. 仕事には標準の週労働時間 \(\bar{H} = 1\) (40時間) があり, 何もしなければ人は \(\bar{H}\) で働きます. 希望労働時間 \(h^*_i\) で働くには, 時間を自由に選べる仕事を探すか, 雇い主と交渉するかが必要で, そのコスト \(\kappa_i\) は人によって異なるとします. 希望時間で働くことの利得を \[ g(w_i, \psi_i) := U\left(h^*_i; w_i, \psi_i\right) - U\left(\bar{H}; w_i, \psi_i\right) \geq 0, \tag{4.14}\] と書くと, 人は \(g(w_i, \psi_i) > \kappa_i\) のときだけ標準労働時間から外れます. コスト \(\kappa_i\) を平均 \(\bar\kappa\) の指数分布とすれば, 外れる確率は 式 4.15 です. \[ \Pr\left(h_i = h^*_i \mid w_i, \psi_i\right) = 1 - \exp\left(-g(w_i, \psi_i) / \bar\kappa\right) \tag{4.15}\] 利得が大きい人ほど, つまり希望時間が40時間から遠い人ほど外れやすく, 希望時間が40時間に近い人は外れる価値がないのでそのまま40時間で働きます. \(\mu_w = 0\) に基準化し, 推定するのは \(\boldsymbol{\theta} = \left(\mu_\psi, \sigma_w, \sigma_\psi, \bar\kappa\right)\) の4つで, 曲率 \(\gamma\) と Frisch 弾力性 \(\eta\) は文献から値を借りてくることにします.
なお, 40時間への集中をモデルに入れる方法はいくつかあり, 代表的な研究として Bick et al. (2022) や Ma (2024) があり, Yanagimoto (2026) も凸型の所得関数を使って40時間の集中を説明しています.
γ_ls, η_ls = 1.5, 0.5 # curvature and Frisch elasticity
H̄ = 1.0 # the standard workweek (40 hours)
hours(w; ψ, γ=γ_ls, η=η_ls) = (w^(1 - γ) / ψ)^(1 / (1 / η + γ))
utility(h, w; ψ, γ=γ_ls, η=η_ls) =
(w * h)^(1 - γ) / (1 - γ) - ψ * h^(1 + 1 / η) / (1 + 1 / η)\(\eta = 0.5\) は Chetty et al. (2011) の値です. 彼らは準実験的なミクロ研究のメタ分析から, intensive margin の Frisch 弾力性を 0.54 と集計し, マクロモデルのカリブレーションには 0.5 を使うよう勧めています. \(\gamma = 1.5\) は前節と同じ値です.
Code
w_grid = range(0.3, 3.0, length=120)
plot(w_grid, [hours(w; ψ=0.5) for w in w_grid], lw=1.8, label=L"\psi = 0.5")
plot!(w_grid, [hours(w; ψ=1.0) for w in w_grid], lw=1.8, ls=:dash,
label=L"\psi = 1.0")
hline!([H̄], c=:black, ls=:dot, label="40 hours")
plot!(xlabel=L"wage $w$", title=L"Desired Hours $h^*(w, \psi)$ (40h = 1)")注: Hours from 式 4.12 as a function of the wage, for two values of the disutility of work. The dotted line is the standard workweek.
賃金が上がると希望労働時間が減ることがわかります. \(\psi\) を上げると, 働くことが苦になる分, 曲線全体が下に移動します. 曲線が点線と交わるあたりの人は希望時間が40時間に近く, 外れる利得 式 4.14 が小さいので, ほぼ確実に40時間で働きます.
この構造が最小距離推定を要求します. 前節では一階条件を個人ごとに評価できました. ここでは, 40時間で働いている人について分かるのは「外れる利得がコストを下回った」ことだけで, その人の希望時間も, したがって選好も, 点としては復元できません. 個人ごとの残差が書けないので, 残るのは分布どうしを突き合わせる最小距離推定です. たとえば平均労働時間のモデル側の対応物は 式 4.16 のような二重積分です. \[ \mathbb{E}\left[h\right] = \iint \left[\left(1 - e^{-g(w, \psi) / \bar\kappa}\right) h^*\left(w, \psi\right) + e^{-g(w, \psi) / \bar\kappa}\, \bar{H}\right] dF_w\left(w\right) dF_\psi\left(\psi\right) \tag{4.16}\] \(h^*\) は閉じた形でも, 混合の重み \(e^{-g(w, \psi)/\bar\kappa}\) が \(\left(w, \psi\right)\) の非線形な関数なので, 積分は数値的にしか計算できません.
データ
個票を使います. Blackburn and Neumark (1992) が用いた全米縦断調査 (National Longitudinal Survey) の若年男性コホートの1980年時点で, 就業していた男性 935人について週平均労働時間と月収が記録されています.2
労働時間はモデルに合わせて, 40時間を1とする単位で測ります. 標本は全員が就業しているので, 参加の意思決定を持たない 式 4.10 のモデルと対象が揃っています.
nls = CSV.read(
@projectroot("lesson", "numerical", "static", "data", "nls80.csv"),
DataFrame)
hs_ls = nls.hours ./ 40 # hours, in 40-hour weeks
ys_ls = log.(nls.earnings) # log monthly earnings
ds_ls = Float64.(nls.hours .== 40) # at the standard workweek
N_ls = nrow(nls)
data_moments(hs, ys, ds) = [mean(hs), std(hs), std(ys), cor(hs, ys), mean(ds)]
m̂s = data_moments(hs_ls, ys_ls, ds_ls)Code
histogram(nls.hours, bins=(20:2:82) .- 1, label="", lc=:white, lw=0.3,
xlabel="usual weekly hours", ylabel="men")注: 935 working men, NLS 1980. More than half report exactly 40 hours.
図 4.6 がこの節のモデルを決めた事実です. 標本の 53% がちょうど40時間で, 40時間未満は 6% しかいません. 残りの 41% は40時間より長く, 45, 50, 55, 60 時間といった切りのよい値に集まっています.
ターゲットとするモーメントは 表 4.3 にまとめました.
Code
moment_labels_ls = ["週間就業時間の平均 \$\\bar{h}\$",
"週間就業時間の標準偏差 \$s_h\$",
"対数月収の標準偏差 \$s_y\$",
"労働時間と対数月収の相関 \$r_{hy}\$",
"ちょうど40時間の割合 \$\\bar{d}\$"]
m̂s_raw = [m̂s[1] * 40, m̂s[2] * 40, m̂s[3], m̂s[4], m̂s[5]]
pretty_table(hcat(moment_labels_ls, round.(m̂s_raw, digits=4));
column_labels=["モーメント", "データの値"],
backend=:markdown, allow_markdown_in_cells=true)| モーメント | データの値 |
|---|---|
| 週間就業時間の平均 \(\bar{h}\) | 43.9294 |
| 週間就業時間の標準偏差 \(s_h\) | 7.2243 |
| 対数月収の標準偏差 \(s_y\) | 0.4211 |
| 労働時間と対数月収の相関 \(r_{hy}\) | -0.0472 |
| ちょうど40時間の割合 \(\bar{d}\) | 0.5316 |
| 注: 標本は NLS 若年男性コホートの1980年時点の就業者 935人. | |
目を引くのは4本目で, 労働時間と対数月収の相関は -0.047, つまりほぼゼロです. 長く働く人ほど稼ぐわけではない, ということです. 散らばりの源が賃金だけならこうはなりません. 賃金が高い人は (このモデルでは) 働く時間が短く, かつ所得は高いので, 相関は明確に負になるはずです.
モーメント条件
式 4.16 の積分は2次元です. 賃金も選好も対数正規なので, どちらも \(\exp\left(\mu + \sqrt{2}\,s\, z\right)\) と変換すればガウス・エルミート求積 (2.3節) がそのまま使えます. 2つの積分は独立なので, グリッドの直積を取るだけです. 各節点では希望時間 \(h^*\) と利得 \(g\) を計算し, 式 4.15 の確率で「希望時間で働く」と「40時間で働く」の2つの枝に重みを分けます. 被積分関数は滑らかなので, 各次元10点で足ります.3
μw_ls = 0.0 # the wage unit, normalised
function model_moments(θs; μw=μw_ls, γ=γ_ls, η=η_ls, n_nodes=10)
μψ, σw, σψ, κ̄ = θs
xs, ωs = gausshermite(n_nodes)
ωs = ωs ./ sqrt(π) # weights sum to 1
mean_h = mean_h2 = mean_y = mean_y2 = mean_hy = mean_d = 0.0
for i in 1:n_nodes, j in 1:n_nodes # wage node i, taste node j
w = exp(μw + sqrt(2) * σw * xs[i])
ψ = exp(μψ + sqrt(2) * σψ * xs[j]) # μψ is the mean of log ψ
h★ = hours(w; ψ=ψ, γ=γ, η=η)
gain = utility(h★, w; ψ=ψ, γ=γ, η=η) -
utility(H̄, w; ψ=ψ, γ=γ, η=η)
prob_deviate = 1 - exp(-gain / κ̄) # @eq-ls-deviate
for (prob, h) in ((prob_deviate, h★), (1 - prob_deviate, H̄))
weight = ωs[i] * ωs[j] * prob
y = log(w * h)
mean_h += weight * h
mean_h2 += weight * h^2
mean_y += weight * y
mean_y2 += weight * y^2
mean_hy += weight * h * y
end
mean_d += ωs[i] * ωs[j] * (1 - prob_deviate)
end
var_h, var_y = mean_h2 - mean_h^2, mean_y2 - mean_y^2
corr_hy = (mean_hy - mean_h * mean_y) / sqrt(var_h * var_y)
return [mean_h, sqrt(var_h), sqrt(var_y), corr_hy, mean_d]
end4.2.2 推定
重み行列 \(\hat{V}\)
命題 4.4 の通り, 最適な重みはデータモーメントの分散共分散行列の逆行列です. GMM では \(\bar{\boldsymbol{g}}_N\) が個人ごとの \(\boldsymbol{g}\left(\boldsymbol{x}_i, \hat{\boldsymbol{\theta}}\right)\) の標本平均だったので, その標本分散共分散 \(\hat{S}\) がそのまま \(V\) の推定量になりました. 最小距離推定でも, モーメントが標本平均なら同じです. 5本のうち1本目と5本目はそうなっていて, \(h_i\) と \(d_i\) の標本分散を \(N\) で割れば済みます. しかし残りの3本, 2つの標準偏差と相関は標本平均の非線形な関数で, 定義 4.3 のように \(\frac{1}{N}\sum_i \boldsymbol{h}\left(\boldsymbol{x}_i\right)\) の形には書けません. この場合, ブートストラップを用いるか影響関数を使ったデルタ法で分散共分散行列を計算します. 表 4.4 を使った結果です.
function bootstrap_moments(R, hs, ys, ds) # keep the loop inside a function
N = length(hs)
moments_resampled = Matrix{Float64}(undef, R, 5)
for r in 1:R
resampled = rand(1:N, N) # resample men, with replacement
moments_resampled[r, :] =
data_moments(hs[resampled], ys[resampled], ds[resampled])
end
return moments_resampled
end
Random.seed!(1234)
R_moments = 1000
V̂_bootstrap = cov(bootstrap_moments(R_moments, hs_ls, ys_ls, ds_ls))Code
moment_symbols_ls = [raw"$\bar{h}$", raw"$s_h$", raw"$s_y$", raw"$r_{hy}$",
raw"$\bar{d}$"]
pretty_table(hcat(moment_symbols_ls, round.(1e4 .* V̂_bootstrap, digits=2));
column_labels=["", moment_symbols_ls...],
backend=:markdown, allow_markdown_in_cells=true)| \(\bar{h}\) | \(s_h\) | \(s_y\) | \(r_{hy}\) | \(\bar{d}\) | |
|---|---|---|---|---|---|
| \(\bar{h}\) | 0.37 | 0.31 | 0.09 | -0.29 | -0.58 |
| \(s_h\) | 0.31 | 0.57 | 0.09 | -0.4 | -0.37 |
| \(s_y\) | 0.09 | 0.09 | 1.23 | -0.05 | -0.12 |
| \(r_{hy}\) | -0.29 | -0.4 | -0.05 | 13.39 | 0.27 |
| \(\bar{d}\) | -0.58 | -0.37 | -0.12 | 0.27 | 2.61 |
| 注: 935人を復元抽出して5本のモーメントを計算し直すことを1000回繰り返し, その分散共分散を取ったもの. 値は \(10^4\) 倍して表示している. 労働時間は40時間を1とする単位で測っている. | |||||
影響関数を使った方法は次のデルタ法の考え方を使います. 個人ごとの量を \(\boldsymbol{a}_i = \left(h_i, h_i^2, y_i, y_i^2, h_i y_i, d_i\right)^\top\) とかける場合を考えます. この時, 5本のモーメントは個人ごとの量の標本平均 \(\bar{\boldsymbol{a}} = \frac{1}{N}\sum_i \boldsymbol{a}_i\) の滑らかな関数 \(\hat{\boldsymbol{m}} = \boldsymbol{\phi}\left(\bar{\boldsymbol{a}}\right)\) と書けます. 例えば労働時間の標準偏差は \(\sqrt{\bar{a}_2 - \bar{a}_1^2}\) です. 標本平均 \(\bar{\boldsymbol{a}}\) の誤差は中心極限定理で扱えるので, あとは \(\phi(\cdot)\) のヤコビアンを考えればデルタ法 (命題 4.5) から \(\phi(\bar{\boldsymbol{a}})\) の誤差の分散共分散行列が計算できます.
命題 4.5 (デルタ法) 確率ベクトルの列 \(\hat{\boldsymbol{x}}_N\) が \(\sqrt{N}\left(\hat{\boldsymbol{x}}_N - \boldsymbol{\mu}\right) \xrightarrow{\ d\ } \mathcal{N}\left(\boldsymbol{0},\ \Sigma\right)\) を満たすとする. 関数 \(\boldsymbol{\phi}\) が \(\boldsymbol{\mu}\) で微分可能で, そのヤコビアンを \(D = \partial \boldsymbol{\phi}\left(\boldsymbol{\mu}\right) / \partial \boldsymbol{x}^\top\) とすると, \[ \sqrt{N}\left(\boldsymbol{\phi}\left(\hat{\boldsymbol{x}}_N\right) - \boldsymbol{\phi}\left(\boldsymbol{\mu}\right)\right) \xrightarrow{\ d\ } \mathcal{N}\left(\boldsymbol{0},\ D \Sigma D^\top\right). \]
これを個人ごとの量 \(\hat{\boldsymbol{\xi}}_i = \hat{D}\left(\boldsymbol{a}_i - \bar{\boldsymbol{a}}\right)\) の分散として計算します (アルゴリズム 4.2).
アルゴリズム 4.2 (影響関数による \(\hat{V}\) の計算)
- 個人ごとの量 \(\boldsymbol{a}_i\) を, モーメントがその標本平均の滑らかな関数 \(\hat{\boldsymbol{m}} = \boldsymbol{\phi}\left(\bar{\boldsymbol{a}}\right)\) になるように並べる.
- \(\boldsymbol{\phi}\) のヤコビアンを標本平均で評価した \(\hat{D} = \partial \boldsymbol{\phi}\left(\bar{\boldsymbol{a}}\right) / \partial \boldsymbol{a}^\top\) を求める.
- 個人ごとの影響関数 \(\hat{\boldsymbol{\xi}}_i = \hat{D}\left(\boldsymbol{a}_i - \bar{\boldsymbol{a}}\right)\) を計算する.
- その標本分散共分散を \(N\) で割る. \[ \hat{V} = \frac{1}{N}\, \widehat{\operatorname{Var}}\left(\hat{\boldsymbol{\xi}}_i\right) \tag{4.17}\]
手順2のヤコビアンは, GMM の \(\hat{G}\) と同じく, \(\boldsymbol{\phi}\) を Julia の関数として書けば ForwardDiff.jacobian で正確に求まるので, 微分を手で計算する必要はありません.
# φ: the five moments as a function of the sample means of a_i
moments_of_means(a) = [a[1], sqrt(a[2] - a[1]^2), sqrt(a[4] - a[3]^2),
(a[5] - a[1] * a[3]) / sqrt((a[2] - a[1]^2) * (a[4] - a[3]^2)), a[6]]
function influence(hs, ys, ds) # row i: ξ̂_i' = (a_i - ā)' D̂'
a_obs = hcat(hs, hs .^ 2, ys, ys .^ 2, hs .* ys, ds) # row i: a_i'
ā = vec(mean(a_obs, dims=1))
D̂ = ForwardDiff.jacobian(moments_of_means, ā) # 5×6
return (a_obs .- ā') * D̂'
end
V̂ = cov(influence(hs_ls, ys_ls, ds_ls)) ./ N_ls # variance of the moments
σs = sqrt.(diag(V̂)) # s.e. of each momentCode
moment_corr = V̂ ./ (σs * σs') # quoted in the text below
pretty_table(hcat(moment_symbols_ls, round.(1e4 .* V̂, digits=2));
column_labels=["", moment_symbols_ls...],
backend=:markdown, allow_markdown_in_cells=true)| \(\bar{h}\) | \(s_h\) | \(s_y\) | \(r_{hy}\) | \(\bar{d}\) | |
|---|---|---|---|---|---|
| \(\bar{h}\) | 0.35 | 0.28 | 0.07 | -0.25 | -0.56 |
| \(s_h\) | 0.28 | 0.54 | 0.09 | -0.32 | -0.36 |
| \(s_y\) | 0.07 | 0.09 | 1.19 | 0.0 | -0.14 |
| \(r_{hy}\) | -0.25 | -0.32 | 0.0 | 14.94 | 0.19 |
| \(\bar{d}\) | -0.56 | -0.36 | -0.14 | 0.19 | 2.67 |
| 注: アルゴリズム 4.2 で計算した \(\hat{V}\). 値は \(10^4\) 倍して表示している. 労働時間は40時間を1とする単位で測っている. | |||||
平均労働時間と労働時間の標準偏差の相関が 0.64 と高いのは, どちらも同じ 935人の労働時間から作られているからです. 40時間の割合は労働時間の散らばりと負に相関します. 40時間に人が集まれば散らばりは小さくなるからです. この非対角成分が使えることが, 個票を持っていることの利点です.
点推定と標準誤差
重み行列を \(W = \hat{V}^{-1}\) として, 定義 4.3 の目的関数を最小化します. モデルのモーメントは求積の和なので \(\boldsymbol{\theta}\) について滑らかで, GMM のとき (4.1.3節) と同じく, 自動微分で求めた勾配を使う L-BFGS で最小化できます. 目的関数は \(J\) 統計量と同じ尺度の量なので, その変化が \(10^{-6}\) を下回ったら止めます (abstol). 大域識別の間接チェックを兼ねて, 離れた3つの初期値から解きます.
W_md = inv(V̂) # efficient weight, fixed
md_loss(m̂s, θs, W) = (gs=model_moments(θs) .- m̂s; dot(gs, W, gs))
# quasi-Newton with box constraints; the gradient comes from ForwardDiff
minimize4(f, θs0) = solve(
OptimizationProblem(
OptimizationFunction((θs, _) -> f(θs), Optimization.AutoForwardDiff()),
θs0; lb=[-5.0, 0.02, 0.01, 1e-4], ub=[5.0, 5.0, 5.0, 10.0]),
NLopt.LD_LBFGS(), maxiters=20_000, abstol=1e-6)
initial_values = ([0.0, 0.45, 0.5, 0.05], [-1.2, 0.3, 0.8, 0.005],
[1.1, 0.9, 0.2, 0.5])
solutions = [minimize4(θs -> md_loss(m̂s, θs, W_md), θs0)
for θs0 in initial_values]
θ̂s_ls = solutions[1].u
param_labels_ls = [raw"$\mu_\psi$", raw"$\sigma_w$", raw"$\sigma_\psi$",
raw"$\bar{\kappa}$"]Code
losses_initial = [md_loss(m̂s, θs0, W_md) for θs0 in initial_values]
start_row_labels = ["初期値 " .* param_labels_ls; "初期値での目的関数";
"収束先 " .* param_labels_ls; "収束先での目的関数"]
start_columns = [vcat(θs0, round(loss0, digits=1), round.(solution.u, digits=6),
round(solution.objective, sigdigits=6))
for (θs0, loss0, solution) in
zip(initial_values, losses_initial, solutions)]
pretty_table(hcat(start_row_labels, start_columns...);
column_labels=["", "初期値 1", "初期値 2", "初期値 3"],
backend=:markdown, allow_markdown_in_cells=true)| 初期値 1 | 初期値 2 | 初期値 3 | |
|---|---|---|---|
| 初期値 \(\mu_\psi\) | 0.0 | -1.2 | 1.1 |
| 初期値 \(\sigma_w\) | 0.45 | 0.3 | 0.9 |
| 初期値 \(\sigma_\psi\) | 0.5 | 0.8 | 0.2 |
| 初期値 \(\bar{\kappa}\) | 0.05 | 0.005 | 0.5 |
| 初期値での目的関数 | 326.2 | 4628.3 | 4529.6 |
| 収束先 \(\mu_\psi\) | -0.367811 | -0.367811 | -0.367811 |
| 収束先 \(\sigma_w\) | 0.456239 | 0.456239 | 0.456239 |
| 収束先 \(\sigma_\psi\) | 0.519161 | 0.519161 | 0.519161 |
| 収束先 \(\bar{\kappa}\) | 0.046842 | 0.046842 | 0.046842 |
| 収束先での目的関数 | 0.00911583 | 0.00911583 | 0.00911583 |
| 注: 目的関数は重み行列を \(W = \hat{V}^{-1}\) とした 定義 4.3 の距離. | |||
表 4.6 のとおり, 目的関数の値が 326 から 4628 の離れた点から出発しても, 3つとも同じ点に収束しました. 収束先は表示した小数第6位まで一致し, そこでの目的関数の値も同じです.
推定量の漸近分散は 式 4.8 で与えられます. 材料のうち \(\hat{V}\) はすでに手元にあるので, 新しく要るのはヤコビアンを推定値で評価した \(\hat{G} = \partial \boldsymbol{m}\left(\hat{\boldsymbol{\theta}}\right) / \partial \boldsymbol{\theta}^\top\) だけです. モデルのモーメントは Gauss–Hermite 求積の和として閉じた式で計算しているので, GMM のときと同じく ForwardDiff.jacobian で求まります. \(\hat{V}\) は \(1/N\) を含んだ量なので, 推定値の分散共分散行列は \(\hat{\Lambda} \hat{V} \hat{\Lambda}^\top\) で, \(\hat{\Lambda} = \left(\hat{G}^\top W \hat{G}\right)^{-1} \hat{G}^\top W\) です. \(W = \hat{V}^{-1}\) なので, これは \(\left(\hat{G}^\top \hat{V}^{-1} \hat{G}\right)^{-1}\) に簡約されます.
Ĝ = ForwardDiff.jacobian(model_moments, θ̂s_ls)
Λ = inv(Ĝ' * W_md * Ĝ) * (Ĝ' * W_md) # how θ̂s_ls moves with each moment
ses_analytic = sqrt.(diag(Λ * V̂ * Λ'))Code
pretty_table(hcat(param_labels_ls, round.(θ̂s_ls, digits=4),
round.(ses_analytic, sigdigits=2));
column_labels=["パラメータ", "推定値", "標準誤差"],
backend=:markdown, allow_markdown_in_cells=true)| パラメータ | 推定値 | 標準誤差 |
|---|---|---|
| \(\mu_\psi\) | -0.3678 | 0.02 |
| \(\sigma_w\) | 0.4562 | 0.012 |
| \(\sigma_\psi\) | 0.5192 | 0.018 |
| \(\bar{\kappa}\) | 0.0468 | 0.004 |
| 注: 重み行列は \(W = \hat{V}^{-1}\). 標準誤差は 式 4.8 による. | ||
対数不効用の平均 \(\hat{\mu}_\psi =\) -0.37 (重みの中央値は \(e^{\hat{\mu}_\psi} =\) 0.69), 対数賃金の標準偏差 \(\hat\sigma_w =\) 0.46, 対数不効用の標準偏差 \(\hat\sigma_\psi =\) 0.52, 逸脱コストの平均 \(\hat{\bar\kappa} =\) 0.047 です. 選好の散らばりは無視できる大きさではなく, 上下1標準偏差の人の間で不効用の重みが 2.8倍違います. 式 4.12 によれば, これは希望労働時間の 1.35倍の違いにあたります. \(\bar\kappa\) は効用の単位なので, それ自体を読むより, 含意される「標準から外れる人の割合」47% の方が分かりやすいでしょう.
モデルのフィット
モデルへの当てはまりは, モーメントの標準誤差で測れます.
Code
moment_plot_labels_ls = ["mean hours", "sd hours", "sd log earnings",
"corr(h, log earnings)", "share at 40 hours"]
deviations = (model_moments(θ̂s_ls) .- m̂s) ./ σs # fit, in moment s.e. units
scatter(deviations, 1:5, ms=6, label="")
vline!([0], c=:black, lw=1, label="")
plot!(xlabel="(model − data) / s.e. of the moment",
yticks=(1:5, moment_plot_labels_ls),
ylims=(0.5, 5.5), xlims=(-1, 1), size=(560, 280))注: Targeted moments, model minus data, in standard errors of the moment.
図 4.7 で見るように, 5本とも 0.1 標準誤差の内側に収まっています. パラメータ4つで5本を狙っているので完全には合いませんが, 符号が入り混じっているのはモデルがデータをうまくとらえていることの示唆です.
重みが \(\hat{V}^{-1}\) なので, 目的関数の最小値がそのまま Hansen の \(J\) 統計量になります (命題 4.2). \(\hat{V}\) がすでに \(1/N\) を含んでいるので, \(N\) を掛ける必要はありません.
J_ls = solutions[1].objective
pval_ls = ccdf(Chisq(5 - 4), J_ls)過剰識別の検定をします. \(J =\) 0.01 で, 自由度 \(5 - 4 = 1\) の \(\chi^2\) の 5% 点 3.84 を下回ります (\(p =\) 0.92). 5本のモーメントは, この4つのパラメータで整合的に説明できる範囲にあります. \(J\) がほとんどゼロなのは偶然ではありません. 4.2.4.2節で見るとおり, 持ち込んだ \(\eta = 0.5\) がこのデータの好む値とほぼ一致していることも関係しています.
4.2.3 標準誤差の別の求め方
4.2.2節の標準誤差は, データのモーメントの分散共分散行列 \(\hat{V}\) とモデル側のヤコビアン \(\hat{G}\) を 式 4.8 に入れたものでした. ここでは, 推定そのものをやり直すブートストラップと, 個票がなく \(\hat{V}\) が手に入らない場合の最悪ケースの標準誤差を見ます.
ブートストラップ
表 4.4 のブートストラップは, データのモーメントだけを作り直すものでした. 影響関数が使いにくいモーメントでも, そうして \(\hat{V}\) さえ作れれば, 標準誤差は 式 4.8 で解析的に求まります.4 標準誤差のためだけに推定をやり直す必要はありません. 推定ごとやり直すブートストラップが必要になる, あるいは有利になるのは, 次のような場面です.
- ヤコビアン \(\hat{G}\). モデル側の都合によってモーメントが滑らかでない場合, \(\hat{G}\) の推定が不安定になり, 構成した解析的な標準誤差が信頼できなくなります.
- 多段階推定. 最小距離推定自体の計算が重いので, モデルを解く前にデータから一部のパラメータを推定する場合がよくあります. この段階の標準誤差が 式 4.8 に含まれないため, 両段階をやり直すブートストラップが必要になります.
- 推定値. 推定が終わったあと, 反実仮想などを行なった結果に標準誤差をつけたい場合があります (厚生変化など). その場合, デルタ法で標準誤差を与えるにはモデル自体を関数と見なしてそのヤコビアンを求める必要があります. しかし, この数値微分は不安定になりがちです.
今回の授業の例ではどれにも当たらないので, ブートスラップを解析的な標準誤差の検算として行います. 935人を復元抽出で引き直し, 5本のモーメントを作り直し, 最適化までやり直す. これを繰り返して得られる \(\hat{\boldsymbol{\theta}}^{\left(r\right)}\) の散らばりが標準誤差です. 重み行列も, 引き直した標本で影響関数から作り直します. 推定の手順を丸ごと再現するのがブートストラップの考え方だからです. 重み行列の推定誤差は1次の漸近分布には効かないので, 全標本の \(\hat{V}^{-1}\) に固定しても漸近的には同じ結果になりますが, 小標本での性質までとらえたいなら作り直すのが筋です.
function bootstrap_md(R, hs, ys, ds, θs0) # keep the loop inside a function
N = length(hs)
estimates = Matrix{Float64}(undef, R, length(θs0))
for r in 1:R
resampled = rand(1:N, N) # resample men, with replacement
hs_r, ys_r, ds_r = hs[resampled], ys[resampled], ds[resampled]
m̂s_r = data_moments(hs_r, ys_r, ds_r)
W_r = inv(cov(influence(hs_r, ys_r, ds_r)) ./ N) # weight rebuilt
estimates[r, :] = minimize4(θs -> md_loss(m̂s_r, θs, W_r), θs0).u
end
return estimates
end
Random.seed!(1234)
R_boot = 250
boot_estimates = bootstrap_md(R_boot, hs_ls, ys_ls, ds_ls, θ̂s_ls)
ses_boot = vec(std(boot_estimates, dims=1))Code
pretty_table(hcat(param_labels_ls, round.(θ̂s_ls, digits=4),
round.(ses_analytic, sigdigits=3), round.(ses_boot, sigdigits=3),
round.(ses_boot ./ ses_analytic, digits=2));
column_labels=["パラメータ", "推定値", "解析的", "ブートストラップ", "比"],
backend=:markdown, allow_markdown_in_cells=true)| パラメータ | 推定値 | 解析的 | ブートストラップ | 比 |
|---|---|---|---|---|
| \(\mu_\psi\) | -0.3678 | 0.0198 | 0.0201 | 1.02 |
| \(\sigma_w\) | 0.4562 | 0.0121 | 0.0126 | 1.04 |
| \(\sigma_\psi\) | 0.5192 | 0.0181 | 0.0178 | 0.98 |
| \(\bar{\kappa}\) | 0.0468 | 0.00401 | 0.00395 | 0.99 |
| 注: ブートストラップ標準誤差は, 935人を復元抽出し, モーメントと重み行列の計算から最適化までやり直すことを250回繰り返したときの推定値の標準偏差. 比はブートストラップ標準誤差を解析的な標準誤差で割ったもの. | ||||
4つとも 4% 以内で一致していルコとが分かります. 一致するはずのものが一致した, というだけの結果ですが, 確かめる価値はあります. モーメントの関数 \(\boldsymbol{\phi}\) の書き方やデルタ法の適用のどこかを間違えていれば, ここでずれます.
Code
# Gaussian kernel density with Silverman's rule-of-thumb bandwidth
function kernel_density(xs; n_points=200)
interquartile_range = quantile(xs, 0.75) - quantile(xs, 0.25)
bandwidth = 0.9 * min(std(xs), interquartile_range / 1.34) *
length(xs)^(-1 / 5)
x_grid = range(minimum(xs) - 3 * bandwidth, maximum(xs) + 3 * bandwidth,
length=n_points)
densities = [mean(pdf.(Normal(x, bandwidth), xs)) for x in x_grid]
return x_grid, densities
end
ps_boot = []
for k in 1:4
x_grid, densities = kernel_density(boot_estimates[:, k])
p = plot(x_grid, densities, lw=1.8, label="",
title=[L"\mu_\psi", L"\sigma_w", L"\sigma_\psi", L"\bar{\kappa}"][k])
vline!(p, [θ̂s_ls[k]], c=:black, lw=2, label="")
vline!(p, θ̂s_ls[k] .+ [-1.96, 1.96] .* ses_analytic[k], c=:black, ls=:dash,
lw=1.5, label="")
push!(ps_boot, p)
end
plot(ps_boot..., layout=(2, 2), size=(700, 440), yticks=false,
bottom_margin=3Plots.mm)注: Kernel density of the 250 bootstrap estimates (Gaussian kernel, Silverman’s rule-of-thumb bandwidth). Solid vertical line: point estimate. Dashed lines: analytic 95% interval.
図 4.8 はブートストラップ推定値の分布をカーネル密度で描いたものです. 破線が解析的な95%区間で, 曲線がブートストラップ推定値のカーネル密度です. 分布はほぼ対称で, 区間の幅も揃っています.
最悪ケース
解析的な方法もブートストラップ法も, 個票が手元にあったから計算できました. しかし, 構造推定ではそうでない場合もよくあります. ターゲットが先行研究や公表統計の点推定で, 手元にあるのは各モーメントの標準誤差 (対角成分) だけ, というのがよくある状況です. もし, 複数の資料からモーメントを集めてくれば, それらの相関はまず手に入りません (Cocci and Plagborg-Møller 2024). この例でいえば, 労働時間は労働力調査から, 所得の分散は賃金構造基本統計調査から, というような組み合わせです.
対角成分しか分からないなら, 重み行列も対角にするしかありません. まずその推定量を作ります.
W_diag = Diagonal(1 ./ σs .^ 2) # all we can build from the s.e.
θ̂s_diag = minimize4(θs -> md_loss(m̂s, θs, W_diag), θ̂s_ls).u
Ĝ_diag = ForwardDiff.jacobian(model_moments, θ̂s_diag)
Λ_diag = inv(Ĝ_diag' * W_diag * Ĝ_diag) * (Ĝ_diag' * W_diag)
ses_indep = sqrt.(diag(Λ_diag * Diagonal(σs .^ 2) * Λ_diag'))Cocci and Plagborg-Møller (2024) は, この状況でも有効な「最悪ケース」の標準誤差を提案しました. 鍵になるのは, 推定量が漸近的にはモーメントの線形結合になるという事実です.
命題 4.6 (最悪ケース標準誤差 (Cocci and Plagborg-Møller 2024)) \(p\) 本のモーメント \(\hat{\boldsymbol{m}}\) をターゲットにした最小距離推定量を考える. \[ \hat{\boldsymbol{\theta}} = \arg\min_{\boldsymbol{\theta}} \left(\boldsymbol{m}\left(\boldsymbol{\theta}\right) - \hat{\boldsymbol{m}}\right)^\top W \left(\boldsymbol{m}\left(\boldsymbol{\theta}\right) - \hat{\boldsymbol{m}}\right) \] ここで, \(W\) は \(p \times p\) の半正定値な重み行列である.
\(G = \partial \boldsymbol{m}\left(\boldsymbol{\theta}_0\right) / \partial \boldsymbol{\theta}^\top\), \(V\) を \(\sqrt{N}\left(\hat{\boldsymbol{m}} - \boldsymbol{m}\left(\boldsymbol{\theta}_0\right)\right)\) の漸近分散共分散行列とする. この時, 関心のある量 \(r\left(\boldsymbol{\theta}\right)\) について, \(\boldsymbol{\lambda} = \partial r\left(\boldsymbol{\theta}_0\right) / \partial \boldsymbol{\theta}\), \(\boldsymbol{x} = W G \left(G^\top W G\right)^{-1} \boldsymbol{\lambda}\) とおくと, デルタ法 (命題 4.5) よりつぎの漸近分布が得られる. \[ \sqrt{N}\left(r\left(\hat{\boldsymbol{\theta}}\right) - r\left(\boldsymbol{\theta}_0\right)\right) \xrightarrow{d} \mathcal{N}\left(0,\ \boldsymbol{x}^\top V \boldsymbol{x}\right) \] \(V\) の対角成分 \(\left(\sigma_1^2, \ldots, \sigma_p^2\right)\) だけが既知のとき, \[ \max_{V \succeq 0,\ V_{jj} = \sigma_j^2}\ \boldsymbol{x}^\top V \boldsymbol{x} = \left(\sum_{j=1}^{p} \sigma_j \left|x_j\right|\right)^2 \] であり, 最大はすべてのモーメントが完全相関する時に達成される.
上限の直感は素朴です. \(V\) が半正定値なら相関係数は1を超えられない, すなわち \(\left|V_{ij}\right| \leq \sigma_i \sigma_j\) なので, \(\boldsymbol{x}^\top V \boldsymbol{x} = \sum_{i,j} x_i x_j V_{ij}\) は各項を \(\left|x_i\right| \left|x_j\right| \sigma_i \sigma_j\) で置き換えたものより大きくなれません. すべてのモーメントが完全相関するとき, その上限がちょうど達成されます.
実装では \(\sigma_j\) をモーメント \(j\) の標準誤差に置き換えれば, \(\sum_j \sigma_j \left|x_j\right|\) がそのまま \(r\left(\hat{\boldsymbol{\theta}}\right)\) の標準誤差になります. 必要なのは各モーメントの標準誤差とヤコビアンだけなので, 先行研究の公表値をターゲットにしたカリブレーションにもそのまま使えます.5 この最悪ケース標準誤差に基づく信頼区間は, どんな相関構造のもとでも漸近的に名目以上の被覆率を持ちます (完全相関のケースで等号). その代償として, 真の相関構造によっては保守的 (広すぎる区間) になります.
Code
ses_full = sqrt.(diag(Λ_diag * V̂ * Λ_diag')) # full V̂, same estimator
ses_worst = [sum(σs .* abs.(Λ_diag[k, :])) for k in 1:4]
pretty_table(hcat(param_labels_ls, round.(θ̂s_diag, digits=4),
round.(ses_full, sigdigits=2), round.(ses_indep, sigdigits=2),
round.(ses_worst, sigdigits=2),
round.(ses_worst ./ ses_indep, digits=2));
column_labels=["パラメータ", "推定値", "分散共分散", "無相関仮定",
"最悪ケース", "最悪 / 無相関"],
backend=:markdown, allow_markdown_in_cells=true)| パラメータ | 推定値 | 分散共分散 | 無相関仮定 | 最悪ケース | 最悪 / 無相関 |
|---|---|---|---|---|---|
| \(\mu_\psi\) | -0.3688 | 0.022 | 0.028 | 0.043 | 1.54 |
| \(\sigma_w\) | 0.4563 | 0.012 | 0.012 | 0.018 | 1.48 |
| \(\sigma_\psi\) | 0.5187 | 0.019 | 0.022 | 0.041 | 1.87 |
| \(\bar{\kappa}\) | 0.0468 | 0.004 | 0.005 | 0.0087 | 1.75 |
| 注: 推定量は, \(\hat{V}\) の対角成分の逆数を重み行列とした最小距離推定量. 「分散共分散」は \(\hat{V}\) を非対角成分まで使った標準誤差, 「無相関仮定」は非対角成分をゼロとみなした標準誤差, 「最悪ケース」は 命題 4.6 の上限. | |||||
表 4.9 は, 同じ推定量 \(\hat{\boldsymbol{\theta}}\) の標準誤差を, \(V\) の分散共分散行列について知っていることを変えて計算したものです. 左の列は \(\hat{V}\) を非対角成分まで使った通常の値で, 個票があるからこそ計算できます. 中の列は非対角成分をゼロとみなした値, 右の列は非対角成分が最悪だったとき (命題 4.6 の完全相関) の値です. 最悪ケースは無相関仮定の 1.5倍から 1.9倍になっています. これが相関を知らないことの代償です. 信頼区間をこれだけ広げておけば, 真の相関がどうであっても被覆率は名目以上になります. そしてこの代償には天井があります. Cocci and Plagborg-Møller (2024) が示すとおり, Hölder の不等式から最悪ケースは無相関仮定の高々 \(\sqrt{p}\) 倍で, モーメント5本なら \(\sqrt{5} \approx\) 2.24 倍です.
ここで注意したいのは無相関仮定の列です. この例では左の解析的な列と比べて 0.98倍から 1.24倍で, 大きくは外れていません. しかし相関を無視した標準誤差が正しい値を上回るのか下回るのかは, 一般には決まりません. Cocci and Plagborg-Møller (2024) のメニューコストモデルを使ったモンテカルロ実験では逆に, 無相関を仮定した信頼区間の被覆率が名目95%を下回り (パラメータによっては89%程度), 検定が過剰に棄却します (Cocci and Plagborg-Møller 2024, sec. 5.1). 相関の符号と感応度 \(\Lambda\) の符号の組み合わせ次第で, どちらにも転びうるということです. 相関が本当に手に入らないなら, 「独立とみなす」のではなく最悪ケースを併記するのが安全です.
4.2.4 識別
設計を点検する
前節と同じ道具で設計を点検します. 標準誤差を計算するときに求めたヤコビアン \(\hat{G}\) を, 目的関数と同じ距離で測り直します. 重みが \(\hat{V}^{-1}\) なので, 行はモーメントの標準誤差, 列はパラメータの値でスケールします.
G_scaled = (Ĝ ./ σs) .* θ̂s_ls'
svals_md = svdvals(G_scaled)
(singular_values=round.(svals_md, sigdigits=3),)(singular_values = [38.7, 34.4, 14.8, 8.64],)4つとも0から離れていて, 最大と最小の比は 4.5倍です. 全てのパラメータがよく識別されていることが示唆されています.
Code
G_percent = G_scaled ./ 100
param_plot_labels_ls = [L"\mu_\psi", L"\sigma_w", L"\sigma_\psi",
L"\bar{\kappa}"]
color_limit = maximum(abs, G_percent)
heatmap(1:4, 1:5, G_percent, c=:RdBu, clims=(-color_limit, color_limit),
colorbar=false, xticks=(1:4, param_plot_labels_ls),
yticks=(1:5, moment_plot_labels_ls), yflip=true, size=(600, 360))
# PrettyTables also exports annotate!, so name the module
Plots.annotate!(vec([(l, k, text(round(G_percent[k, l], digits=2), 8,
abs(G_percent[k, l]) > 0.6 * color_limit ? :white : :black))
for k in 1:5, l in 1:4]))
plot!(xlabel="Parameter", ylabel="Moment")注: Change of each moment (in s.e. units) per 1% change of each parameter.
図 4.9 も素直に解釈することができます. 平均労働時間の行は \(\mu_\psi\) に強く反応し, 対数所得の散らばりは \(\sigma_w\) によく反応しており, 労働時間の散らばりはほぼ \(\sigma_\psi\) だけに反応します.
決定的なのは相関の行です. 相関だけが, \(\sigma_w\) と \(\sigma_\psi\) に符号の違う反応を示します. 賃金の散らばりを増やすと相関は下がり, 選好の散らばりを増やすと相関は上がる. 他の行はどちらも同じ向きに動かすので, 散らばりの源が賃金なのか選好なのかを区別できません. 相関を1本足したことが, \(\sigma_w\) と \(\sigma_\psi\) を分けているということです. 最後の行では, 40時間の割合が \(\bar\kappa\) と \(\mu_\psi\) に反応しています. 外れるコストが高ければ40時間に留まる人が増え, 不効用の重みが動けば希望時間の分布全体がずれて, 40時間の近くにいる人の数が変わるからです.
感応度と頑健性
推定結果を報告するとき, 推定値と標準誤差のほかに何を添えるべきか. 構造推定の実務で定着している答えが2つあります. モーメントに対する感応度と, 外生的に定めた値に対する頑健性です.
第一の感応度は, 標準誤差の計算で使った \(\Lambda = \left(\hat{G}^\top W \hat{G}\right)^{-1} \hat{G}^\top W\) そのものです. Andrews et al. (2017) はこれを感応度 (sensitivity) と名づけ, 報告することを勧めています.6 意味は 命題 4.4 の線形化そのもので, モーメント \(j\) が \(\Delta\) だけずれたとき推定値 \(\hat\theta_k\) は \(\Lambda_{kj}\Delta\) だけ動く, という一次近似です. モーメントの1本が誤って測られている, あるいはモデルの仮定の1つが破れてそのモーメントだけがずれる, という状況を読者が想定したとき, それが推定値にどう伝わるかを \(\Lambda\) から読み取れます. 単位を揃えるために, 列はモーメントの標準誤差, 行は推定値の標準誤差でスケールしておきます (Andrews et al. 2017, sec. IV.A).
Code
Λ_se = (Λ .* σs') ./ ses_analytic # rows: parameters, cols: moments
sensitivity = permutedims(Λ_se) # same layout as the Jacobian plot
color_limit_λ = maximum(abs, sensitivity)
heatmap(1:4, 1:5, sensitivity, c=:RdBu, clims=(-color_limit_λ, color_limit_λ),
colorbar=false, xticks=(1:4, param_plot_labels_ls),
yticks=(1:5, moment_plot_labels_ls), yflip=true, size=(600, 360))
Plots.annotate!(vec([(l, k, text(round(sensitivity[k, l], digits=2), 8,
abs(sensitivity[k, l]) > 0.6 * color_limit_λ ? :white : :black))
for k in 1:5, l in 1:4]))
plot!(xlabel="Parameter", ylabel="Moment")注: Sensitivity \(\hat\Lambda\): change of each estimate (in its s.e. units) when each moment moves by one s.e. of the moment.
図 4.10 の読み方は 図 4.9 の逆向きです. 縦軸と横軸をそろえてあるので, 同じマスどうしを見比べられます. ヤコビアンは「パラメータを動かすとモーメントがどう動くか」, \(\Lambda\) は「モーメントが動くと推定値がどう動くか」を並べています. たとえば \(\hat\sigma_\psi\) は週間就業時間の標準偏差に, \(\hat{\bar\kappa}\) はちょうど40時間の割合にもっとも強く依存しています. \(\Lambda\) の計算に追加の費用はかからないので, 推定値の表と一緒に載せておく価値があります.
第二の頑健性は, 外生的に定めた値についてです. ここでは消費の曲率 \(\gamma = 1.5\) と Frisch 弾力性 \(\eta = 0.5\) の2つで, どちらも先行研究の値です. 外から定めた値を動かして結果がどう変わるかを示すのは, 構造推定の論文では付録の定番で, Bick et al. (2022) は所得関数の傾きを区間で動かして主結果が変わらないことを確かめ (Bick et al. 2022, sec. III.C), Ma (2024) は曲率 \(\gamma\) の代替値で弾力性を再計算しています (Ma 2024, sec. 5.2.2). 同じことをします.
Code
scenarios = ((raw"基準", (;)), (raw"$\gamma = 1.2$", (γ=1.2,)),
(raw"$\gamma = 1.8$", (γ=1.8,)), (raw"$\eta = 0.25$", (η=0.25,)),
(raw"$\eta = 1$", (η=1.0,)))
scenario_solutions = [
minimize4(θs -> (gs=model_moments(θs; overrides...) .- m̂s;
dot(gs, W_md, gs)), θ̂s_ls)
for (_, overrides) in scenarios]
scenario_rows = [Any[label, round.(sol.u[1:3], digits=3)...,
round(sol.u[4], digits=4), round(sol.objective, digits=2)]
for ((label, _), sol) in zip(scenarios, scenario_solutions)]
pretty_table(permutedims(hcat(scenario_rows...));
column_labels=["持ち込んだ値", param_labels_ls..., raw"$J$"],
backend=:markdown, allow_markdown_in_cells=true)| 持ち込んだ値 | \(\mu_\psi\) | \(\sigma_w\) | \(\sigma_\psi\) | \(\bar{\kappa}\) | \(J\) |
|---|---|---|---|---|---|
| 基準 | -0.368 | 0.456 | 0.519 | 0.0468 | 0.01 |
| \(\gamma = 1.2\) | -0.32 | 0.439 | 0.453 | 0.0381 | 31.47 |
| \(\gamma = 1.8\) | -0.404 | 0.459 | 0.572 | 0.0574 | 19.18 |
| \(\eta = 0.25\) | -0.587 | 0.448 | 0.795 | 0.0647 | 10.36 |
| \(\eta = 1\) | -0.262 | 0.459 | 0.375 | 0.0385 | 11.46 |
| 注: 各行は, 持ち込んだ値のうち1つだけを変えて推定し直した結果. 基準は \(\gamma = 1.5\), \(\eta = 0.5\). | |||||
どちらも当てはまりまで動かします. \(\gamma\) をどちらに動かしても \(J\) は悪化し, \(\eta\) は 0.5 から離れるほど悪化します. \(\eta = 0.25\) で 10.4, \(\eta = 1\) で 11.5 と, どちらも 5% 点 3.84 を超えて棄却されます. 持ち込んだ Frisch 弾力性について, データははっきり意見を持っているということです.
それなら \(\eta\) も推定してしまえばよい, と思うのは自然です. パラメータが5つ, モーメントが5本で丁度識別になるので, 過剰識別の検定は使えなくなりますが, 推定自体はできます.
minimize5(f, θs0) = solve(
OptimizationProblem(
OptimizationFunction((θs, _) -> f(θs), Optimization.AutoForwardDiff()),
θs0; lb=[-5.0, 0.02, 0.01, 1e-4, 0.05], ub=[5.0, 5.0, 5.0, 10.0, 5.0]),
NLopt.LD_LBFGS(), maxiters=20_000, abstol=1e-6)
model_moments_η(θs) = model_moments(θs[1:4]; η=θs[5])
θ̂s_η = minimize5(θs -> (gs=model_moments_η(θs) .- m̂s; dot(gs, W_md, gs)),
[θ̂s_ls; η_ls]).u
Ĝ_η = ForwardDiff.jacobian(model_moments_η, θ̂s_η)
ses_η = sqrt.(diag(inv(Ĝ_η) * V̂ * inv(Ĝ_η)')) # just identified: Λ = G^{-1}推定された Frisch 弾力性は \(\hat\eta =\) 0.51 (標準誤差 0.1) です. Chetty et al. (2011) のメタ分析の 0.54 と, 持ち込んだ 0.5 は, どちらも95%信頼区間 [0.31, 0.71] の中にあります. 準実験からの推定値と, 横断面の分布の形からの推定値が一致するのは, それ自体が1つの結果です. ただし標準誤差は他のパラメータより大きく, 丁度識別なので当てはまりの良し悪しでモデルを疑う手段を失います. 外から信頼できる値が持ち込めるなら, それを使って過剰識別の余地を残す方が, 教科書的には筋が通っています.
最後に大きさを見比べておきます. \(\hat\sigma_\psi\) の標準誤差は 0.018 で, 推定値の 3.5% でした. 一方 表 4.10 では, 外生的に定めた値を動かすと \(\hat\sigma_\psi\) が 0.37 から 0.79 まで動きます. 標本誤差より一桁大きい幅です. そのためこの表は, 標準誤差と切り離さずに, 推定結果の一部として並べて報告するのが望ましいでしょう.
Mazzocco (2007) は逆に, 家計のオイラー方程式を使って家計内のコミットメントを検定します. Lise and Yamada (2019) が期内条件のみを使うのは, その検定が流動性制約と区別できないという問題を避けるためです.↩︎
R の
wooldridgeパッケージ (GPL-3) のwage2として配布されているものです.↩︎8点から20点まで変えても, 以下で報告する \(J\) 統計量は小数第2位まで動きません.↩︎
分位点がその例です. \(\tau\) 分位点 \(q_\tau\) の影響関数は \(\left(\tau - 1\left\{y_i \leq q_\tau\right\}\right) / f\left(q_\tau\right)\) と解析的に書けますが, 分位点での密度 \(f\left(q_\tau\right)\) をカーネル密度推定で求める必要があり, 結果がバンド幅の選び方に左右されます. また, この節の労働時間のように値が1点に集中していると (中央値はちょうど40時間です), その点で密度が存在せず, 分位点の漸近正規性そのものが崩れます.↩︎
逆に, 先行研究が点推定だけを報告していて, モーメントの標準誤差すら手に入らない場合には, \(\sigma_j\) が分からないので 命題 4.6 の上限も計算できず, 推定値の標準誤差は求まりません.↩︎
Andrews et al. (2017) はモーメント条件を \(\boldsymbol{g}\left(\boldsymbol{\theta}\right) = \hat{\boldsymbol{s}} - \boldsymbol{s}\left(\boldsymbol{\theta}\right)\) と「データ引くモデル」の向きで書くので, 定義は \(\Lambda = -\left(G^\top W G\right)^{-1} G^\top W\) と符号がつきます (Andrews et al. 2017, def. 1). この節では \(\boldsymbol{m}\left(\boldsymbol{\theta}\right) - \hat{\boldsymbol{m}}\) の向きで書いているので符号が消え, 同じ行列になります.↩︎