附录 I — 双重机器学习 Stata 实操

本附录配合第七章《控制变量既多又非线性时：从 FWL 到双重机器学习》使用。正文讲清楚问题、适用场景、前提和解决思路，本附录讲怎么落地：把 FWL 的「两边同时剔除」手算验证一遍，用模拟对比天真估计与正交估计的偏误，再在官方示例数据上跑通 ddml 的部分线性、交互、IV 三种模型，最后给面板固定效应一个最小示意与局限提醒。

运行环境：Jupyter Notebook (.ipynb)，可在 VS Code 中用 nbstata 内核执行，或通过 Stata-MCP 调用本地 Stata。习惯直接在 Stata 里操作的读者，用配套的 A07_fwl_ddml.do (纯代码版)。

数据：H.2、H.3 用内置模拟数据 (set seed 保证可复现)；H.4 用 401(k) 数据 sipp1991.dta，H.5 用 cattaneo2，H.6 用 AJR.dta，H.7 用 nlswork——全部在线读取，不占本地仓库。

需要的外部命令：ddml、pystacked、lassopack (提供 rlasso)、rforest，H.1 会在线安装 (已装则跳过)。其中 pystacked 依赖 Stata 的 Python 集成 (需已装 Python + scikit-learn)；若未配置，相关行用 capture 跳过，主线改用 reg / rlasso / rforest 等纯 Stata 学习器仍可运行。

稳健性设计：核心命令用官方规范语法；环境敏感的命令 (尤其 pystacked)用 capture noisily 包裹，失败不影响主线；version 17 固定。

关于输出：本附录的代码单元格尚未在本仓库执行，输出为空，等待在装有 Stata 的环境中运行填充。每个模块都写明了「看哪一行、如何解读」；模拟模块的真值由数据生成过程 (DGP)设定。

一句话边界 (与正文一致)：DDML 治的是「控制变量既多又非线性」带来的模型误设，前提是选择性可观测；它不解决不可观测的遗漏变量问题。

I.1 环境准备

固定 Stata 版本，在线安装本附录用到的命令，已装的会跳过。pystacked 若因缺少 Python 集成而报错，不影响后面用纯 Stata 学习器跑主线。

version 17
set more off

* DDML 及相关命令(在线安装，已装则跳过)
cap which ddml
if _rc ssc install ddml, replace
cap which pystacked
if _rc ssc install pystacked, replace
cap which rlasso                      // 由 lassopack 提供 rlasso / cvlasso / lasso2
if _rc ssc install lassopack, replace
cap which rforest
if _rc ssc install rforest, replace

* 注：pystacked 需要 Stata 的 Python 集成(已装 Python + scikit-learn)。
*     若本机未配置，本附录中 pystacked 的行会被 capture 跳过，
*     主线改用 reg / rlasso / rforest 等纯 Stata 学习器仍可运行。

Running D:\stata19/profile.do ...

I.2 FWL 手算验证：三种回归给出同一个系数

对应正文「FWL：partial out 是基本功」一节，以及韦恩图 图 7.1 和流程图 图 7.2。

我们造一份 x1、x2 相关的模拟数据，真值 x1 的系数是 2。分三种做法估它：(1) x1、x2 一起进的完整回归；(2) 把 x2 从 y 和 x1 两边都 partial out，再残差对残差；(3) 只从 y 一边剔除 x2、x1 用原始值的错误做法。看前两者是否给出同一个系数，第三种为什么偏。

clear
set seed 20260710
set obs 500
gen x2 = rnormal()
gen x1 = 0.7*x2 + rnormal()             // X1 与 X2 相关，制造需要 partial out 的场景
gen y  = 1 + 2*x1 + 3*x2 + rnormal()    // 真值：beta1 = 2

* (1) 完整回归：X1 与 X2 一起进模型
reg y x1 x2
scalar b_full = _b[x1]

* (2) FWL：把 X2 从 Y 和 X1 两边同时 partial out，再残差对残差
reg y x2
predict yt,  resid                      // Y 对 X2 的残差
reg x1 x2
predict x1t, resid                      // X1 对 X2 的残差
reg yt x1t                              // 残差对残差，系数应等于 b_full
scalar b_fwl = _b[x1t]

* (3) 只对 Y 一边剔除 X2(错误做法：漏了 X1 那一边)
reg yt x1                               // 用原始 x1，不是残差 x1t
scalar b_wrong = _b[x1]

di as txt "b_full = " as res %6.4f b_full ///
   as txt "   b_fwl = " as res %6.4f b_fwl ///
   as txt "   b_wrong = " as res %6.4f b_wrong




Number of observations (_N) was 0, now 500.





      Source |       SS           df       MS      Number of obs   =       500
-------------+----------------------------------   F(2, 497)       =   5897.04
       Model |  12228.8592         2   6114.4296   Prob > F        =    0.0000
    Residual |  515.321531       497  1.03686425   R-squared       =    0.9596
-------------+----------------------------------   Adj R-squared   =    0.9594
       Total |  12744.1807       499  25.5394404   Root MSE        =    1.0183

------------------------------------------------------------------------------
           y | Coefficient  Std. err.      t    P>|t|     [95% conf. interval]
-------------+----------------------------------------------------------------
          x1 |      2.070      0.043    48.02   0.000        1.986       2.155
          x2 |      2.850      0.054    52.71   0.000        2.744       2.956
       _cons |      0.985      0.046    21.62   0.000        0.895       1.074
------------------------------------------------------------------------------



      Source |       SS           df       MS      Number of obs   =       500
-------------+----------------------------------   F(1, 498)       =   1685.95
       Model |  9838.16692         1  9838.16692   Prob > F        =    0.0000
    Residual |  2906.01382       498  5.83536912   R-squared       =    0.7720
-------------+----------------------------------   Adj R-squared   =    0.7715
       Total |  12744.1807       499  25.5394404   Root MSE        =    2.4157

------------------------------------------------------------------------------
           y | Coefficient  Std. err.      t    P>|t|     [95% conf. interval]
-------------+----------------------------------------------------------------
          x2 |      4.329      0.105    41.06   0.000        4.121       4.536
       _cons |      0.956      0.108     8.85   0.000        0.744       1.169
------------------------------------------------------------------------------



      Source |       SS           df       MS      Number of obs   =       500
-------------+----------------------------------   F(1, 498)       =    239.18
       Model |  267.884165         1  267.884165   Prob > F        =    0.0000
    Residual |  557.765137       498  1.12001031   R-squared       =    0.3245
-------------+----------------------------------   Adj R-squared   =    0.3231
       Total |  825.649301       499  1.65460782   Root MSE        =    1.0583

------------------------------------------------------------------------------
          x1 | Coefficient  Std. err.      t    P>|t|     [95% conf. interval]
-------------+----------------------------------------------------------------
          x2 |      0.714      0.046    15.47   0.000        0.624       0.805
       _cons |     -0.014      0.047    -0.29   0.772       -0.107       0.079
------------------------------------------------------------------------------



      Source |       SS           df       MS      Number of obs   =       500
-------------+----------------------------------   F(1, 498)       =   2310.33
       Model |  2390.69229         1  2390.69229   Prob > F        =    0.0000
    Residual |  515.321531       498  1.03478219   R-squared       =    0.8227
-------------+----------------------------------   Adj R-squared   =    0.8223
       Total |  2906.01382       499    5.823675   Root MSE        =    1.0172

------------------------------------------------------------------------------
          yt | Coefficient  Std. err.      t    P>|t|     [95% conf. interval]
-------------+----------------------------------------------------------------
         x1t |      2.070      0.043    48.07   0.000        1.986       2.155
       _cons |     -0.000      0.045    -0.00   1.000       -0.089       0.089
------------------------------------------------------------------------------



      Source |       SS           df       MS      Number of obs   =       500
-------------+----------------------------------   F(1, 498)       =    623.00
       Model |  1615.02567         1  1615.02567   Prob > F        =    0.0000
    Residual |  1290.98816       498   2.5923457   R-squared       =    0.5558
-------------+----------------------------------   Adj R-squared   =    0.5549
       Total |  2906.01382       499    5.823675   Root MSE        =    1.6101

------------------------------------------------------------------------------
          yt | Coefficient  Std. err.      t    P>|t|     [95% conf. interval]
-------------+----------------------------------------------------------------
          x1 |      1.399      0.056    24.96   0.000        1.289       1.509
       _cons |      0.005      0.072     0.07   0.943       -0.136       0.147
------------------------------------------------------------------------------


b_full = 2.0703   b_fwl = 2.0703   b_wrong = 1.3986

输出怎么读：

  • b_full 与 b_fwl 应当完全相等 (都约等于真值 2)，这正是 FWL 定理：完整回归里 x1 的系数，等于「两边都对 x2 取残差后，残差对残差」的斜率。
  • b_wrong 明显偏离 2。它只把 x2 从 y 一边剔除，x1 仍带着与 x2 的相关，于是 x2 的影响漏进了系数。这对应正文强调的那句：要得到干净的 \(\theta\)，必须把控制变量从处理变量和结果变量两边都剔除。
  • 教学连接：把这里的 x1 换成处理变量 D、x2 换成一堆控制变量 X，第 (3) 种「只剔一边」的错误，正是下一节天真机器学习失败的根源。

I.3 天真估计 vs 正交估计：偏误对比模拟

对应正文「补救：对处理变量也 partial out」一节，以及乘积偏误 式 7.8。

数据生成过程 (DGP)设为 \[ D = m_0(X) + V,\qquad Y = \theta_0 D + g_0(X) + U,\qquad \theta_0 = 1, \] 其中 \(m_0\)、\(g_0\) 是二次加交互的非线性函数，线性设定不足以吸收。我们用带正则化的灵活学习器 (rlasso 配二次字典)拟合两个辅助函数，然后比较两种估计：

  • 天真：只把 \(X\) 从 \(Y\) 一边扣掉 (残差 yr)，再用原始 d 回归；
  • 正交：对 d 也扣掉 \(X\) (残差 dr)，用残差对残差。

重复 200 次，看两者的抽样分布中心离真值 1 有多远。

clear all
set seed 20260710

cap program drop ddmlsim
program define ddmlsim, rclass
    drop _all
    quietly set obs 500
    forvalues j = 1/10 {
        quietly gen x`j' = rnormal()
    }
    * 非线性混淆
    quietly gen d = 0.8*x1 + 0.5*x2^2 - 0.4*x3*x4 + rnormal()
    quietly gen y = 1*d + 1.0*x1^2 + 0.6*x2 - 0.5*x3*x5 + rnormal()   // theta_0 = 1

    * 带正则化的灵活学习器(二次字典)，拟合两个辅助函数
    quietly rlasso y c.(x1-x10)##c.(x1-x10)
    quietly predict double yhat, xb
    quietly gen double yr = y - yhat                 // Y 的残差

    quietly rlasso d c.(x1-x10)##c.(x1-x10)
    quietly predict double dhat, xb
    quietly gen double dr = d - dhat                 // D 的残差

    * 天真：只对 Y 一边扣掉 X，再用原始 D 回归
    quietly reg yr d
    return scalar naive = _b[d]

    * 正交：对 D 也扣掉 X，用残差对残差
    quietly reg yr dr
    return scalar ortho = _b[dr]
end

simulate naive = r(naive) ortho = r(ortho), reps(200): ddmlsim
summ naive ortho
* 真值 theta_0 = 1：naive 的均值应明显偏离 1，ortho 的均值应接近 1







      Command: ddmlsim
        naive: r(naive)
        ortho: r(ortho)

Simulations (200): .........10.........20.........30.........40.........50.....
> ....60.........70.........80.........90.........100.........110.........120..
> .......130.........140.........150.........160.........170.........180.......
> ..190.........200 done


    Variable |        Obs        Mean    Std. dev.       Min        Max
-------------+---------------------------------------------------------
       naive |        200    .6406176    .0450327   .5283877   .7854204
       ortho |        200    1.029012    .0468403   .8902555   1.145817

输出怎么读：

  • summ 给出 naive 与 ortho 两列在 200 次重复中的均值。真值是 \(\theta_0 = 1\)。
  • naive 的均值应系统性偏离 1：辅助函数的正则化误差没有被抵消，乘上 \(m_0(X)\) 留在系数里 (正文 式 7.6)。
  • ortho 的均值应接近 1：对 d 也做了 partial out，偏误变成两个估计误差的乘积 (正文 式 7.8)，量级小得多。
  • 想更严格，还可以在此基础上加交叉拟合 (用留出折算残差)，进一步消掉过拟合带来的虚假关联——那正是 ddml 从 H.4 起自动做的事。

I.4 部分线性模型 ddml partial (401k 数据)

对应正文「DDML 的完整流程」一节与表 表 7.1。这是最常用的起点：处理效应大体同质、处理变量与控制变量可线性分离。

官方 401(k) 示例：net_tfa (净金融资产)为结果，e401 (是否符合 401k 计划资格)为处理，X 是一组家庭特征。四步流程——ddml init → 指定 E[Y|X]、E[D|X] 学习器 → ddml crossfit → ddml estimate——正是正文 图 7.7 画的骨架。

use https://github.com/aahrens1/ddml/raw/master/data/sipp1991.dta, clear

global Y net_tfa
global D e401
global X tw age inc fsize educ db marr twoearn pira hown
set seed 42

ddml init partial, kfolds(2)

* 结果侧 E[Y|X]：先加一个纯 Stata 线性学习器(必定可跑)，再加 pystacked 随机森林
ddml E[Y|X]: reg $Y $X
capture noisily ddml E[Y|X]: pystacked $Y $X, type(reg) method(rf)

* 处理侧 E[D|X]：同样两个学习器
ddml E[D|X]: reg $D $X
capture noisily ddml E[D|X]: pystacked $D $X, type(reg) method(rf)

ddml crossfit
ddml estimate, robust allcombos








Learner Y1_reg added successfully.

Learner Y2_pystacked added successfully.

Learner D1_reg added successfully.

Learner D2_pystacked added successfully.

Cross-fitting E[y|X] equation: net_tfa
Cross-fitting fold 1 2 ...completed cross-fitting
Cross-fitting E[D|X] equation: e401
Cross-fitting fold 1 2 ...completed cross-fitting



Model:                  partial, crossfit folds k=2, resamples r=1
Mata global (mname):    m0
Dependent variable (Y): net_tfa
 net_tfa learners:      Y1_reg Y2_pystacked
D equations (1):        e401
 e401 learners:         D1_reg D2_pystacked

DDML estimation results:
spec  r     Y learner     D learner         b        SE 
   1  1        Y1_reg        D1_reg  5397.208 (1130.776)
   2  1        Y1_reg  D2_pystacked  6715.427  (879.778)
*  3  1  Y2_pystacked        D1_reg  7099.006 (1153.524)
   4  1  Y2_pystacked  D2_pystacked  7097.568  (775.359)
* = minimum MSE specification for that resample.

Min MSE DDML model
y-E[y|X]  = y-Y2_pystacked_1                       Number of obs   =      9915
D-E[D|X]  = D-D1_reg_1 
------------------------------------------------------------------------------
             |               Robust
     net_tfa | Coefficient  std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
        e401 |   7099.006   1153.524     6.15   0.000     4838.140    9359.871
       _cons |   -237.128    357.443    -0.66   0.507     -937.703     463.448
------------------------------------------------------------------------------

命令逐项解读：

  • ddml init partial, kfolds(2)：声明部分线性模型，2 折交叉拟合 (对应流程第 1 步「选模型」)。
  • ddml E[Y|X]: reg … 与 ddml E[D|X]: reg …：分别为结果侧和处理侧指定学习器 (第 2 步)。这里每侧加了两个学习器——一个纯 Stata 的 reg (必定可跑)，一个 pystacked 随机森林 (灵活，但需 Python，故用 capture noisily 包裹)。
  • ddml crossfit：在各折上估计辅助函数、算出 out-of-fold 残差 (第 3 步)。
  • ddml estimate, robust allcombos：用正交矩条件汇总 (第 4 步)；allcombos 会把不同学习器组合的结果都列出来。

输出怎么读：

  • 看 E[Y|X] 与 E[D|X] 各学习器的交叉验证 MSE，判断哪个辅助模型拟合更好。
  • 主表里 e401 那一行的系数就是处理效应 \(\hat\theta\)，配 robust 标准误和置信区间。
  • 若 pystacked 被跳过 (没配 Python)，allcombos 只会显示基于 reg 的组合——主线仍然完整。

I.5 交互模型 ddml interactive (cattaneo2 数据)

对应表 表 7.1 的 interactive 行：处理为二元、效应完全异质、不假设与 \(X\) 可线性分离，目标参数是 ATE。

官方 cattaneo2 示例：bweight (出生体重)为结果，mbsmoke (孕期是否吸烟)为二元处理。注意交互模型对处理的两种取值分别建结果模型，故写作 E[Y|X,D]；处理侧 E[D|X] 是个分类问题，用 logit 或分类学习器。

webuse cattaneo2, clear

global Y bweight
global D mbsmoke
global X mage prenatal1 mmarried fbaby medu
set seed 42

ddml init interactive, kfolds(5) reps(5)

* 结果侧 E[Y|X,D](交互模型对处理分别建模)
ddml E[Y|X,D]: reg $Y $X
capture noisily ddml E[Y|X,D]: pystacked $Y $X, type(reg) method(gradboost)

* 处理侧 E[D|X]：二元处理用 logit / 分类学习器
ddml E[D|X]: logit $D $X
capture noisily ddml E[D|X]: pystacked $D $X, type(class) method(gradboost)

ddml crossfit
ddml estimate


(Excerpt from Cattaneo (2010) Journal of Econometrics 155: 138–154)





warning - model m0 already exists
all existing model results and variables will
be dropped and model m0 will be re-initialized

Learner Y1_reg added successfully.

Learner Y2_pystacked added successfully.

Learner D1_logit added successfully.

Learner D2_pystacked added successfully.

Cross-fitting E[y|X,D] equation: bweight
Resample 1...
Cross-fitting fold 1 2 3 4 5 ...completed cross-fitting
Resample 2...
Cross-fitting fold 1 2 3 4 5 ...completed cross-fitting
Resample 3...
Cross-fitting fold 1 2 3 4 5 ...completed cross-fitting
Resample 4...
Cross-fitting fold 1 2 3 4 5 ...completed cross-fitting
Resample 5...
Cross-fitting fold 1 2 3 4 5 ...completed cross-fitting
Cross-fitting E[D|X] equation: mbsmoke
Resample 1...
Cross-fitting fold 1 r(7102);

输出怎么读：

  • ddml init interactive, kfolds(5) reps(5)：交互模型，5 折、重复 5 次 (reps 多次分折取平均，降低随机性)。
  • 结果侧用 E[Y|X,D] (对处理组、控制组分别拟合)，处理侧 E[D|X] 用 logit (纯 Stata)加 pystacked 分类器 (type(class)，capture 包裹)。
  • 主表报告的是 ATE (平均处理效应)。
  • 教学连接：与 H.4 对照，能看出「选模型」这一步如何改变矩条件和目标参数——正文说的「第二步选模型决定后面一切」。

I.6 部分线性 IV 模型 ddml iv (AJR 数据)

对应表 表 7.1 的 iv 行：处理变量内生、且有有效工具变量时使用。

IV 需要一个有效工具，401k 和 cattaneo2 都没有，所以这里换用官方 AJR.dta (Acemoglu-Johnson-Robinson 的制度与经济发展数据)：logpgp95 (人均 GDP 对数)为结果，avexpr (征收保护指数，内生)为处理，logem4 (殖民者死亡率对数)为工具，X 是一组地理与制度控制。

这一节的学习器用 rforest (纯 Stata 插件随机森林，不需要 Python)，演示没有 pystacked 时怎样用灵活学习器。多出一条 E[Z|X]：工具变量也要对 \(X\) 做 partial out。

use https://statalasso.github.io/dta/AJR.dta, clear

global Y logpgp95
global D avexpr
global Z logem4
global X lat_abst edes1975 avelf temp* humid* steplow-oilres
set seed 42

ddml init iv, kfolds(30)

ddml E[Y|X]: reg $Y $X
capture noisily ddml E[Y|X], vtype(none): rforest $Y $X, type(reg)
ddml E[D|X]: reg $D $X
capture noisily ddml E[D|X], vtype(none): rforest $D $X, type(reg)
ddml E[Z|X]: reg $Z $X
capture noisily ddml E[Z|X], vtype(none): rforest $Z $X, type(reg)

qui ddml crossfit
ddml estimate, robust








warning - model m0 already exists
all existing model results and variables will
be dropped and model m0 will be re-initialized

Learner Y1_reg added successfully.

Learner Y2_rforest added successfully.

Learner D1_reg added successfully.

Learner D2_rforest added successfully.

Learner Z1_reg added successfully.

Learner Z2_rforest added successfully.




Model:                  iv, crossfit folds k=30, resamples r=1
Mata global (mname):    m0
Dependent variable (Y): logpgp95
 logpgp95 learners:     Y1_reg Y2_rforest
D equations (1):        avexpr
 avexpr learners:       D1_reg D2_rforest
Z equations (1):        logem4
 logem4 learners:       Z1_reg Z2_rforest

DDML estimation results:
spec  r     Y learner     D learner         b        SE      Z learner
 mse  1    Y2_rforest    D2_rforest     0.772    (0.207)    Z2_rforest
mse = minimum MSE specification for that resample.

Min MSE DDML model
y-E[y|X]  = y-Y2_rforest_1                         Number of obs   =        64
D-E[D|X]  = D-D2_rforest_1 
Z-E[Z|X]  = Z-Z2_rforest_1 
------------------------------------------------------------------------------
             |               Robust
    logpgp95 | Coefficient  std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
      avexpr |      0.772      0.207     3.73   0.000        0.367       1.178
       _cons |     -0.012      0.101    -0.12   0.906       -0.210       0.186
------------------------------------------------------------------------------

输出怎么读：

  • 比 H.4 多了 ddml E[Z|X]：把 \(X\) 从工具 logem4 里也清走，这样识别用的是工具中与 \(X\) 无关的外生变异。
  • vtype(none) 是 rforest 作为学习器时的技术选项 (预测值类型)，照抄官方示例即可。
  • kfolds(30) 折数偏大，是因为 AJR 样本很小 (约 60 多个国家)，多分折能稳住估计。
  • 主表里 avexpr 的系数就是 IV 版处理效应。若 rforest 未安装、被 capture 跳过，主线仍用 reg 学习器给出结果。
  • 边界提醒：IV 版 DDML 处理的是「控制变量既多又非线性」下的工具变量估计；工具本身是否满足排他性，仍是 DDML 之外的识别假设 (正文 表 7.2)。

I.7 面板固定效应：别直接把 within 变换套截面 DDML (最小示意 + 局限)

对应正文 callout「面板固定效应：别直接套截面 DDML」。标准 ddml 是截面框架。面向政策评估的读者习惯写「固定效应 + 一堆控制变量」，容易顺手把面板 within 变换后直接套 ddml——这里演示这个做法，同时讲清它为什么只能当示意、不能当推荐。

用 nlswork (女性劳动收入面板)：ln_wage 为结果，union (是否工会)为处理，做去个体均值的 within 变换后，套截面部分线性 ddml。

webuse nlswork, clear
xtset idcode year
keep if !missing(ln_wage, union, age, ttl_exp, tenure, hours)

* 手动 within(去个体均值)变换
foreach v in ln_wage union age ttl_exp tenure hours {
    bysort idcode: egen double m_`v' = mean(`v')
    gen double w_`v' = `v' - m_`v'
}

global Y w_ln_wage
global D w_union
global X w_age w_ttl_exp w_tenure w_hours
set seed 42

ddml init partial, kfolds(2)
ddml E[Y|X]: reg $Y $X
capture noisily ddml E[Y|X]: pystacked $Y $X, type(reg) method(rf)
ddml E[D|X]: reg $D $X
capture noisily ddml E[D|X]: pystacked $D $X, type(reg) method(rf)
ddml crossfit
ddml estimate, robust


(National Longitudinal Survey of Young Women, 14-24 years old in 1968)


Panel variable: idcode (unbalanced)
 Time variable: year, 68 to 88, but with gaps
         Delta: 1 unit

(9,558 observations deleted)






warning - model m0 already exists
all existing model results and variables will
be dropped and model m0 will be re-initialized

Learner Y1_reg added successfully.

Learner Y2_pystacked added successfully.

Learner D1_reg added successfully.

Learner D2_pystacked added successfully.

Cross-fitting E[y|X] equation: w_ln_wage
Cross-fitting fold 1 2 ...completed cross-fitting
Cross-fitting E[D|X] equation: w_union
Cross-fitting fold 1 2 ...completed cross-fitting



Model:                  partial, crossfit folds k=2, resamples r=1
Mata global (mname):    m0
Dependent variable (Y): w_ln_wage
 w_ln_wage learners:    Y1_reg Y2_pystacked
D equations (1):        w_union
 w_union learners:      D1_reg D2_pystacked

DDML estimation results:
spec  r     Y learner     D learner         b        SE 
 mse  1        Y1_reg        D1_reg     0.102    (0.007)
mse = minimum MSE specification for that resample.

Min MSE DDML model
y-E[y|X]  = y-Y1_reg_1                             Number of obs   =     18976
D-E[D|X]  = D-D1_reg_1 
------------------------------------------------------------------------------
             |               Robust
   w_ln_wage | Coefficient  std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
     w_union |      0.102      0.007    14.05   0.000        0.088       0.116
       _cons |      0.000      0.002     0.02   0.983       -0.003       0.003
------------------------------------------------------------------------------

输出怎么读与局限：

  • w_union 的系数是 within 变换后、截面 ddml 给出的处理效应估计。它能跑出来，但不建议照此下结论。
  • 两个问题：(1) within 变换在个体内引入依赖，随机分折不再合理，更稳妥是按个体分折；(2) 未观测个体异质性用相关随机效应 (Mundlak/CRE)处理往往更稳。
  • 面板 DDML 仍在发展中，本模块只作风险提示，不推荐某一种固定方案。系统的面板实现见 R 包 xtdml (Clarke and Polselli 2026)；面板下的分折与 CRE 讨论见 Fuhr and Papies (2024)。

I.8 对口文献 (含复现资料)

  • FWL 与部分回归：Filoso (2013), Regression Anatomy…, Stata Journal 13(1)；中文推导与 Stata 实例见李金桐 (2023, 连享会 No.1221)。
  • DDML 奠基：Chernozhukov et al. (2018), Double/debiased machine learning…, The Econometrics Journal 21(1)。
  • ddml 官方导论与 Stata 实现：Ahrens et al. (2025), arXiv:2504.08324；Ahrens et al. (2024), ddml: Double/debiased machine learning in Stata, Stata Journal 24(1)。官方文档站：https://statalasso.github.io/docs/ddml/。
  • stacking 学习器：Ahrens et al. (2024), Model Averaging and Double Machine Learning，配套 pystacked。
  • 面板 DDML：Clarke and Polselli (2026), The Econometrics Journal 29(1)，R 包 xtdml：https://github.com/POLSEAN/xtdml；Fuhr and Papies (2024), arXiv:2409.01266。
  • 遗漏变量敏感性：Chernozhukov et al. (2022), Long Story Short: Omitted Variable Bias in Causal Machine Learning。

完整参考文献与延伸阅读见第七章正文「参考文献与延伸阅读」一节。本附录只列与各模块代码直接对口的条目。