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）
异质性 Agent	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
因果机器学习 / 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

扰动法（二阶近似）

• 对 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
Wild bootstrap	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}$$

运行 六个渐进式规格（M1–M6）：

模型	单位 FE	时间 FE	协变量	基准×趋势	区域 FE	单位趋势
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，理想情况下 > 104，参见 Lee 等，2022）、排除限制、独立性。

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 诊断清单

检验	检验内容	通过标准
第一阶段 F	工具变量相关性	F > 104（Lee 等）或 > 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 / Shift-Share 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)

---

因果方法选择指南

场景	方法	主要库
含固定效应的 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
基于阈值的分配规则	Sharp / Fuzzy RD	rdrobust
斜率不连续	RKD	rdrobust（deriv=1）
N=1 个处理组单位、长面板	合成控制法	pysynth / 手动实现
合成控制 + 面板结构	合成 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 权重不可靠；明确报告拟合情况
VFI 内层循环未使用 Numba	使用 @njit 装饰
收入使用均匀网格	使用 Tauchen / Rouwenhorst 离散化
资产网格为线性网格	在借款约束附近采用对数 / 指数间距
未检查求解器是否收敛	检查 sol.success 和残差

---

调试清单

DID： 处理前时期的事件研究系数 ≈ 0；使用 Goodman-Bacon 分解检查 TWFE 权重

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)；Jacobian：np.allclose(J_analytic, J_fd, rtol=1e-4)