python-econ-computing

Python による経済数値計算

• 著者：Wenli Xu
• メール： [email protected]
• 2026-03-11

---

概要

Python によるマクロ経済モデリング（DSGE/HANK）、因果推論、データ分析のベストプラクティスです。基本原則：まずベクトル化し、ループは Numba で高速化し、コード構造を経済理論に沿わせます。

---

ライブラリ早見表

用途	推奨ライブラリ
数値計算の中核	numpy、scipy
ループの高速化	numba（@njit、@njit(parallel=True)）
経済学ツールキット	quantecon
HANK／シーケンス空間	sequence_jacobian（SSJ）
異質なエージェント	HARK
固定効果を含む線形モデル	pyfixest（pip install pyfixest）
DID／DD／DDD	diff-diff（pip install diff-diff）
IV／2SLS／GMM	linearmodels（または固定効果を含むパネル IV には pyfixest）
RD／RDD／RKD	rdrobust、rddensity、rdlocrand
合成コントロール法	pysynth、synth_control、sdid
マッチング	causalml、pymatch、econml
因果 ML／DML	econml、dowhy
データ操作	pandas、polars（大規模データセット）
可視化	matplotlib、seaborn

---

DSGE モデル

線形化と解法（Blanchard-Kahn）

import numpy as np
from scipy.linalg import ordqz

def solve_bk(A, B, n_fwd):
    """
    Solve linear DSGE: A E_t[x_{t+1}] = B x_t + C eps_t
    n_fwd: number of forward-looking variables
    Returns decision rule matrix P such that x_t = P x_{t-1} + ...
    """
    AA, BB, alpha, beta, Q, Z = ordqz(A, B, sort='ouc')
    n = A.shape[0]
    Z21 = Z[n - n_fwd:, :n - n_fwd]
    Z22 = Z[n - n_fwd:, n - n_fwd:]
    P = -np.linalg.solve(Z22, Z21)
    return P

摂動法（2 次近似）

• LQ 問題には quantecon.lqcontrol を使用します
• 高次の摂動法には perturbpy または手動実装を使用します
• 定常状態の求解には scipy.optimize.fsolve／root を使用します

---

HANK モデル

シーケンス空間ヤコビアン法（SSJ）

import sequence_jacobian as sj

# 1. Define steady-state blocks
@sj.simple
def household_ss(r, w, beta, sigma):
    # Return steady-state aggregates
    ...

# 2. Build DAG
model = sj.create_model([household_block, firm_block, market_clearing],
                         name='HANK')

# 3. Solve steady state
ss = model.solve_steady_state(calibration, unknowns, targets)

# 4. Compute Jacobians → solve transition dynamics
G = model.solve_jacobian(ss, unknowns, targets, T=300)

価値関数反復 — Numba による高速化

from numba import njit
import numpy as np

@njit
def vfi(V0, a_grid, y_grid, r, beta, sigma, tol=1e-8, max_iter=1000):
    """Heterogeneous agent VFI over asset grid × income grid"""
    n_a, n_y = len(a_grid), len(y_grid)
    V = V0.copy()
    policy = np.zeros((n_a, n_y))

    for it in range(max_iter):
        V_new = np.empty_like(V)
        for ia in range(n_a):
            for iy in range(n_y):
                best_val = -1e10
                best_a = 0
                for ia2 in range(n_a):
                    c = (1 + r) * a_grid[ia] + y_grid[iy] - a_grid[ia2]
                    if c <= 0:
                        continue
                    u = c ** (1 - sigma) / (1 - sigma)
                    val = u + beta * V[:, iy].mean()  # use transition matrix in practice
                    if val > best_val:
                        best_val = val
                        best_a = ia2
                V_new[ia, iy] = best_val
                policy[ia, iy] = a_grid[best_a]
        if np.max(np.abs(V_new - V)) < tol:
            break
        V = V_new
    return V, policy

分布の反復（Young 2010）

def iterate_distribution(policy_idx, trans_mat, dist0, T=500):
    """Iterate joint distribution to steady state given policy indices and income transition matrix"""
    dist = dist0.copy()
    n_a, n_y = dist.shape
    for _ in range(T):
        dist_new = np.zeros_like(dist)
        for iy in range(n_y):
            for iy2 in range(n_y):
                dist_new[policy_idx[:, iy], iy2] += dist[:, iy] * trans_mat[iy, iy2]
        dist = dist_new
    return dist

---

固定効果を含む線形モデル（pyfixest）

ルール：固定効果を含む OLS／Poisson／Logit には、すべて pyfixest を使用します。R の fixest と同じ構文です。

import pyfixest as pf

# OLS with unit + time FE, cluster-robust SEs
fit = pf.feols("y ~ treat_post | unit + year",
               data=df, vcov={"CRV1": "id"})
fit.summary()

# Multiple high-dimensional FE (Frisch-Waugh absorbed)
fit = pf.feols("y ~ x1 + x2 | unit + year + industry",
               data=df, vcov={"CRV1": "id"})

# Wild cluster bootstrap (few clusters, <50)
fit = pf.feols("y ~ treat_post | unit + year",
               data=df, vcov={"CRV1": "id"})
fit.wildboottest(param="treat_post", B=9999, seed=42)

# Event study via i() syntax
fit = pf.feols("y ~ i(rel_year, ref=-1) | unit + year",
               data=df, vcov={"CRV1": "id"})
pf.iplot(fit)  # event study plot

# Poisson (count / log-linear) with FE
fit_pois = pf.fepois("y ~ treat_post | unit + year",
                     data=df, vcov={"CRV1": "id"})

# Access results
fit.coef()           # coefficient estimates
fit.se()             # standard errors
fit.pvalue()         # p-values
fit.confint()        # confidence intervals
fit._N               # number of observations

pyfixest と statsmodels の比較

用途	使用するライブラリ
固定効果を含む任意の OLS／WLS	pyfixest
固定効果を含む Poisson／logit	pyfixest
ワイルドブートストラップ	pyfixest
時系列 ARIMA、VAR	statsmodels
固定効果なしの MLE／GLM	statsmodels

---

因果推論：DID／DD／DDD 法

ルール：DID、DD、DDD、または段階的な差の差デザインには、すべて diff-diff を使用し、以下の一般的な実証分析ワークフローに従います。

ソース：https://github.com/igerber/diff-diff | https://github.com/wenddymacro/A-General-Empirical-Workflow-for-DID

diff-diff 推定量のリファレンス

別名	クラス	使用場面
DiD	DifferenceInDifferences	基本的な 2×2 DiD
TWFE	TwoWayFixedEffects	標準的なパネル DiD
EventStudy	MultiPeriodDiD	動学的効果／イベントスタディ
CS	CallawaySantAnna	段階的導入、異質な効果
SA	SunAbraham	段階的導入、負の重みを回避
BJS	ImputationDiD	Borusyak らの補完アプローチ
SDiD	SyntheticDiD	合成 DiD
DDD	TripleDifference	三重差分

主な API パラメーター

from diff_diff import DiD, TWFE, EventStudy, CS, SA, BJS, DDD

# Common fit() arguments
results = estimator.fit(
    data,
    outcome='y',           # dependent variable
    treatment='treated',   # binary treatment indicator
    time='post',           # binary post-period (or period var for panel)
    unit='id',             # unit identifier (panel)
    covariates=['x1','x2'],# control variables
    absorb=['region'],     # high-dim fixed effects (within-transform)
    cluster='id',          # clustered standard errors
    robust=True,           # HC1 robust SEs
    inference='wild_bootstrap',  # for few clusters (<50)
    n_bootstrap=999,
)

# Results
results.att           # ATT estimate
results.se            # standard error
results.p_value
results.conf_int      # confidence interval tuple
results.print_summary()
results.to_dataframe()

DDD（三重差分）

from diff_diff import DDD

ddd = DDD()
results = ddd.fit(
    data,
    outcome='y',
    treatment='treated',
    time='post',
    third_diff='group_var',  # third differencing dimension
    cluster='id',
)

---

DID の実証分析ワークフロー（11 ステップ）

すべての DID／DD／DDD 論文または分析で、このワークフローに従ってください。

ステップ 1 — データと記述統計

• パネルを構築します。分析単位、期間、処置変数を定義します
• データソース、欠損値、ウィンザー化について記録します
• コホート構造（コホートごとの処置群の単位数、政策開始日）を記述します
• 表： 全サンプルの統計量と、t 検定を用いた処置群と対照群の比較
• 処置前の共変量バランス検定

共変量の種類：

種類	形式	目的
共変量	処置前、時間不変	平行トレンドを条件付ける
制御変数	ベースライン × 時間トレンド	残差の異質性を吸収する

ステップ 2 — 識別戦略

以下を述べ、根拠を示します。

1. 平行トレンド — 事前期間のイベントスタディで部分的に検証可能です
2. SUTVA／干渉なし — 単位間の波及がないことを論じます
3. 予期なし — 事前期間の係数はゼロであるべきです

政策の割り当てメカニズムを記録し、外生性の根拠として政策文書を引用します。

ステップ 3 — ベースライン回帰

モデル： $$Y_{it} = \alpha + \beta(\text{Treat}i \times \text{Post}{it}) + \gamma W_i + \delta(Z_i^{pre} \times t) + \mu_i + \lambda_t + \varepsilon_{it}$$

6つの段階的な仕様（M1–M6）を推定します。

モデル	単位固定効果	時間固定効果	共変量	ベースライン×トレンド	地域固定効果	単位トレンド
M1	✓	—	—	—	—	—
M2	✓	✓	—	—	—	—
M3	✓	✓	✓	✓	—	—
M4	✓	✓	✓	✓	✓	—
M5	✓	✓	✓	✓	✓	✓
M6	✓	✓	✓	✓	産業×年	✓

M1→M6で係数が安定していれば、識別が支持されます。Oster（2019）の δ*（選択バイアス比）を報告します。|δ*| > 1 は基本的な頑健性、|δ*| > 2 は強い頑健性を示します。

import pyfixest as pf

# M1: unit FE only
pf.feols("y ~ treat_post | unit", data=df, vcov={"CRV1": "id"}).summary()

# M2: unit + time FE
pf.feols("y ~ treat_post | unit + year", data=df, vcov={"CRV1": "id"}).summary()

# M3–M4: add covariates and baseline×trend
pf.feols("y ~ treat_post + x1 + x2 + baseline:year | unit + year",
         data=df, vcov={"CRV1": "id"}).summary()

# M5: + regional FE
pf.feols("y ~ treat_post + x1 + x2 + baseline:year | unit + year + region",
         data=df, vcov={"CRV1": "id"}).summary()

# M6: + industry×year FE
pf.feols("y ~ treat_post + x1 + x2 + baseline:year | unit + industry^year",
         data=df, vcov={"CRV1": "id"}).summary()

ステップ 4 — 平行トレンド：イベントスタディのプロット

from diff_diff import EventStudy

es = EventStudy()
res = es.fit(data, outcome='y', treatment='treated',
             unit='id', time='year', base_period=-1,
             cluster='id')
res.plot()  # shows pre/post coefficients with CIs

プロットの基準：

• Y軸：「係数（ATT）」、X軸：「政策からの相対年」
• 0の基準線を表示し、基準期間 k = −1 は省略します
• 本文にはM6を示し、M1–M5は付録に示します
• 事前期間の係数（k ∈ {−4,−3,−2}）は統計的に有意でないことが望まれます

段階的導入のデザインでは、すべての推定量を実行します。

from diff_diff import CS, SA, BJS

# Callaway & Sant'Anna (2021)
cs = CS().fit(data, outcome='y', unit='id', time='year',
              cohort='treat_year', control_group='never_treated', cluster='id')

# Sun & Abraham (2021)
sa = SA().fit(data, outcome='y', unit='id', time='year',
              cohort='treat_year', cluster='id')

# Borusyak et al. (2024) imputation
bjs = BJS().fit(data, outcome='y', unit='id', time='year',
                cohort='treat_year', horizons=range(5), cluster='id')

ステップ 5 — 平行トレンドの感度分析：HonestDiD

Rambachan-Roth（2023）の境界を用いて、平行トレンドの違反に対する頑健性を定量化します。

# After event study, extract pre/post coefficients and covariance
# Pass to HonestDiD (R package via rpy2, or use diff-diff's built-in honest DiD)
from diff_diff import EventStudy

es = EventStudy()
res = es.fit(data, ..., honest_did=True,
             sensitivity_constraint='smoothness')
res.plot_honest_did()  # shows identified set under relaxed PT assumption

ステップ 6 — 代替説明を排除する

脅威	検定
空間的波及	地理的プラセボ、処置対象単位からの距離別の効果
予期効果	事前期間のイベントスタディ係数 ≈ 0
政策の重複	同時期の政策を除外するか、統制します

ステップ 7 — 頑健性検証

from diff_diff import DiD
from diff_diff.diagnostics import PlaceboTest, GoodmanBaconDecomposition

# Goodman-Bacon decomposition (TWFE bias diagnosis)
gb = GoodmanBaconDecomposition().fit(data, outcome='y', treatment='treated',
                                      unit='id', time='year')
gb.plot()

# Placebo tests
placebo = PlaceboTest(method='fake_timing').fit(data, ...)
placebo = PlaceboTest(method='permutation', n_permutations=500).fit(data, ...)

# Subsample / specification robustness
for subsample_mask in subsamples:
    res = CS().fit(data[subsample_mask], ...)

標準的な頑健性検証項目：

• プラセボ／偽のタイミング
• 改革前期間に対する反証テスト
• 代替のアウトカム変数
• 事前に定められた特性によるサブサンプル分割
• ランダム化推論

ステップ 8 — 異質な処置効果

# Triple difference for effect heterogeneity by subgroup Z
from diff_diff import DDD

ddd = DDD().fit(data, outcome='y', treatment='treated',
                time='post', third_diff='high_exposure', cluster='id')

# Interaction-based heterogeneity in TWFE
from diff_diff import TWFE
res = TWFE().fit(data, outcome='y',
                 treatment='treat_post',
                 interactions=['treat_post:firm_size'],
                 absorb=['id', 'year'], cluster='id')

ステップ 9 — メカニズム分析

1. アウトカムの段階： 直接的 → 中間的 → 最終的なアウトカム
2. 媒介分析： 提案されたメカニズムを統制変数として含め、ATTを含める場合と含めない場合で比較します
3. メカニズムとしての異質性： 初期条件別にサブグループのイベントスタディを行います

ステップ 10 — 福祉と政策への示唆

• ATTを集計して政策レベルの影響を算出します
• 費用便益比
• 分配効果（受益者と損失を被る人）

ステップ 11 — ワークフロー全体の概要

1. Data & balance
2. Identification assumptions
3. Baseline specs M1–M6 + Oster δ*
4. Event study (TWFE + CS + SA + BJS)
5. HonestDiD sensitivity
6. Alternative explanations
7. Robustness battery
8. Heterogeneous effects (DDD / interactions)
9. Mechanisms
10. Welfare implications

---

因果推論：手法の選択

What is your identification strategy?
├── Policy/treatment with parallel trends → DID (see above)
├── Exogenous instrument for endogenous X → IV
├── Discontinuity in assignment rule → RD / RKD
├── Control units that can be reweighted → Synthetic Control
├── Selection on observables → Matching / IPW
└── High-dimensional / ML setting → DML / Causal Forest

---

操作変数法（IV／2SLS／GMM）

ライブラリ： linearmodels（パネルIVではstatsmodelsより推奨）

主な仮定： 関連性（F > 10、理想的にはLee et al. 2022に従い > 104）、排除制約、独立性。

from linearmodels.iv import IV2SLS, IVGMM, IVLIML

# Basic 2SLS: y ~ X_exog + [X_endog ~ Z_instruments]
res = IV2SLS(dependent=y,
             exog=X_exog,        # included exogenous (+ constant)
             endog=X_endog,      # endogenous regressors
             instruments=Z).fit(cov_type='robust')

# Panel IV with fixed effects
from linearmodels import PanelOLS, BetweenOLS
from linearmodels.iv import IV2SLS
# absorb FE first (within transform), then IV on residuals
# or use linearmodels.panel with IV support

# GMM (efficient with heteroskedasticity)
res = IVGMM(y, X_exog, X_endog, Z).fit(cov_type='robust')

# LIML (less biased with weak instruments)
res = IVLIML(y, X_exog, X_endog, Z).fit(cov_type='robust')

# Key diagnostics
print(res.first_stage)          # first-stage F-statistic
print(res.wooldridge_score)     # endogeneity test (H0: OLS consistent)
print(res.sargan)               # overidentification test (J-stat, requires overid)

IV診断チェックリスト

検定	確認する内容	合格基準
第1段階F	操作変数の関連性	F > 104（Lee et al.）または > 10（経験則）
Cragg-Donald／Kleibergen-Paap	弱い操作変数（内生変数が複数の場合）	Stock-Yogoの臨界値を上回る
Sargan-Hansen J検定	過剰識別（排除制約）	p > 0.1（妥当性を棄却できない）
Hausman／Wooldridge	Xの内生性	p < 0.05 → IVが必要
誘導形	操作変数がアウトカムに影響するか	有意であること

# Anderson-Rubin confidence set (robust to weak instruments)
from linearmodels.iv import IV2SLS
res = IV2SLS(y, X_exog, X_endog, Z).fit(cov_type='robust')
print(res.anderson_rubin)  # AR test, valid even with weak instruments

# Conley spatial HAC SEs (geographic instruments)
res = IV2SLS(y, X_exog, X_endog, Z).fit(cov_type='kernel', bandwidth=5)

Bartik／シフトシェアIV

# Bartik instrument: Z_i = sum_k s_{ik} * g_k
# s_{ik}: industry share of unit i; g_k: national industry growth
import numpy as np

def bartik_instrument(shares, growth):
    """
    shares: (n_units, n_industries)
    growth: (n_industries,)
    returns: (n_units,) Bartik instrument
    """
    return shares @ growth

---

回帰不連続デザイン（RD／RKD／ファジーRD）

ライブラリ： rdrobust（RパッケージのPython移植版）

主な仮定： 単位がカットオフ付近のランニング変数を正確に操作できないこと。

from rdrobust import rdrobust, rdbwselect, rdplot

# Sharp RD
res = rdrobust(y, x, c=cutoff)          # default: MSE-optimal bandwidth, local linear
res = rdrobust(y, x, c=0,
               kernel='triangular',     # triangular (default) / uniform / epanechnikov
               bwselect='mserd',        # MSE-optimal (default); 'cerrd' for coverage
               vce='hc1',              # robust SEs
               cluster=cluster_var)
print(res)

# Fuzzy RD (instrument = 1[x >= c])
res_fuzzy = rdrobust(y, x, c=0,
                     fuzzy=treatment_var)  # IV-style, estimates LATE

# Bandwidth selection
bw = rdbwselect(y, x, c=0, bwselect='mserd')
print(bw.bws)    # optimal bandwidth

# Visualization
rdplot(y, x, c=0)   # binned scatter with polynomial fit

RD診断チェックリスト

from rddensity import rddensity
from rdrobust import rdrobust

# 1. McCrary density test (H0: no manipulation at cutoff)
den = rddensity(x, c=cutoff)
print(den.test)   # p > 0.05: no evidence of manipulation

# 2. Covariate smoothness (placebo on pre-determined covariates)
for cov in baseline_covariates:
    res = rdrobust(cov, x, c=cutoff)
    print(f'{cov}: {res.coef[0]:.3f} (p={res.pv[2]:.3f})')  # should be insignificant

# 3. Placebo cutoffs (should find no effect at fake cutoffs)
for fake_c in [cutoff - 0.5, cutoff + 0.5]:
    res = rdrobust(y, x, c=fake_c)
    print(f'Placebo c={fake_c}: {res.coef[0]:.3f}')

# 4. Sensitivity to bandwidth
for h in [bw_opt * 0.5, bw_opt * 0.75, bw_opt, bw_opt * 1.25, bw_opt * 1.5]:
    res = rdrobust(y, x, c=cutoff, h=h)
    print(f'h={h:.2f}: {res.coef[0]:.3f}')

# 5. Donut hole (exclude units very close to cutoff)
mask = np.abs(x - cutoff) > donut_radius
res_donut = rdrobust(y[mask], x[mask], c=cutoff)

回帰キンクデザイン（RKD）

# RKD: identifies effect via kink (slope discontinuity) rather than level jump
res_rkd = rdrobust(y, x, c=cutoff, deriv=1)  # deriv=1 estimates slope discontinuity

---

合成コントロール法

使用する場面： 処置ユニットが少数（多くの場合 N=1）で、処置前のパネルデータが長く、明確なコントロール群がない場合。

ライブラリ： pysynth、synth_control（pip）、または scipy.optimize を使った手動実装。

# --- Option 1: pysynth ---
from pysynth import Synth

sc = Synth()
sc.fit(
    dataprep={
        'foo_table': df,
        'predictors': ['gdp', 'trade', 'invest'],
        'predictors_op': 'mean',
        'time_predictors_prior': list(range(1980, 1990)),
        'special_predictors': [('gdp', [1985, 1988], 'mean')],
        'dependent': 'gdp',
        'unit_variable': 'country',
        'time_variable': 'year',
        'treatment_identifier': 'basque',
        'controls_identifier': control_countries,
        'time_optimize_ssr': list(range(1960, 1990)),
        'time_plot': list(range(1960, 1998)),
    }
)
sc.plot(['trends', 'weights', 'gaps'])

# --- Option 2: manual (scipy) ---
from scipy.optimize import minimize
import numpy as np

def synth_loss(w, Y_pre_control, Y_pre_treated):
    """Minimize pre-treatment fit: ||Y_treated - Y_control @ w||^2"""
    return np.sum((Y_pre_treated - Y_pre_control @ w) ** 2)

n_controls = Y_pre_control.shape[1]
w0 = np.ones(n_controls) / n_controls
constraints = [{'type': 'eq', 'fun': lambda w: w.sum() - 1}]
bounds = [(0, 1)] * n_controls

res = minimize(synth_loss, w0,
               args=(Y_pre_control, Y_pre_treated),
               method='SLSQP',
               bounds=bounds,
               constraints=constraints)
w_opt = res.x
Y_synth = Y_post_control @ w_opt
gap = Y_post_treated - Y_synth

合成コントロール法の診断

# Pre-treatment fit (RMSPE)
rmspe_pre = np.sqrt(np.mean((Y_pre_treated - Y_pre_control @ w_opt)**2))

# Placebo tests: apply SC to each control unit, compute distribution of gaps
placebo_gaps = []
for ctrl in control_units:
    Y_treated_placebo = Y_pre[:, ctrl_idx]
    Y_control_placebo = np.delete(Y_pre, ctrl_idx, axis=1)
    # ... fit and store gap
    placebo_gaps.append(gap_placebo)

# Ratio: treated RMSPE_post / RMSPE_pre vs. controls (Abadie et al. 2010)
ratio_treated = rmspe_post / rmspe_pre
# Inference: fraction of placebos with ratio >= ratio_treated → p-value

# In-time placebo: apply SC using period before actual treatment as fake treatment
# In-space placebo: already done above

合成 DiD（SDiD）

# Combines SC weights with DiD — robust to both parallel trends violations and
# imperfect pre-treatment fit
from diff_diff import SDiD

sdid = SDiD()
res = sdid.fit(data, outcome='y', treatment='treated',
               unit='id', time='year', cluster='id')
res.print_summary()

---

マッチングと再重み付け

使用する場面： 観測可能な共変量に基づく選択があり、ベースラインの共変量データが豊富な場合。

推定対象： ATT（処置群とマッチングしたコントロール群の比較）、ATE（母集団平均）。

傾向スコア法

from sklearn.linear_model import LogisticRegression
from sklearn.preprocessing import StandardScaler
import numpy as np

# 1. Estimate propensity score
X_scaled = StandardScaler().fit_transform(X_covariates)
ps_model = LogisticRegression(C=1.0, max_iter=1000)
ps_model.fit(X_scaled, treatment)
p_score = ps_model.predict_proba(X_scaled)[:, 1]

# 2. Check overlap / common support
import matplotlib.pyplot as plt
plt.hist(p_score[treatment==1], alpha=0.5, label='Treated', bins=30)
plt.hist(p_score[treatment==0], alpha=0.5, label='Control', bins=30)
plt.legend(); plt.xlabel('Propensity Score')
# Trim tails: drop obs with p_score outside [0.05, 0.95]
mask = (p_score >= 0.05) & (p_score <= 0.95)

IPW / AIPW（二重ロバスト）

from econml.dr import LinearDRLearner
from sklearn.linear_model import LassoCV, LogisticRegressionCV

# Doubly robust (AIPW) — consistent if either outcome or propensity model correct
dr = LinearDRLearner(
    model_regression=LassoCV(),      # outcome model
    model_propensity=LogisticRegressionCV(),  # propensity model
    featurizer=None
)
dr.fit(Y, T, X=X_het, W=X_controls)  # X: effect modifiers, W: controls
ate = dr.ate(X_het)
print(dr.ate_interval(X_het))         # confidence interval

最近傍マッチング

from causalml.match import NearestNeighborMatch
from causalml.propensity import ElasticNetPropensityModel

# Propensity score matching
pm = ElasticNetPropensityModel()
ps = pm.fit_predict(X_covariates, treatment)

matcher = NearestNeighborMatch(replace=False, ratio=1, random_state=42)
matched = matcher.match(data=df, treatment_col='treated', score_cols=['ps'])

# ATT on matched sample
att = matched[matched.treated==1]['y'].mean() - matched[matched.treated==0]['y'].mean()

# OR: Mahalanobis distance matching (better for low-dimensional X)
from pymatch.Matcher import Matcher
m = Matcher(test=df[df.treated==1], control=df[df.treated==0],
            yvar='y', exclude=['id'])
m.fit_scores(balance=True, nmodels=10)
m.predict_scores()
m.match(method='min', nmatches=1, threshold=0.001)
m.assess_balance(actual=True)

共変量バランスの診断

# Standardized mean differences (SMD) before/after matching
def smd(x_treat, x_control):
    return (x_treat.mean() - x_control.mean()) / np.sqrt(
        (x_treat.var() + x_control.var()) / 2
    )

for col in covariates:
    before = smd(df[df.treated==1][col], df[df.treated==0][col])
    after  = smd(matched[matched.treated==1][col], matched[matched.treated==0][col])
    print(f'{col}: SMD before={before:.3f}, after={after:.3f}')
# Target: |SMD| < 0.1 after matching

# Love plot
import matplotlib.pyplot as plt
smds_before = [...]
smds_after  = [...]
plt.scatter(smds_before, covariates, label='Before', marker='o')
plt.scatter(smds_after,  covariates, label='After',  marker='s')
plt.axvline(0, color='k', lw=0.5); plt.axvline(0.1, color='r', ls='--')
plt.legend(); plt.xlabel('Standardized Mean Difference')

エントロピーバランシング

# Reweight controls to exactly match treated means (and optionally variances)
# Install: pip install ebal
from ebal import ebal

# Balances moments of X_controls exactly — no propensity model needed
weights = ebal(X_control=X[treatment==0],
               X_treated=X[treatment==1],
               moments=1)   # 1=means, 2=means+variances

# Use weights in weighted regression
import pyfixest as pf
df["w"] = np.where(treatment == 1, 1.0, weights)
res = pf.feols("y ~ treated", data=df, weights="w", vcov={"CRV1": "id"})

---

二重機械学習（DML）と因果フォレスト

使用する場面： 高次元のコントロール変数がある場合、処置効果に異質性がある場合、関数形に柔軟性を持たせたい場合。

from econml.dml import LinearDML, CausalForestDML, NonParamDML
from econml.dr import ForestDRLearner
from sklearn.ensemble import GradientBoostingRegressor, GradientBoostingClassifier
from sklearn.linear_model import LassoCV, LogisticRegressionCV

# --- Linear DML (Partially Linear Robinson model) ---
dml = LinearDML(
    model_y=LassoCV(),              # outcome residualization
    model_t=LassoCV(),              # treatment residualization
    discrete_treatment=False,
    cv=5,
)
dml.fit(Y, T, X=X_het, W=X_controls)
print(dml.ate(), dml.ate_interval())

# --- Causal Forest (nonparametric CATE) ---
cf = CausalForestDML(
    model_y=GradientBoostingRegressor(),
    model_t=GradientBoostingRegressor(),
    n_estimators=1000,
    min_samples_leaf=5,
    max_depth=5,
    discrete_treatment=False,
    cv=5,
)
cf.fit(Y, T, X=X_het, W=X_controls)

# Heterogeneous effects
tau_hat = cf.effect(X_het)            # CATE for each unit
lb, ub = cf.effect_interval(X_het)   # 95% CI

# Feature importance for heterogeneity
cf.feature_importances_  # which X drives heterogeneity

# Best linear predictor of CATE
blp = cf.ate_inference(X_het)
blp.summary_frame()

# --- IV + DML (DRIV for endogenous treatment) ---
from econml.iv.dr import LinearDRIV
driv = LinearDRIV(
    model_y_xw=LassoCV(),
    model_t_xw=LassoCV(),
    model_z=LogisticRegressionCV(),  # instrument model
    discrete_instrument=True,
)
driv.fit(Y, T, Z=Z_instrument, X=X_het, W=X_controls)

---

因果推論手法の選択ガイド

設定	手法	主なライブラリ
FE を含む OLS / WLS / Poisson	線形モデル	pyfixest
パネルデータ＋政策ショック、平行トレンド	DID / TWFE / CS / SA	diff-diff + pyfixest
段階的な導入	CS、SA、BJS	diff-diff
外生的な操作変数	2SLS / GMM / LIML	linearmodels（またはパネル IV には pyfixest）
操作変数の弱さが懸念される場合	AR 信頼集合、LIML	linearmodels
カットオフによる割り当てルール	シャープ / ファジー RD	rdrobust
傾きの不連続	RKD	rdrobust（deriv=1）
処置群 N=1、長期パネル	合成コントロール法	pysynth / 手動
SC＋パネル構造	合成 DiD	diff-diff（SDiD）
観測可能な変数に基づく選択	PSM / IPW / EB	causalml、ebal
高次元のコントロール変数、二値 T	AIPW / DR-Learner	econml
異質な効果	因果フォレスト	econml
内生的な T＋異質性	DRIV	econml

---

一般的な数値計算パターン

根の探索

from scipy.optimize import brentq, root

# Scalar: prefer brentq (robust)
r_star = brentq(lambda r: asset_market_clearing(r, params), -0.05, 0.1)

# Multivariate
sol = root(equilibrium_system, x0=initial_guess, method='hybr', tol=1e-10)

所得過程の離散化

import quantecon as qe

# Tauchen: AR(1) log y' = rho log y + sigma_e * eps
mc = qe.tauchen(rho, sigma_e, n=7)
y_grid = np.exp(mc.state_values)
trans_mat = mc.P

# Rouwenhorst (better for high persistence)
mc = qe.rouwenhorst(n=7, rho=rho, sigma=sigma_e)

パフォーマンスの優先順位

# 1. Vectorize with numpy first
# 2. Must loop → @njit
# 3. Parallelizable outer loop → @njit(parallel=True) + prange
# 4. Sparse structure → scipy.sparse

from numba import njit, prange

@njit(parallel=True)
def parallel_vfi(V, a_grid, y_grid, beta, sigma):
    n_a = len(a_grid)
    V_new = np.empty_like(V)
    for ia in prange(n_a):
        ...
    return V_new

---

よくある間違い

間違い	正しい方法
段階的な処置に TWFE を使用する	負の重みによるバイアスを避けるため、CS / SA / BJS を使用する
クラスター化標準誤差なしで DID を行う	diff_diff で cluster='id' を指定する
クラスター数が少ない（<50）	diff-diff で inference='wild_bootstrap' を指定する
IV：第一段階の F を確認しない	常に res.first_stage を出力する。F > 104 が望ましい
IV：過剰識別の場合に J 検定の p < 0.05	操作変数が無効である可能性が高いため、除外制約を再検討する
RD：単一のバンド幅のみを選択する	複数のバンド幅にわたる頑健性を示す
RD：カットオフで密度を検定しない	必ず McCrary / rddensity 検定を実行する
マッチング：バランスを確認しない	マッチング前後の SMD を報告する。目標は |SMD| < 0.1
マッチング：共通サポートを無視する	p スコアが [0.05, 0.95] の範囲外にあるサンプルを除外する
SC：処置前の適合度が低い	RMSPE_pre が高い → SC の重みの信頼性が低い。適合度を明示的に報告する
Numba を使わず VFI の内側ループを実行する	@njit を付ける
所得に一様グリッドを使う	Tauchen / Rouwenhorst による離散化を使う
資産に線形グリッドを使う	借入制約の近くでは対数／指数間隔を使う
ソルバーの収束を確認しない	sol.success と残差を確認する

---

デバッグチェックリスト

DID： 処置前期間のイベントスタディ係数 ≈ 0。TWFE の重みを確認するには Goodman-Bacon 分解を実施する

IV： 第一段階の F > 104、誘導形が有意、J 検定の p > 0.1（過剰識別の場合）。操作変数が弱い場合は AR 信頼集合を使用する

RD： 密度検定の p > 0.05、カットオフで共変量が滑らか、バンド幅の選択に対して頑健

SC： 処置前 RMSPE が小さい。推論にはプラセボ RMSPE 比（処置後／処置前）を使用する

マッチング： マッチング後の |SMD| < 0.1。Love プロット。共通サポートの重なりを確認する

DSGE/HANK： すべての市場清算残差が < 1e-8。VFI：max|V_{n+1} - V_n| をプロットする。分布：assert np.isclose(dist.sum(), 1.0)。ヤコビアン：np.allclose(J_analytic, J_fd, rtol=1e-4)