双重机器学习之Python 实战：从理论到代码一文通关

本人专门研究经管和计算机交叉领域，欢迎大家评论关注,完整代码在文章末尾！

摘要：双重机器学习（Double Machine Learning, DML）是近年来经管实证研究中最热门的因果推断方法之一。本文从直觉出发，系统介绍 DML 的理论框架、核心步骤，并提供完整的 Python 代码实现（含模拟数据、手动实现、EconML 库调用、蒙特卡洛验证），适合经管、社科、商业分析方向的同学快速上手。

一、为什么需要双重机器学习？

在经管实证研究中，我们经常面临一个核心问题：估计某个政策/干预/变量对结果的因果效应（处理效应）。

传统的做法是 OLS 回归，但当控制变量与处理变量、结果变量之间存在复杂的非线性关系时，线性模型的函数形式设定错误（misspecification）会导致遗漏变量偏差（omitted variable bias），估计结果有偏。

一个自然的想法是：能不能用机器学习来拟合这些复杂的非线性关系？

答案是可以，但直接把 ML 塞进回归框架会产生正则化偏差（regularization bias）。Chernozhukov et al. (2018) 提出的双重机器学习（Double/Debiased Machine Learning）框架，通过巧妙的"正交化 + 交叉拟合"策略，完美解决了这个问题。

DML 的核心优势

二、DML 理论框架

2.1 部分线性模型（Partially Linear Model）

DML 最经典的设定是部分线性模型：

$$ Y = \\theta_0 D + g_0(X) + \\epsilon, \\quad E[\\epsilon|D, X] = 0 $$

$$ D = m_0(X) + \\eta, \\quad E[\\eta|X] = 0 $$

其中：

• $Y$ 是结果变量
• $D$ 是处理变量（我们关心的因果变量）
• $X$ 是高维控制变量（协变量）
• $\\theta_0$ 是我们要估计的处理效应
• $g_0(X)$ 是控制变量对 $Y$ 的未知函数（可以是任意非线性形式）
• $m_0(X)$ 是控制变量对 $D$ 的未知函数（产生内生性/选择偏差）

2.2 Frisch-Waugh-Lovell 定理的 ML 推广

DML 的核心思想可以理解为 FWL 定理的机器学习版本：

$$ \\hat{\\theta} = \\frac{\\sum_i \\tilde{D}_i \\tilde{Y}_i}{\\sum_i \\tilde{D}_i^2} $$

2.3 为什么需要"双重"和"交叉拟合"？

"双重"（Double） 指的是对 $Y$ 和 $D$ 都进行正交化（去偏），这保证了即使一阶段的 ML 估计有偏差，对最终 $\\theta$ 的影响也是二阶小量（Neyman 正交性）。

"交叉拟合"（Cross-Fitting） 类似于交叉验证。将样本分为 $K$ 折：

• 用第 $k$ 折之外的数据训练 ML 模型
• 用第 $k$ 折的数据做预测

这避免了 ML 的过拟合偏差，保证了 $\\hat{\\theta}$ 的 $\\sqrt{n}$ 收敛速率和渐近正态性。

三、数据生成过程（DGP）

为了验证 DML 的效果，我们构造如下模拟数据：

def generate_dml_data(n=2000, p=10, seed=42):
np.random.seed(seed)
TRUE_THETA = 2.0 # 真实处理效应

X = np.random.normal(0, 1, size=(n, p))

# 非线性干扰函数 g(X)
g_X = (2 * np.sin(X[:, 0] * X[:, 1])
+ X[:, 2] ** 2
– 1.5 * np.abs(X[:, 3])
+ 0.8 * X[:, 4] * X[:, 5]
+ np.exp(0.3 * X[:, 6]))

# 非线性倾向函数 m(X)
m_X = (0.8 * X[:, 0]
+ 0.5 * X[:, 1] ** 2
– 0.6 * X[:, 2]
+ 0.4 * np.sin(X[:, 3]))

D = m_X + np.random.normal(0, 1, n)
Y = TRUE_THETA * D + g_X + np.random.normal(0, 1, n)

return X, D, Y, TRUE_THETA

关键设计：$g_0(X)$ 和 $m_0(X)$ 包含交互项、平方项、三角函数、指数函数等强非线性成分，线性模型无法正确捕捉，因此 OLS 估计会有偏。

四、完整代码实现

4.1 手动实现 DML

这是最核心的部分，三步走：

from sklearn.ensemble import GradientBoostingRegressor
from sklearn.model_selection import cross_val_predict, KFold

def dml_manual(X, D, Y, n_splits=5, seed=42):
kf = KFold(n_splits=n_splits, shuffle=True, random_state=seed)

# Step 1 & 2: 交叉拟合获取残差
ml_y = GradientBoostingRegressor(n_estimators=200, max_depth=4,
learning_rate=0.05, random_state=seed)
ml_d = GradientBoostingRegressor(n_estimators=200, max_depth=4,
learning_rate=0.05, random_state=seed)

Y_hat = cross_val_predict(ml_y, X, Y, cv=kf)
D_hat = cross_val_predict(ml_d, X, D, cv=kf)

Y_tilde = Y – Y_hat # 结果残差
D_tilde = D – D_hat # 处理残差

# Step 3: 残差回归
theta_hat = np.sum(D_tilde * Y_tilde) / np.sum(D_tilde ** 2)

# 异方差稳健标准误
residuals = Y_tilde – theta_hat * D_tilde
V = np.mean(D_tilde**2 * residuals**2) / (np.mean(D_tilde**2)**2)
se = np.sqrt(V / len(Y))

return theta_hat, se

代码要点：

• cross_val_predict 自动完成交叉拟合，避免手写循环
• 标准误使用异方差稳健（Heteroscedasticity-robust）公式
• ML 学习器可任意替换（随机森林、Lasso、XGBoost 等）

4.2 使用 EconML 库（推荐生产环境）

微软开源的 econml 库封装了完整的 DML 实现：

from econml.dml import LinearDML

dml = LinearDML(
model_y=GradientBoostingRegressor(n_estimators=200, max_depth=4),
model_t=GradientBoostingRegressor(n_estimators=200, max_depth=4),
cv=5,
random_state=42,
)
dml.fit(Y, D.reshape(-1, 1), X=None, W=X)

# 获取估计结果
theta = dml.const_marginal_effect()
inference = dml.const_marginal_effect_inference()
print(inference.summary())

注意：EconML 中 X 参数是异质性效应的特征（如果估计 CATE），W 是控制变量。对于 ATE 估计，传入 X=None, W=控制变量。

4.3 不同 ML 学习器对比

DML 的一大灵活性在于一阶段 ML 学习器的选择。以下对比三种常用学习器：

实验结论：当 DGP 包含强非线性时，GBRT 和随机森林的 DML 估计接近真实值，而 Lasso 的估计仍有偏差。

五、蒙特卡洛模拟验证

为了验证 DML 的无偏性，我们进行 200 次蒙特卡洛模拟：

# 每次模拟：
# 1. 重新生成数据
# 2. 分别用 OLS 和 DML 估计 theta
# 3. 收集所有估计值，画分布图

典型结果：

可以看到：

• OLS 估计系统性偏高（因为无法处理非线性干扰）
• DML 估计以真实值为中心，且近似正态分布

六、实际应用建议

6.1 什么时候用 DML？

适合以下场景：

• 你关心的是某个特定变量的因果效应（不是预测）
• 有大量控制变量，且控制变量与处理/结果之间可能存在非线性关系
• 传统 OLS 的函数形式设定可能不正确

6.2 实操 Checklist

6.3 常见误区

七、核心参考文献

八、完整代码获取

本文所有代码已整理为完整可运行的 Python 脚本，包含数据生成、手动 DML、EconML 调用、学习器对比、蒙特卡洛模拟和可视化等模块，可以直接在 Jupyter Notebook 或命令行运行。

运行环境要求：

Python >= 3.8
numpy, pandas, scikit-learn, matplotlib, scipy
econml (可选，用于 EconML 模块)

安装命令：

pip install numpy pandas scikit-learn matplotlib scipy econml

"""
双重机器学习（Double Machine Learning, DML）完整实战代码
============================================================
基于 Chernozhukov et al. (2018) 的理论框架
适用场景：因果推断中处理效应的无偏估计（高维控制变量）

依赖安装：
pip install numpy pandas scikit-learn econml matplotlib statsmodels
"""

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.ensemble import GradientBoostingRegressor, GradientBoostingClassifier
from sklearn.ensemble import RandomForestRegressor, RandomForestClassifier
from sklearn.linear_model import LassoCV
from sklearn.model_selection import cross_val_predict, KFold
from scipy import stats
import warnings
warnings.filterwarnings("ignore")

# ============================================================
# 第一部分：模拟数据生成（DGP）
# ============================================================
def generate_dml_data(n=2000, p=10, seed=42):
"""
生成用于双重机器学习的模拟数据

数据生成过程（DGP）：
Y = theta * D + g(X) + epsilon
D = m(X) + eta

其中：
– Y: 结果变量
– D: 处理变量（核心解释变量）
– X: 高维控制变量
– theta: 真实处理效应（我们要估计的参数）
– g(X): 控制变量对Y的非线性影响
– m(X): 控制变量对D的非线性影响（内生性来源）
"""
np.random.seed(seed)

# 真实处理效应
TRUE_THETA = 2.0

# 生成高维控制变量 X
X = np.random.normal(0, 1, size=(n, p))

# 非线性干扰函数 g(X) —— 控制变量对 Y 的复杂影响
g_X = (
2 * np.sin(X[:, 0] * X[:, 1])
+ X[:, 2] ** 2
– 1.5 * np.abs(X[:, 3])
+ 0.8 * X[:, 4] * X[:, 5]
+ np.exp(0.3 * X[:, 6])
)

# 非线性倾向函数 m(X) —— 控制变量对 D 的影响（产生内生性/选择偏差）
m_X = (
0.8 * X[:, 0]
+ 0.5 * X[:, 1] ** 2
– 0.6 * X[:, 2]
+ 0.4 * np.sin(X[:, 3])
+ 0.3 * X[:, 4]
)

# 误差项
epsilon = np.random.normal(0, 1, n)
eta = np.random.normal(0, 1, n)

# 生成处理变量 D（连续型）
D = m_X + eta

# 生成结果变量 Y
Y = TRUE_THETA * D + g_X + epsilon

# 整理为 DataFrame
col_names = [f"X{i+1}" for i in range(p)]
df = pd.DataFrame(X, columns=col_names)
df["D"] = D
df["Y"] = Y

print("=" * 60)
print("数据生成完成")
print(f" 样本量: {n}, 控制变量维度: {p}")
print(f" 真实处理效应 θ = {TRUE_THETA}")
print(f" Y 均值: {Y.mean():.3f}, D 均值: {D.mean():.3f}")
print("=" * 60)

return df, TRUE_THETA

# ============================================================
# 第二部分：朴素 OLS 估计（有偏，用于对照）
# ============================================================
def naive_ols(df):
"""朴素OLS回归（不控制非线性关系，结果有偏）"""
from sklearn.linear_model import LinearRegression

X_cols = [c for c in df.columns if c.startswith("X")]
X_ols = df[X_cols + ["D"]].values
Y_ols = df["Y"].values

reg = LinearRegression().fit(X_ols, Y_ols)
theta_ols = reg.coef_[-1] # D 的系数

return theta_ols

# ============================================================
# 第三部分：手动实现 DML（Partially Linear Model）
# ============================================================
def dml_manual(df, n_splits=5, seed=42):
"""
手动实现双重机器学习（DML）

核心步骤：
Step 1: 用 ML 模型拟合 Y ~ X，得到残差 Y_tilde = Y – E[Y|X]
Step 2: 用 ML 模型拟合 D ~ X，得到残差 D_tilde = D – E[D|X]
Step 3: 用 OLS 回归 Y_tilde ~ D_tilde，系数即为处理效应 θ

关键技术点：
– 使用交叉拟合（Cross-Fitting）避免过拟合偏差
– 使用梯度提升树作为一阶段的 ML 学习器
"""
X_cols = [c for c in df.columns if c.startswith("X")]
X = df[X_cols].values
D = df["D"].values
Y = df["Y"].values
n = len(Y)

print("\\n" + "=" * 60)
print("手动实现 DML（Partially Linear Model）")
print("=" * 60)

# ———————————————————-
# Step 1 & 2: 交叉拟合（Cross-Fitting）获取残差
# ———————————————————-
kf = KFold(n_splits=n_splits, shuffle=True, random_state=seed)

# ML 学习器（可替换为任意 ML 模型）
ml_y = GradientBoostingRegressor(
n_estimators=200, max_depth=4, learning_rate=0.05,
subsample=0.8, random_state=seed
)
ml_d = GradientBoostingRegressor(
n_estimators=200, max_depth=4, learning_rate=0.05,
subsample=0.8, random_state=seed
)

# 交叉拟合预测
Y_hat = cross_val_predict(ml_y, X, Y, cv=kf)
D_hat = cross_val_predict(ml_d, X, D, cv=kf)

# 计算残差（正交化/去偏）
Y_tilde = Y – Y_hat # 结果残差
D_tilde = D – D_hat # 处理残差

print(f" Y 残差均值: {Y_tilde.mean():.6f}")
print(f" D 残差均值: {D_tilde.mean():.6f}")

# ———————————————————-
# Step 3: 二阶段 OLS 回归残差
# ———————————————————-
theta_hat = np.sum(D_tilde * Y_tilde) / np.sum(D_tilde ** 2)

# 标准误估计（异方差稳健）
residuals = Y_tilde – theta_hat * D_tilde
V = np.mean(D_tilde ** 2 * residuals ** 2) / (np.mean(D_tilde ** 2) ** 2)
se = np.sqrt(V / n)

# 置信区间
ci_lower = theta_hat – 1.96 * se
ci_upper = theta_hat + 1.96 * se

# t 统计量和 p 值
t_stat = theta_hat / se
p_value = 2 * (1 – stats.norm.cdf(np.abs(t_stat)))

print(f"\\n ── DML 估计结果 ──")
print(f" θ_hat = {theta_hat:.4f}")
print(f" 标准误 = {se:.4f}")
print(f" t 统计量 = {t_stat:.4f}")
print(f" p 值 = {p_value:.6f}")
print(f" 95% CI = [{ci_lower:.4f}, {ci_upper:.4f}]")

return theta_hat, se, ci_lower, ci_upper

# ============================================================
# 第四部分：使用 EconML 库实现 DML
# ============================================================
def dml_econml(df, seed=42):
"""使用微软 EconML 库的 LinearDML 实现"""
try:
from econml.dml import LinearDML
except ImportError:
print("\\n[提示] 请安装 econml: pip install econml")
return None, None, None, None

X_cols = [c for c in df.columns if c.startswith("X")]
X = df[X_cols].values
D = df["D"].values.reshape(-1, 1)
Y = df["Y"].values

print("\\n" + "=" * 60)
print("EconML 库实现 DML（LinearDML）")
print("=" * 60)

# 配置 DML 模型
dml = LinearDML(
model_y=GradientBoostingRegressor(
n_estimators=200, max_depth=4,
learning_rate=0.05, subsample=0.8, random_state=seed
),
model_t=GradientBoostingRegressor(
n_estimators=200, max_depth=4,
learning_rate=0.05, subsample=0.8, random_state=seed
),
cv=5,
random_state=seed,
)

# 拟合模型
dml.fit(Y, D, X=None, W=X)

# 获取结果
theta_hat = dml.const_marginal_effect().flatten()[0]
inference = dml.const_marginal_effect_inference()
ci = inference.conf_int(alpha=0.05)
se = inference.stderr.flatten()[0]
p_value = inference.pvalue().flatten()[0]

ci_lower = ci[0].flatten()[0]
ci_upper = ci[1].flatten()[0]

print(f"\\n ── EconML 估计结果 ──")
print(f" θ_hat = {theta_hat:.4f}")
print(f" 标准误 = {se:.4f}")
print(f" p 值 = {p_value:.6f}")
print(f" 95% CI = [{ci_lower:.4f}, {ci_upper:.4f}]")

return theta_hat, se, ci_lower, ci_upper

# ============================================================
# 第五部分：不同 ML 学习器对比实验
# ============================================================
def compare_ml_learners(df, seed=42):
"""对比不同一阶段 ML 学习器对 DML 估计的影响"""
X_cols = [c for c in df.columns if c.startswith("X")]
X = df[X_cols].values
D = df["D"].values
Y = df["Y"].values
n = len(Y)

kf = KFold(n_splits=5, shuffle=True, random_state=seed)

learners = {
"Lasso (线性)": (LassoCV(cv=5, random_state=seed), LassoCV(cv=5, random_state=seed)),
"随机森林": (
RandomForestRegressor(n_estimators=200, max_depth=6, random_state=seed, n_jobs=-1),
RandomForestRegressor(n_estimators=200, max_depth=6, random_state=seed, n_jobs=-1),
),
"梯度提升树 (GBRT)": (
GradientBoostingRegressor(n_estimators=200, max_depth=4, learning_rate=0.05, random_state=seed),
GradientBoostingRegressor(n_estimators=200, max_depth=4, learning_rate=0.05, random_state=seed),
),
}

print("\\n" + "=" * 60)
print("不同 ML 学习器对比实验")
print("=" * 60)

results = []
for name, (ml_y, ml_d) in learners.items():
Y_hat = cross_val_predict(ml_y, X, Y, cv=kf)
D_hat = cross_val_predict(ml_d, X, D, cv=kf)

Y_tilde = Y – Y_hat
D_tilde = D – D_hat

theta = np.sum(D_tilde * Y_tilde) / np.sum(D_tilde ** 2)
resid = Y_tilde – theta * D_tilde
V = np.mean(D_tilde ** 2 * resid ** 2) / (np.mean(D_tilde ** 2) ** 2)
se = np.sqrt(V / n)

results.append({
"学习器": name,
"θ_hat": theta,
"标准误": se,
"CI_lower": theta – 1.96 * se,
"CI_upper": theta + 1.96 * se,
})

print(f" {name:20s} | θ = {theta:.4f} ± {se:.4f} "
f"95% CI = [{theta – 1.96 * se:.4f}, {theta + 1.96 * se:.4f}]")

return pd.DataFrame(results)

# ============================================================
# 第六部分：可视化
# ============================================================
def plot_results(true_theta, theta_ols, theta_dml, se_dml,
theta_econml=None, se_econml=None,
learner_df=None):
"""可视化估计结果对比"""

fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# — 图1：方法对比 —
ax1 = axes[0]
methods = ["OLS (有偏)", "DML (手动)"]
estimates = [theta_ols, theta_dml]
ses = [0, se_dml]
colors = ["#e74c3c", "#2ecc71"]

if theta_econml is not None:
methods.append("DML (EconML)")
estimates.append(theta_econml)
ses.append(se_econml)
colors.append("#3498db")

y_pos = np.arange(len(methods))
ax1.barh(y_pos, estimates, height=0.4, color=colors, alpha=0.8, edgecolor="gray")

for i, (est, se) in enumerate(zip(estimates, ses)):
if se > 0:
ax1.errorbar(est, i, xerr=1.96 * se, fmt="none", color="black",
capsize=5, linewidth=1.5)

ax1.axvline(x=true_theta, color="black", linestyle="–", linewidth=2, label=f"真实值 θ={true_theta}")
ax1.set_yticks(y_pos)
ax1.set_yticklabels(methods, fontsize=12)
ax1.set_xlabel("处理效应估计值", fontsize=12)
ax1.set_title("方法对比：OLS vs DML", fontsize=14, fontweight="bold")
ax1.legend(fontsize=11, loc="lower right")
ax1.grid(axis="x", alpha=0.3)

# — 图2：不同学习器对比 —
ax2 = axes[1]
if learner_df is not None:
y_pos2 = np.arange(len(learner_df))
ax2.barh(y_pos2, learner_df["θ_hat"], height=0.4,
color=["#9b59b6", "#1abc9c", "#f39c12"], alpha=0.8, edgecolor="gray")

for i, row in learner_df.iterrows():
ax2.errorbar(row["θ_hat"], i, xerr=1.96 * row["标准误"],
fmt="none", color="black", capsize=5, linewidth=1.5)

ax2.axvline(x=true_theta, color="black", linestyle="–", linewidth=2, label=f"真实值 θ={true_theta}")
ax2.set_yticks(y_pos2)
ax2.set_yticklabels(learner_df["学习器"], fontsize=12)
ax2.set_xlabel("处理效应估计值", fontsize=12)
ax2.set_title("不同 ML 学习器对比", fontsize=14, fontweight="bold")
ax2.legend(fontsize=11, loc="lower right")
ax2.grid(axis="x", alpha=0.3)

plt.tight_layout()
plt.savefig("dml_results.png", dpi=150, bbox_inches="tight")
plt.close()
print("\\n[图表已保存] dml_results.png")

# ============================================================
# 第七部分：蒙特卡洛模拟（验证无偏性）
# ============================================================
def monte_carlo_simulation(n_sim=200, n=1000, p=10):
"""蒙特卡洛模拟验证 DML 的无偏性"""
print("\\n" + "=" * 60)
print(f"蒙特卡洛模拟（{n_sim} 次重复）")
print("=" * 60)

TRUE_THETA = 2.0
ols_estimates = []
dml_estimates = []

for sim in range(n_sim):
np.random.seed(sim)
X = np.random.normal(0, 1, size=(n, p))
g_X = 2 * np.sin(X[:, 0] * X[:, 1]) + X[:, 2] ** 2 – 1.5 * np.abs(X[:, 3])
m_X = 0.8 * X[:, 0] + 0.5 * X[:, 1] ** 2 – 0.6 * X[:, 2]

epsilon = np.random.normal(0, 1, n)
eta = np.random.normal(0, 1, n)
D = m_X + eta
Y = TRUE_THETA * D + g_X + epsilon

# OLS
from sklearn.linear_model import LinearRegression
X_ols = np.column_stack([X, D])
ols_theta = LinearRegression().fit(X_ols, Y).coef_[-1]
ols_estimates.append(ols_theta)

# DML
kf = KFold(n_splits=3, shuffle=True, random_state=sim)
ml_y = GradientBoostingRegressor(n_estimators=100, max_depth=3, random_state=sim)
ml_d = GradientBoostingRegressor(n_estimators=100, max_depth=3, random_state=sim)
Y_hat = cross_val_predict(ml_y, X, Y, cv=kf)
D_hat = cross_val_predict(ml_d, X, D, cv=kf)
Y_tilde = Y – Y_hat
D_tilde = D – D_hat
dml_theta = np.sum(D_tilde * Y_tilde) / np.sum(D_tilde ** 2)
dml_estimates.append(dml_theta)

if (sim + 1) % 50 == 0:
print(f" 已完成 {sim + 1}/{n_sim} 次模拟…")

ols_arr = np.array(ols_estimates)
dml_arr = np.array(dml_estimates)

print(f"\\n ── 蒙特卡洛结果 ──")
print(f" 真实值: θ = {TRUE_THETA}")
print(f" OLS 均值: {ols_arr.mean():.4f} 偏差: {ols_arr.mean() – TRUE_THETA:.4f} 标准差: {ols_arr.std():.4f}")
print(f" DML 均值: {dml_arr.mean():.4f} 偏差: {dml_arr.mean() – TRUE_THETA:.4f} 标准差: {dml_arr.std():.4f}")

# 画分布图
fig, ax = plt.subplots(figsize=(10, 5))
ax.hist(ols_arr, bins=30, alpha=0.5, color="#e74c3c", label=f"OLS (均值={ols_arr.mean():.3f})", density=True)
ax.hist(dml_arr, bins=30, alpha=0.5, color="#2ecc71", label=f"DML (均值={dml_arr.mean():.3f})", density=True)
ax.axvline(x=TRUE_THETA, color="black", linestyle="–", linewidth=2, label=f"真实值 θ={TRUE_THETA}")
ax.set_xlabel("处理效应估计值", fontsize=12)
ax.set_ylabel("密度", fontsize=12)
ax.set_title(f"蒙特卡洛模拟：OLS vs DML（{n_sim} 次）", fontsize=14, fontweight="bold")
ax.legend(fontsize=12)
ax.grid(alpha=0.3)
plt.tight_layout()
plt.savefig("dml_monte_carlo.png", dpi=150, bbox_inches="tight")
plt.close()
print("[图表已保存] dml_monte_carlo.png")

return ols_arr, dml_arr

# ============================================================
# 主程序
# ============================================================
if __name__ == "__main__":

# 1. 生成模拟数据
df, true_theta = generate_dml_data(n=2000, p=10, seed=42)

# 2. 朴素 OLS 估计
theta_ols = naive_ols(df)
print(f"\\n 朴素 OLS 估计: θ_OLS = {theta_ols:.4f} (真实值: {true_theta})")

# 3. 手动 DML 估计
theta_dml, se_dml, ci_l, ci_u = dml_manual(df)

# 4. EconML 库 DML 估计
theta_econml, se_econml, _, _ = dml_econml(df)

# 5. 不同 ML 学习器对比
learner_df = compare_ml_learners(df)

# 6. 可视化
plot_results(true_theta, theta_ols, theta_dml, se_dml,
theta_econml, se_econml, learner_df)

# 7. 蒙特卡洛模拟
ols_mc, dml_mc = monte_carlo_simulation(n_sim=200, n=1000, p=10)

print("\\n" + "=" * 60)
print("全部分析完成！")
print("=" * 60)

作者寄语：DML 是将机器学习引入因果推断的一座桥梁。它不是万能的（仍需满足条件独立假设/无混淆假设），但在控制变量维度高、非线性关系复杂的场景下，相比传统 OLS 有本质的提升。建议结合自己的研究问题理解其适用条件，而非盲目套用。

如果觉得有帮助，欢迎点赞收藏关注，后续会更新更多因果推断方法（CATE、DID+ML、因果森林等）的实战教程！