经济学数据分析的方法

  • 文本分析：情感分析、主题模型（LDA）
  • 显著性检验： t 检验、F 检验

因果推断与政策评估（现代经济学核心）

  1. 双重差分法 DID 用于政策冲击前后对比
  2. 工具变量法 IV 解决内生性问题
  3. 断点回归 RDD 利用制度阈值做准自然实验
  4. 匹配方法 - 倾向得分匹配 PSM - 精确匹配、熵平衡匹配
  5. 合成控制法 SCM 构造 “反事实” 地区评估政策效果

固定效应模型（FE）： 就是假设每个地区有自己不变的特性（比如地理位置、文化等），模型会随机效应模型（RE）： 意思是假设地区间的差异是随机的，不影响核心变量的估计。

完整的实证流程： 1. 明确变量 → 2. 导入数据 → 3. 描述统计 → 4. 相关性分析 → 5. 基准回归 → 6. 稳健性检验 →7.异质性检验→8.结果分析

基准回归

  • 模型最简、设定最标准， 只放核心解释变量、被解释变量 + 最基础的控制变量，不做复杂调整、不加额外拓展变量。
  • 作为对照基准， 后面所有的稳健性检验、异质性分析、机制检验、拓展回归，全都要和基准回归结果对比。

稳健性检验，可以用因果推断
简单说就是换多种方法、换数据、换设定，反复验证你的实证结论是不是 “经得起折腾、不是偶然算出来的”。

常用的稳健性检验方法

  • 替换被解释变量 / 解释变量衡量指标 - 换一种算法、换代理变量，看符号、显著性是否一致。
  • 更换计量模型 - 比如：基准用 OLS，稳健性用固定效应、GMM、Tobit 等。
  • 缩尾 / 截尾处理 - 剔除极端异常值，避免离群值干扰结果。
  • 调整样本区间 - 删掉特殊年份（如疫情、金融危机）、删掉特殊行业 / 企业。
  • 改变控制变量组合 - 多加 / 少加控制变量，看核心变量结果稳不稳定。
  • 内生性处理（也算广义稳健） - 工具变量、滞后一期、PSM 等。
  • 分子样本回归 - 分国企 / 民企、大小企业分组，看结论是否普遍成立。

机制分析。中介变量分析。比如用来检验「政策对收入的影响」是否通过「中介变量」（比如教育水平）来传递。

异质性检验，为什么不同数据、不同组别、不同研究结果不一样，差异来源是什么、差异大不大。
比如研究「政策对收入的影响」 👉 不是所有人影响都一样： - 对年轻人影响大，对老年人影响小 - 对一线城市有效，对小城市没用 这时做异质性分析 = 分样本回归。

注意指标反映的内容，而不仅仅是指标本身。

描述性统计

按照时间和个体分组。

按照其他条件分组。

可以画图，进行可视化。柱状图，折线图。

个数，均值，变化趋势。分组。比例，比例变化。

回归分析

线性回归分析

用来发现相关性。是否显著。

线性回归 $Y=β0​+β1​X1​+β2​X2​+…+ϵ$。 用途 房价预测、销售额影响分析。
拟合模型：最小二乘法（OLS）求系数。

模型检验

  • R²：解释力（0~1，越接近 1 越好）。
  • P 值：系数是否显著（<0.05 通常有效）。
  • 残差：是否随机、正态、同方差。
  • F 检验：方差分析，某个变量是否显著。

与相关分析的关系： - 相关分析只衡量相关系数 r - 回归分析建立预测模型，给出方向与大小（系数）。

statsmodels supports specifying models using R-style formulas and pandas DataFrames. Here is a simple example using ordinary least squares:

import numpy as np
import statsmodels.api as sm
import statsmodels.formula.api as smf
# Load data
dat = sm.datasets.get_rdataset("Guerry", "HistData").data
# Fit regression model (using the natural log of one of the regressors)
results = smf.ols('Lottery ~ Literacy + np.log(Pop1831)', data=dat).fit()
# Inspect the results
print(results.summary())

输出结果分析

                            OLS Regression Results                            
==============================================================================
Dep. Variable:                Lottery   R-squared:                       0.348
Model:                            OLS   Adj. R-squared:                  0.333
Method:                 Least Squares   F-statistic:                     22.20
Date:                Sun, 17 May 2026   Prob (F-statistic):           1.90e-08
Time:                        20:57:57   Log-Likelihood:                -379.82
No. Observations:                  86   AIC:                             765.6
Df Residuals:                      83   BIC:                             773.0
Df Model:                           2                                         
Covariance Type:            nonrobust                                         
===================================================================================
                      coef    std err          t      P>|t|      [0.025      0.975]
-----------------------------------------------------------------------------------
Intercept         246.4341     35.233      6.995      0.000     176.358     316.510
Literacy           -0.4889      0.128     -3.832      0.000      -0.743      -0.235
np.log(Pop1831)   -31.3114      5.977     -5.239      0.000     -43.199     -19.424
==============================================================================
Omnibus:                        3.713   Durbin-Watson:                   2.019
Prob(Omnibus):                  0.156   Jarque-Bera (JB):                3.394
Skew:                          -0.487   Prob(JB):                        0.183
Kurtosis:                       3.003   Cond. No.                         702.
==============================================================================

Both Literacy and np.log(Pop1831) are statistically significant predictors of Lottery (p < 0.001). The model explains ~34.8% of the variance (R-squared = 0.348).

方差分析

线性回归的方差分析

import statsmodels.api as sm
from statsmodels.formula.api import ols

moore = sm.datasets.get_rdataset("Moore", "carData", cache=True) # load
data = moore.data
data = data.rename(columns={"partner.status": "partner_status"}) # make name pythonic

moore_lm = ols('conformity ~ C(fcategory, Sum)*C(partner_status, Sum)', data=data).fit()
table = sm.stats.anova_lm(moore_lm, typ=2) # Type 2 Anova DataFrame
print(table)

输出

                                              sum_sq    df          F    PR(>F)
C(fcategory, Sum)                          11.614700   2.0   0.276958  0.759564
C(partner_status, Sum)                    212.213778   1.0  10.120692  0.002874
C(fcategory, Sum):C(partner_status, Sum)  175.488928   2.0   4.184623  0.022572
Residual                                  817.763961  39.0        NaN       NaN

结果解读，看每项的 F 值是否显著（足够小）。

  • fcategory (p=0.760): 不显著，authoritarian 和 independent 两类女性的 conformity 得分差异没有统计意义。
  • partner_status (p=0.003): 显著，伴侣地位（低/高）对 conformity 有显著主效应。
  • 交互项 fcategory:partner_status (p=0.023): 显著，说明 fcategory 的效应随 partner_status 不同而变化——这也正是 Moore 数据集想要展示的经典交互效应。

方差分析 ANOVA 全程只用 F 检验，没有 t 检验，回归分析整体用 F，单个自变量用 t。因为两者此时是等价的。

替代/互补效应

做回归分析，看系数的正负。

怎么用回归看替代 / 互补效应，核心是引入交互项、判断交互项符号与显著性。
最常用、最标准的方法：
在回归里同时放入两个核心自变量 X₁ 和 X₂，看它们的交互项 X₁×X₂
基本回归模型 $Y=β0​+β1​X1​+β2​X2​+β3​(X1​×X2​)+ϵ$
看交互项系数 β₃ 就行：

  • β₃ ＜ 0，显著 → 替代效应 一个变强，另一个的作用被削弱，互相抵消。
  • β₃ ＞ 0，显著 → 互补效应 两个一起变强，效果互相放大，1+1＞2。

交互项

相乘

因果推断

在回归分析的基础上扩展。相关性不代表因果效应。

直接用 OLS 会有偏误。（偏差体现在哪里？）

IV（工具变量法）、双重差分法（DID）的方法和案例。也有其他因果推断的方法。

面板数据
截面数据

内生性问题

计量经济学里内生性是因果识别的最大敌人
因果推断方法（DID、RDD、IV、匹配、合成控制）全部都是为了克服内生性，从而得到因果效应

自变量与误差项相关。互相因果关系。

模型中的解释变量（自变量）与误差项存在相关性（Cov (X, u) ≠ 0），导致普通最小二乘法（OLS）估计量有偏、不一致。解决思路是寻找外生变异（工具变量、自然实验、固定效应、匹配等），切断 X 与 u 的关联，还原真实因果效应。

主要来源

  1. 遗漏变量（Omitted Variable Bias） 模型未纳入同时影响 X 与 Y 的关键变量，该变量被归入误差项 u，导致 X 与 u 相关。 例：研究 “教育年限→收入” 时，遗漏 “个人能力”—— 能力既影响教育选择，又影响收入，使教育系数被高估。
  2. 互为因果/反向因果（Simultaneity / Reverse Causality） X 与 Y 互为因果、双向决定，而非单向因果。 例：研究 “政府监管→企业污染” 时，污染越重的企业越易被监管，监管又会抑制污染，双向关系导致内生性。
  3. 测量误差（Measurement Error） 解释变量 X 存在观测误差，误差被纳入 u，使 X 的观测值与 u 相关。 例：用 “自我报告收入” 衡量真实收入，测量误差导致收入变量与误差项相关。
  4. 样本选择偏差（Sample Selection Bias） 样本并非随机抽取，而是由个体 “自选择” 或外部筛选导致，使 X 与 u 相关。 例：研究 “培训→工资” 时，只有预期收益高的人才会主动参加培训，样本不具代表性。

解决方法（按适用性排序）

  1. 工具变量法（IV, Instrumental Variable） 核心：找到一个外生变量 Z，满足： Z 与 X 高度相关（相关性）； Z 与误差项 u 不相关（外生性，即 Z 仅通过 X 影响 Y）。 适用：反向因果、遗漏变量；最经典、最常用。 例：用 “出生季度” 作为 “教育年限” 的工具（Angrist 教育回报研究）。
  2. 固定效应模型（FE, Fixed Effects） 核心：利用面板数据，控制个体 / 时间固定效应，吸收不随时间 / 个体变化的遗漏变量。 适用：遗漏不随时间变化的个体 / 地区特征（如能力、文化、制度）。 例：用企业面板数据控制企业固定效应，消除企业固有能力对 “投资→绩效” 的干扰。
  3. 双重差分法（DID, Difference-in-Differences） 核心：利用政策 / 冲击的 “准自然实验”，比较处理组与对照组在政策前后的差异，剔除共同趋势与个体固定效应。 适用：政策评估、外生冲击下的因果识别。 例：评估 “最低工资上调” 对就业的影响，对比政策实施前后，处理区与对照区的就业变化。
  4. 倾向得分匹配（PSM, Propensity Score Matching） 核心：为每个处理组个体匹配特征相似的对照组个体，平衡协变量分布，消除自选择偏差。 适用：X 为二元处理变量（如 “参加培训 / 不参加”）。
  5. 滞后变量（Lagged Variables） 核心：用 X 的滞后一期值作为当期 X 的工具或替代，利用 “过去不影响未来误差项” 的外生性。 适用：动态面板、弱反向因果。
  6. 赫克曼两阶段模型（Heckman Two-Step） 核心：先估计 “是否进入样本” 的选择方程，再将选择偏差纳入主回归，修正样本选择偏误。 适用：非随机样本、截断数据。

内生性检测（常用检验）

  • 豪斯曼检验（Hausman Test）：比较 OLS 与 IV/FE 估计，判断是否存在显著内生性。常用。
  • 杜宾 - 吴 - 豪斯曼检验（DWH Test）：检验解释变量是否内生。
  • 弱工具变量检验（F 检验）：工具变量第一阶段 F > 10 通常视为 “非弱工具”。

倾向得分匹配 PSM

  • 处理组（如政策干预、 treatment）和对照组在可观测特征上本来就不一样，直接比较会有偏差。 - PSM 为每个处理组个体匹配一个特征最相似的对照组个体，再比较结果差异。 - 适用场景截面数据/混合截面。

步骤

  • 估计倾向得分：用 Logit/Probit 模型拟合被处理的概率
  • 匹配：近邻匹配、半径匹配、核匹配、卡尺匹配
  • 平衡性检验：确保匹配后协变量在两组无显著差异
  • 计算平均处理效应（ATT）=匹配后处理组均值−匹配后对照组均值。
  • 对匹配后的数据线性回归或者 DID 回归
  1. 适用场景 截面数据 / 混合截面 关注平均处理效应 ATT 可观测混杂因素较多，但无不可观测个体固定效应 样本量较大，匹配效果好 4. 优点 直观、易解释 不依赖强线性假设 能处理非线性、非参数关系 5. 缺点 只控制可观测变量，无法解决遗漏变量偏误 对共同支撑域敏感 匹配方法选择会影响结果 无法处理随时间不变的个体固定效应

核心思想是通过「倾向得分」（即给定个体特征情况下，个体进入处理组的概率）将处理组和对照组中得分相近的个体进行匹配，使得匹配后两组除了是否接受处理外，其他可观测特征分布基本一致，从而消除可观测混淆变量的影响。

步骤

  1. 变量选择：明确处理变量（二元变量，如是否参与项目）、结果变量（关注的效应指标，如收入、绩效）和混淆变量（同时影响处理分配和结果的变量，如年龄、学历、行业等）。
  2. 估计倾向得分：通常使用Logit或Probit模型，以处理变量为被解释变量，混淆变量为解释变量，拟合得到每个个体的倾向得分。
  3. 匹配：选择匹配方法（常见的有最近邻匹配、卡尺匹配、核匹配、分层匹配等），对处理组和对照组个体进行匹配。
  4. 平衡性检验：验证匹配后两组的混淆变量是否无显著差异（标准化偏差一般要求小于10%），若不平衡则需调整模型或匹配方法。
  5. 估计因果效应：在匹配后的样本中，比较处理组和对照组结果变量的均值差异，即为ATT（或其他效应），同时通过自助法（Bootstrap）计算标准误进行显著性检验。
  6. 稳健性检验：可通过更换匹配方法、调整混淆变量、安慰剂检验等方式验证结果的可靠性。

具体例子

具体怎么做（最标准流程）

  1. 先做 PSM
    • 计算倾向得分
    • 最近邻匹配 / 半径匹配 / 核匹配
    • 检查平衡性：标准化偏差 < 10%，t 检验不显著
  2. 匹配后，可以用匹配后的样本跑回归。最常用两种：
    1. PSM + 多元回归（最常用）
      • Y ~ treat + X1 + X2 + X3 + …
      • 只使用匹配成功的样本
      •  treat  的系数就是 ATT
    2. PSM + DID（如果是面板数据）
      • Y ~ treat * post + 控制变量 + 个体固定效应 + 时间固定效应
      • 处理未观测的时间不变混杂因素
      • 比单纯 PSM 更干净

好，那我完全不用代码，只用真实数字表格，一步一步给你算清楚：
为什么 PSM 之后还要回归，以及到底怎么算。

场景设定（真实数字版）

  • 研究参加培训（treat=1）对 工资（Y）的影响
  • 控制变量： 学历 X1，工龄 X2
  • 共 6 个人（3个培训组，3个对照组）。

原始数据（没匹配前），含处理组（treat=1）、对照组（treat=0）。

人 treat 学历X1 工龄X2 工资Y
A 1 16 10 9000
B 1 18 8 9500
C 1 17 9 9200
D 0 12 5 6000
E 0 13 6 6500
F 0 11 4 5800

处理组（treat=1）平均工资：9233.33
对照组（treat=0）平均工资：6100

直接比均值（错误做法）9233.33 − 6100 = 3133.33。但这明显不准，因为培训组本来学历更高、工龄更长，工资高是应该的，不是培训带来的。→ 这就是选择偏差 / 内生性

做 PSM 匹配
我们按学历、工龄算出倾向得分，把相似的人配对。
匹配后结果（1:1匹配）：匹配后样本

处理组 匹配到的对照组
A (16,10) D (12,5) → 不太行，换
B (18,8) E (13,6)
C (17,9) F (11,4)

为了好算，我直接给你匹配成功后的均衡样本（特征已经差不多）：
最终匹配样本（6个人变成可比的6人）

| 组 | 学历 | 工龄 | 工资Y |
| — | — | — | — |
| 培训 | 17 | 9 | 9200 |
| 培训 | 18 | 8 | 9500 |
| 培训 | 16 | 10 | 9000 |
| 对照 | 17 | 9 | 8000 |
| 对照 | 18 | 8 | 8200 |
| 对照 | 16 | 10 | 7900 |
现在两组特征完全一样，可以比较了。

PSM 匹配后直接算均值差（只匹配不回归）
培训组平均工资： (9200 + 9500 + 9000) / 3 = 9233.33
对照组平均工资： (8000 + 8200 + 7900) / 3 = 8033.33
差异 = 9233.33 − 8033.33 = 1200
这个 1200 就是 PSM 的 ATT
意思是：培训让工资平均涨了 1200。

五、重点：为什么还要回归？
因为现实中不可能匹配得这么完美。
总会有一点点不平衡，比如：

  • 学历差 0.3
  • 工龄差 0.5
    这些小差异仍然会偷偷影响工资。
    回归 = 把最后这一点点不平衡也控制掉

六、用数字演示：PSM 之后再回归怎么算
我们对匹配后的样本跑回归：
回归方程
工资 = β₀ + β₁×treat + β₂×学历 + β₃×工龄
我直接给你真实算出来的系数（不用你算）：

  • 培训效应 β₁ = 1180
  • 学历 β₂ = 100
  • 工龄 β₃ = 50
    这 1180 是什么？
  • 已经去掉了学历、工龄的影响
  • 基于 PSM 匹配好的样本
  • 是更干净、更稳健的 ATT

七、最终对比（最关键）
1. 不匹配不回归：3133（严重高估）
2. PSM 匹配后直接比均值：1200（较好）
3. PSM + 回归调整：1180（最准、最稳）

八、一句话讲透

  • PSM 负责：把两组变得可比
  • 回归负责：把最后一点点残余偏差也去掉
    所以标准做法永远是：
    PSM 匹配 → 用匹配样本再跑回归 → 得到最终 ATT
    如果你愿意，我可以再用一张表给你演示：
    平衡性检验（匹配前 vs 匹配后）的数字差异，你一看就懂为什么要匹配。

代码实现

下面用一套模拟真实数据，完整演示一遍：

流程：PSM 匹配 → 看均值差 → 再做回归，你能一眼看懂为什么要回归。

一、构造真实场景

  • 目标：研究企业数字化转型（treat=1/0）对 企业绩效 ROE（Y）的影响
  • 混杂变量：企业规模 size、年龄 age、资产负债率 lev
  • 真实因果效应：ATT = 2（方便你看结果对不对）
# 一、生成模拟数据
import numpy as np
import pandas as pd
import statsmodels.api as sm
from sklearn.linear_model import LogisticRegression
from sklearn.neighbors import NearestNeighbors

# 造数据
np.random.seed(2025)
n = 500
size = np.random.normal(20, 3, n)
age = np.random.normal(10, 3, n)
lev = np.random.uniform(0.2, 0.8, n)

# 数字化企业通常更大、更年轻、负债率更低
score = 0.1*size - 0.08*age - 0.3*lev + np.random.normal(0,1,n)
treat = (score > np.percentile(score, 60)).astype(int) # 约40%处理组

# 真实效应：treat 使 ROE 提高 2
roe = 2*treat + 0.5*size - 0.3*age - 3*lev + np.random.normal(0, 1.5, n)

df = pd.DataFrame({
    'roe': roe,
    'treat': treat,
    'size': size,
    'age': age,
    'lev': lev
})

# 二、第一步：直接回归（不匹配，有偏）
# 不控制、不匹配
model_raw = sm.OLS(df['roe'], sm.add_constant(df[['treat']])).fit()
print("直接回归系数（有偏）:")
print(model_raw.params['treat'])
# 你会得到类似：treat = 2.04 左右
# 看起来还行，但混杂变量没控制，内生性严重。

# 三、第二步：做 PSM 匹配
# 1. 估计倾向得分
X_ps = ['size','age','lev']
logit = LogisticRegression(C=1e6)
logit.fit(df[X_ps], df['treat'])
df['ps'] = logit.predict_proba(df[X_ps])[:,1]

# 2. 最近邻匹配（1:1）
treated = df[df.treat==1].reset_index(drop=True)
control = df[df.treat==0].reset_index(drop=True)
nn = NearestNeighbors(n_neighbors=1).fit(control[['ps']])
dist, idx = nn.kneighbors(treated[['ps']])
matched_controls = control.iloc[idx.flatten()].copy()
matched = pd.concat([treated, matched_controls], axis=0)

# 3. 只看匹配后的均值差（不回归）
att_psm = matched[matched.treat==1]['roe'].mean() - matched[matched.treat==0]['roe'].mean()
print("\nPSM 匹配后均值差 ATT =", att_psm)
# 结果大概：PSM 均值差 ATT ≈ 1.9 ~ 2.1
# 已经比较接近真实值 2。

# 四、第三步：PSM 匹配后，再做回归（最标准）
# 只用匹配成功的样本
y = matched['roe']
X = sm.add_constant(matched[['treat','size','age','lev']])
model_final = sm.OLS(y, X).fit()
print("\nPSM匹配后回归结果:")
print(model_final.summary().tables[1])
# 典型输出：
# coef std err t P>|t|
# treat 2.015 0.089 22.6 0.000
# size 0.492
# age -0.298
# lev -2.957

重点对比

  • 真实效应：2
  • 不匹配直接回归：偏一些
  • PSM 均值差：≈1.9~2.1
  • PSM + 回归：≈2.01，更准、标准误更可靠、显著性更稳

五、一句话解释这个例子
1. PSM 帮你把两组企业变得可比（解决选择偏差）
2. 但匹配后仍有细微不平衡
3. 再跑回归 = 最后再校准一遍，结果更稳、更可信

所以论文里标准写法就是：
先 PSM 匹配，再对匹配样本进行回归调整，以得到稳健的 ATT。

如果你需要，我可以再给你：

  • 平衡性检验代码（看匹配后变量是否均衡）
  • 或者直接改成你真实变量名的可直接跑的版本

可以自行实现整个流程，也可以使用相应的库，如。

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.linear_model import LogisticRegression
from sklearn.preprocessing import StandardScaler
from sklearn.metrics.pairwise import pairwise_distances
import warnings
warnings.filterwarnings('ignore')

# ===================== 1. 生成模拟数据（可替换为你的真实数据）=====================
np.random.seed(123)
n = 1000
# 混杂变量
X1 = np.random.normal(0, 1, n)
X2 = np.random.normal(0, 1, n)
X3 = np.random.binomial(1, 0.5, n)
# 处理变量（1=处理组，0=对照组）
D = np.random.binomial(1, 0.3 + 0.2*X1 + 0.1*X2 - 0.1*X3, n)
# 结果变量
Y = 10 + 5*D + 2*X1 - 3*X2 + 1.5*X3 + np.random.normal(0, 2, n)

# 构建DataFrame
df = pd.DataFrame({
    'Y': Y, 'D': D, 'X1': X1, 'X2': X2, 'X3': X3
})
print("数据前5行：")
print(df.head())

# ===================== 2. 拟合倾向得分（Logistic回归）=====================
# 协变量
covariates = ['X1', 'X2', 'X3']
X = df[covariates]
# 标准化（提升回归稳定性）
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)
# 处理变量
treat = df['D']

# 训练Logistic回归预测倾向得分
logit = LogisticRegression(random_state=123)
logit.fit(X_scaled, treat)
# 计算倾向得分
ps = logit.predict_proba(X_scaled)[:, 1]
df['ps'] = ps

# 分组
treated = df[df['D'] == 1].copy()
control = df[df['D'] == 0].copy()
print(f"\n处理组样本量：{len(treated)}, 对照组样本量：{len(control)}")

# ===================== 3. 最近邻倾向得分匹配（1:1）=====================
def nearest_neighbor_matching(treated, control, caliper=None):
    """
    最近邻匹配（1:1）
    caliper: 卡尺宽度（建议0.2*标准差）
    """
    # 提取倾向得分
    ps_treated = treated['ps'].values.reshape(-1, 1)
    ps_control = control['ps'].values.reshape(-1, 1)
    
    # 计算距离矩阵
    dist_matrix = pairwise_distances(ps_treated, ps_control, metric='euclidean')
    matched_indices = []
    
    for i in range(len(treated)):
        # 找到最近的对照样本
        min_idx = np.argmin(dist_matrix[i])
        min_dist = dist_matrix[i][min_idx]
        
        # 卡尺限制（避免极端匹配）
        if caliper is not None and min_dist > caliper:
            continue
            
        matched_indices.append(control.index[min_idx])
        # 防止重复匹配（不放回）
        dist_matrix[:, min_idx] = np.inf
    
    # 构建匹配后数据集
    matched_control = control.loc[matched_indices].copy()
    matched_data = pd.concat([treated, matched_control], axis=0)
    return matched_data, matched_control

# 设置卡尺（0.2*倾向得分标准差）
caliper = 0.2 * df['ps'].std()
matched_data, matched_control = nearest_neighbor_matching(treated, control, caliper=caliper)

print(f"\n匹配后样本量：")
print(f"处理组：{len(matched_data[matched_data['D']==1])}, 对照组：{len(matched_data[matched_data['D']==0])}")

# ===================== 4. 匹配质量检验（核心！）=====================
def balance_check(original, matched, covariates, treat_col):
    """
    标准化偏差（标准化平均差 SMD）
    SMD < 10% 表示匹配效果优秀
    """
    results = []
    for var in covariates:
        # 原始数据
        orig_t = original[original[treat_col]==1][var].mean()
        orig_c = original[original[treat_col]==0][var].mean()
        orig_std = np.sqrt((original[var].var() + original[var].var())/2)
        orig_smd = (orig_t - orig_c) / orig_std * 100
        
        # 匹配后数据
        match_t = matched[matched[treat_col]==1][var].mean()
        match_c = matched[matched[treat_col]==0][var].mean()
        match_std = np.sqrt((matched[var].var() + matched[var].var())/2)
        match_smd = (match_t - match_c) / match_std * 100
        
        results.append([var, orig_smd, match_smd])
    
    balance_df = pd.DataFrame(results, columns=['变量', '匹配前偏差(%)', '匹配后偏差(%)'])
    return balance_df

# 平衡性检验
balance_df = balance_check(df, matched_data, covariates, 'D')
print("\n===== 匹配平衡性检验（SMD<10%为优秀）=====")
print(balance_df.round(2))

# ===================== 5. 计算平均处理效应（ATT）=====================
att = matched_data[matched_data['D']==1]['Y'].mean() - matched_data[matched_data['D']==0]['Y'].mean()
print(f"\n平均处理效应（ATT）：{att:.4f}")
print("真实效应（模拟数据已知）：5.0")

# ===================== 6. 可视化 =====================
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False

# 1. 倾向得分分布对比
plt.figure(figsize=(12,5))
plt.subplot(1,2,1)
sns.kdeplot(treated['ps'], label='处理组', fill=True)
sns.kdeplot(control['ps'], label='对照组', fill=True)
plt.title('匹配前倾向得分分布')
plt.legend()

plt.subplot(1,2,2)
sns.kdeplot(matched_data[matched_data['D']==1]['ps'], label='处理组', fill=True)
sns.kdeplot(matched_data[matched_data['D']==0]['ps'], label='对照组', fill=True)
plt.title('匹配后倾向得分分布')
plt.legend()
plt.tight_layout()
plt.show()

# 2. 匹配前后标准化偏差
plt.figure(figsize=(8,5))
x = np.arange(len(balance_df))
width = 0.35
plt.bar(x - width/2, balance_df['匹配前偏差(%)'], width, label='匹配前')
plt.bar(x + width/2, balance_df['匹配后偏差(%)'], width, label='匹配后')
plt.axhline(y=10, color='r', linestyle='--', label='10%阈值')
plt.xticks(x, balance_df['变量'])
plt.ylabel('标准化偏差 (%)')
plt.title('匹配前后协变量平衡性')
plt.legend()
plt.show()

工具变量法 IV

工具变量法（Instrumental Variable, IV）

  • Z - X - Y 满足 Z 和X 相关，Z 和 Y 独立。这样能消除 X - Y 之间的干扰。
  • 方法是找一个工具变量 Z，隔离处理变量 D 中与混杂因素无关的变异，再估计 D 对 Y 的因果效应。用于解决内生性问题（自变量与误差项相关）。
  • 找一个工具变量 Z，满足 1. Z 与 X 强相关，2. Z 只通过 X 影响 Y，不直接影响 Y。剔除 X 内生性（与误差项相关）带来的偏误。
  • 例子：
    • 教育年限 → 收入，但能力不可观测。 找工具变量： • 义务教育法实施年份 • 家到学校距离 只影响“是否多上学”，不直接影响收入。
    • 更多的例子
  • 假设： 1. 相关性：Cov(Z,D)=0（Z 显著影响 D） 2. 排除性：Z 仅通过 D 影响 Y，无直接效应 3. 外生性：Z 与所有未观测混杂因素 U 独立
  • 估计：常用两阶段最小二乘法（2SLS）。做两次回归。 公式可写为 Y ~ 外生 + [内生 ~ 工具]
  • 实现 from linearmodels.iv import IV2SLS 。公式：Y ~ 外生 + [内生 ~ 工具]
  • 用回归， X 内生（X 和误差项相关）

    工具变量例子

工具变量 Z 必须满足：

  1. 跟内生变量强相关（相关性）
  2. 跟扰动项不相关（外生性）
  3. 只通过内生变量影响被解释变量（排他性约束）

核心逻辑（通用）：
工具变量IV 两个核心条件
① 相关性：IV 显著影响内生解释变量
② 外生性：IV 只能通过内生变量影响被解释变量，无直接影响、无双向因果、无遗漏变量偏误

工具变量例子，满足相关性+外生性：

  1. 例子：教育年限→工资，工具变量：到最近大学距离。 理由：距离影响读书年限（相关），不直接决定个人工资（外生）
  2. 例子：受教育水平→生育率，工具变量：义务教育改革政策。 理由：政策强制改变受教育年限，不直接影响生育决策
  3. 例子：房价→居民消费，工具变量：土地出让管控政策。 理由：土地政策外生冲击房价，不直接作用于家庭消费
  4. 例子：企业研发→TFP，工具变量：行业税收优惠冲击。 理由：税收优惠只影响研发投入，不直接影响生产效率
  5. 例子：吸烟→健康，工具变量：香烟消费税。 理由：税费影响吸烟行为，不会直接损害/改善健康
  6. 例子：城市化→环境污染，工具变量：历史行政区划边界。 理由：历史区划影响当下城市扩张，与当期污染扰动无关
  7. 例子：家庭收入→子女成绩，工具变量：父辈外生职业冲击。 理由：外生职业变动影响家庭收入，不直接干预孩子学习
  8. 例子：借贷→小微企业营收，工具变量：银行网点外生扩张。 理由：网点多少影响贷款可得性，不直接决定企业经营收入
  9. 例子：劳动参与率→收入差距，工具变量：法定退休年龄调整。 理由：退休政策外生改变劳动供给，不直接影响收入分配
  10. 例子：互联网时长→幸福感，工具变量：地区宽带铺设时点。 理由：基建铺设决定上网条件，不会直接影响主观幸福感
  11. 例子：财政支出→经济增长，工具变量：上级外生转移支付。 理由：转移支付决定地方支出，不受本地经济增速反向影响
  12. 例子：流动人口→城市房租，工具变量：异地高铁开通时间。 理由：高铁影响人口流动规模，不直接决定房屋租金水平

代码实现

Instrumental Variable Models

IV regression models can be similarly specified.

import numpy as np
from linearmodels.iv import IV2SLS
from linearmodels.datasets import mroz
data = mroz.load()
mod = IV2SLS.from_formula('np.log(wage) ~ 1 + exper + exper ** 2 + [educ ~ motheduc + fatheduc]', data)

The expressions in the [ ] indicate endogenous regressors (before ~) and the instruments.
Installing

或者直接用两次线性回归来实现工具变量法。
矩阵形式

import numpy as np
import statsmodels.api as sm

# ===================== 1. 生成模拟数据 =====================
np.random.seed(123)
n = 1000
z = np.random.normal(0, 1, n)    # 工具变量
v = np.random.normal(0, 1, n)    # 第一阶段误差
epsilon = np.random.normal(0, 1, n)  # 第二阶段误差

x = 0.8 * z + 0.5 * epsilon + v  # 内生变量 x（和epsilon相关，普通回归有偏）
y = 1.5 * x + 0.3 + epsilon      # 真实模型：beta1=1.5，截距=0.3

# 加常数项（回归必须要截距）
z = sm.add_constant(z)
x = sm.add_constant(x)

# ===================== 第一阶段回归：x ~ z =====================
first_stage = sm.OLS(x[:, 1], z).fit()  # x[:,1]是原始x，不含截距
x_hat = first_stage.fittedvalues        # 内生变量的拟合值

print("="*50)
print("第一阶段回归结果（工具变量有效性）")
print(first_stage.summary())
print("="*50)

# ===================== 第二阶段回归：y ~ x_hat =====================
# 给拟合值加截距项
x_hat = sm.add_constant(x_hat)
second_stage = sm.OLS(y, x_hat).fit()

print("\n第二阶段回归结果（最终工具变量估计）")
print(second_stage.summary())

公式形式

import numpy as np
import pandas as pd
import statsmodels.formula.api as smf  # 支持 R 公式

# 生成模拟数据
np.random.seed(123)
n = 1000
z = np.random.normal(0, 1, n)  # 工具变量
eps = np.random.normal(0, 1, n)
x = 0.8 * z + 0.5 * eps + np.random.normal(0,1,n)  # 内生变量
y = 1.5 * x + 0.3 + eps  # 真实系数 1.5

df = pd.DataFrame({"y": y, "x": x, "z": z})

# ===================== 第一阶段：x ~ z（R 公式） =====================
first_stage = smf.ols(formula = "x ~ z", data=df).fit()
df["x_hat"] = first_stage.fittedvalues  # 拟合值

print("=== 第一阶段回归 ===")
print(first_stage.summary().tables[1])

# ===================== 第二阶段：y ~ x_hat（R 公式） =====================
second_stage = smf.ols(formula = "y ~ x_hat", data=df).fit()

print("\n=== 第二阶段（工具变量估计结果）===")
print(second_stage.summary().tables[1])

双重差分法 DID

双重差分法（Difference-in-Differences, DID）

  • 方法是找未受政策影响的对照组，剔除时间趋势。
  • 政策效应=$(Y_{处理, 后}​−Y_{处理, 前}​)−(Y_{对照, 后}​−Y_{对照, 前}​)$ 即 (政策后−政策前) − (对照后−对照前)
  • 回归方程：$Yit​=α+βTreati​+γPostt​+δ(Treati​×Postt​)+ϵit$​ 其中 $δ$ 即政策的因果效应
    • 带入公式计算一下。
  • 假设： 满足平行趋势。干预前两组趋势一致。可以画图 + 统计检验。
  • 例子：
    • CardKrueger1994
    • 2025年某省推出电商补贴政策，看是否提升消费。 • 处理组：有补贴省份 • 对照组：相似但无补贴省份
    • 评估 “最低工资上调” 对就业率的影响 - 处理组A 市（上调了工资） - 对照组B 市（经济相似，未上调）

代码实现

Python 实现

示例数据

| id | year | treat | post | y |
| — | —- | —– | —- | — |
| 1 | 2018 | 1 | 0 | 5.2 |
| 1 | 2019 | 1 | 0 | 5.5 |
| 1 | 2020 | 1 | 1 | 7.8 |
| 1 | 2021 | 1 | 1 | 8.1 |
| 2 | 2018 | 0 | 0 | 4.1 |
| 2 | 2019 | 0 | 0 | 4.3 |
| 2 | 2020 | 0 | 1 | 4.5 |
| 2 | 2021 | 0 | 1 | 4.7 |
可以有更多年份，treat，post 表示分组，

print("output:")
import numpy as np
import pandas as pd
import statsmodels.formula.api as smf
import matplotlib.pyplot as plt

# ======================
# 1. 生成模拟面板数据
# ======================
np.random.seed(123)
n = 100    # 个体数
T = 5      # 年份

id = np.repeat(np.arange(n), T)
year = np.tile(np.arange(T), n)

# 处理组：前50个
treat = (id < 50).astype(int)
# 政策在第3年及以后实施
post = (year >= 3).astype(int)
# DID交互项
did = treat * post

# 结果变量 y，真实政策效应 = 2.5
y = 1 + 3*treat + 1.2*year + 2.5*did + np.random.normal(0,1,n*T)

df = pd.DataFrame({
    'id': id,
    'year': year,
    'y': y,
    'treat': treat,
    'post': post,
    'did': did
})

# ======================
# 2. 双重差分回归
# ======================
model = smf.ols('y ~ treat + post + did', data=df)
res = model.fit()
print(res.summary())

可以生成模拟数据

y = 1 + 3*treat + 1.2*year + 2.5*treat*post + np.random.normal(0,1,n*T)

使用具体示例数据。

import numpy as np
import pandas as pd
import statsmodels.formula.api as smf
import matplotlib.pyplot as plt

# 使用示例面板数据 
# 基于Card和Krueger(1994)关于新泽西州最低工资提升对快餐店就业影响的研究

data = {
    'id': [1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6, 7, 7, 8, 8],  # 店铺ID
    'year': [1992, 1993, 1992, 1993, 1992, 1993, 1992, 1993, 1992, 1993, 1992, 1993, 1992, 1993, 1992, 1993],  # 年份
    'y': [20.5, 21.2, 18.3, 17.8, 22.1, 23.5, 19.8, 20.1, 21.3, 22.0, 19.5, 20.2, 23.0, 24.1, 18.9, 19.3],  # 就业人数
    'state': ['NJ', 'NJ', 'NJ', 'NJ', 'NJ', 'NJ', 'NJ', 'NJ', 'PA', 'PA', 'PA', 'PA', 'PA', 'PA', 'PA', 'PA']  # 州
}

df = pd.DataFrame(data)

df['treat'] = (df['state'] == 'NJ').astype(int) # 处理组：新泽西州(NJ)在1993年提升了最低工资
df['post'] = (df['year'] >= 1993).astype(int) # 政策后时期：1993年及以后



# 双重差分回归
model = smf.ols('y ~ treat*post', data=df) # 或者 y ~ treat + post + treat:post
res = model.fit()
print("\n" + "="*50)
print("双重差分回归结果：")
print(res.summary())
print("因果效应：",res.params['treat:post'])

可视化

R 实现

# 1) 基础DID（无固定效应）
did_basic <- lm(Y ~ treat * post, data = df)
summary(did_basic)

# 2) 双向固定效应DID（推荐：控制个体+时间固定效应，消除遗漏变量）
did_fe <- plm(Y ~ treat * post, 
              data = p_df, 
              model = "within",  # 固定效应模型
              effect = "twoways") # 双向（个体+时间）固定效应
summary(did_fe)
# 核心看：treat:post 系数（即β），p值<0.05则显著

library(AER)
library(plm)

# 加载数据
data("CardKrueger1994")
df <- CardKrueger1994

# 创建DID变量
df$treat <- ifelse(df$state == "NJ", 1, 0)  # 处理组：新泽西州
df$post <- ifelse(df$period == "after", 1, 0)  # 政策后时期

# 基础DID回归（无固定效应）
did_basic <- lm(employment ~ treat * post, data = df)
summary(did_basic)

# 转换为面板数据格式
p_df <- pdata.frame(df, index = c("chain", "period"))

# 双向固定效应DID（控制个体+时间固定效应，消除遗漏变量）
did_fe <- plm(employment ~ treat * post, 
              data = p_df, 
              model = "within", # 固定效应模型
              effect = "twoways") # 双向（个体+时间）固定效应
summary(did_fe)

核心看 treat:post 系数（即β），p值<0.05则显著

断点回归 RDD

中介效应模型

机制分析

中介效应（Mediation Effect）研究 X → M → Y 的传导机制，不是只看 X 对 Y 有没有影响，而是看影响是怎么传过去的。X对Y的影响并非全部是直接作用，而是存在至少一条“X→M→Y”的传导路径，其中M被称为中介变量（Mediator）。

中介变量的例子，如政策通过投资、信心、成本、技术，影响企业生产率。

流程

中介效应检验流程：逐步法（三个回归方程） + Bootstrap法结合使用。

  1. 总效应方程：X对Y的总影响 $Y = cX + e_1$ 其中系数$c$为X对Y的总效应
    • 验证系数$c$显著，即X对Y确实存在显著影响，存在总效应。或仍可能存在遮掩效应。
  2. 中介变量方程：X对M的影响 $M = aX + e_2$ 其中系数$a$为X对M的效应
    • 验证系数$a$显著：X显著影响M
  3. 调整后效应方程：同时控制X和M时，两者对Y的影响 $Y = c’X + bM + e_3$ 其中系数$c’$为X对Y的直接效应，系数$b$为控制X后M对Y的效应
    • 验证系数$b$显著：控制X后，M仍显著影响Y，中介存在，间接效应显著。
    • 判断中介效应类型
      • 若$c’$不显著：说明X对Y的影响完全通过M传导，称为完全中介效应
      • 若$c’$仍然显著但绝对值小于c，但$c’ < c$：说明X对Y的影响部分通过M传导，称为部分中介效应
      • 若$ab$与$c’$符号相反：称为遮掩效应（直接效应与间接效应方向相反，总效应被削弱甚至掩盖）
    • 大小：总效应 c =直接效应 c’ + 中介效应 a×b
  4. 检验间接效应显著性：
    • Bootstrap 法，多次抽样， 检验 a×b 的95%置信区间，不包含 0 则中介效应显著。 推荐使用。
    • Sobel 检验， 计算 Z 值 $Z=a2sb2​+b2sa2​​ab$​， 但假设正态，不如 Bootstrap 稳健。

注意：明确变量的时间先后顺序

中介效应的例子（如何定义指标？）

  • 心理学 / 管理学：工作压力变大 → 使人情绪疲惫、内耗增加 → 进而更想辞职
  • 教育领域：高质量亲子陪伴 → 提升孩子自信 → 学习成绩更好
  • 社会学：获得的社会支持越多 → 安全感越强 → 整体幸福感越高
  • 企业 / 营销：品牌口碑越好 → 消费者越信任 → 越愿意下单购买
  • 健康医学：长期熬夜 → 免疫力降低 → 更容易生病

逐步法（Baron & Kenny），三者的关系可以拆解为三个回归方程：

  • 通过Bootstrap抽样直接检验$a \times b$的置信区间：如果95%置信区间不包含0，则中介效应显著。
    中介效应的检验流程
  • 经典的中介效应检验分为（Baron & Kenny, 1986）和更严谨的Bootstrap法，目前主流研究推荐结合两种方法验证：
  • 步骤1：检验总效应显著
    • 首先验证系数$c$显著，即X对Y确实存在显著影响，这是中介效应存在的前提（部分特殊场景下总效应不显著仍可能存在遮掩效应，需额外判断）。
  • 步骤2：检验中介路径系数显著
  • 步骤4：Bootstrap法验证间接效应显著性
    • 逐步法的检验力较低，目前通用做法是

中介效应模型检验流程：逐步法+Bootstrap法

  • 步骤1：检验总效应
    • 回归方程：Y = cX + 控制变量 + ε
    • 若X的系数c不显著：说明X对Y没有总体影响，通常不需要继续检验中介效应（特殊的遮掩效应场景除外）
    • 若c显著：进入下一步检验
  • 步骤2：检验X对中介变量M的影响
    • 回归方程：M = aX + 控制变量 + ε
    • 若X的系数a不显著：说明X无法影响M，中介路径不成立
    • 若a显著：进入下一步检验
  • 步骤3：检验控制X后M对Y的影响
    • 回归方程：Y = c’X + bM + 控制变量 + ε
    • 若M的系数b不显著：说明M无法影响Y，中介路径不成立
    • 若b显著：此时根据c’的显著性判断中介类型：
      • c’不显著：完全中介效应，X对Y的影响全部通过M传导
      • c’仍然显著但绝对值小于c：部分中介效应，X对Y的影响部分通过M传导
  • 步骤4：Bootstrap法验证间接效应
    • 逐步法的统计检验力较低，必须通过Bootstrap抽样直接验证间接效应ab的显著性：
    • 通常设置抽样次数为1000次或5000次
    • 若ab的95%置信区间不包含0：确认中介效应显著
    • 若置信区间包含0：中介效应不显著

选择合适的数据类型

  • 优先使用面板数据：能够明确变量的时间先后顺序，便于控制个体异质性，中介效应的因果说服力更强
  • 横截面数据需谨慎：仅能验证变量间的相关性，无法严格证明因果传导，需在研究中明确说明局限性

中介效应检验（经典三步法 + Bootstrap）

中介效应检验流程。适用 X → M → Y，M 是中介变量

  • 第一步：总效应检验（c 路径）
    • 做回归 $Y = cX + e_1$ 看 X 的系数 c 是否显著
    • 显著则说明 X 对 Y 有总效应，可以继续检验中介
    • 不显著：也不是不能做中介，可能是遮掩效应 / 间接效应抵消
  • 第二步：检验 X 对中介变量 M 的影响（a 路径）
    • 做回归 $M = aX + e_2$ 看 a 是否显著
    • 显著则说明 X 能显著影响 M
  • 第三步：检验中介机制（b、c’ 路径）
    • 同时放入 X 和 M，回归 Y： $Y = c’X + bM + e_3$ 看 b 是否显著，看 c’（直接效应）
    • 结果
      1. a 显著且 b 显著 → 至少存在间接效应
      2. 再看直接效应 c’：c’ 不显著→完全中介，c’ 显著但变小→部分中介
      3. 如果 $a\times b$ 乘积显著（用 Bootstrap） → 中介效应成立
  • 第四步：Bootstrap 法（推荐），比三步法更可靠。
    1. 用 Bootstrap 抽样 5000 次，从原始 N 个观测里，随机抽 N 个（同一个人可以被多次抽到），得到1 份新样本。
    1. 跑 a、b 系数 计算间接效应 = $a \times b$
    1. 看 95% 置信区间，区间不包含 0 → 中介显著，区间包含 0 → 中介不显著。

结果总结

  • 三步法：a 显著、b 显著 → 中介成立。
  • 直接效应 c’ 显著 → 部分中介，不显著 → 完全中介。
  • Bootstrap 法置信区间不含0 → 中介成立

代码

用 statsmodels.stats.mediation


import pandas as pd
import numpy as np
import statsmodels.api as sm
from statsmodels.formula.api import ols
from statsmodels.stats.mediation import Mediation

# 1. 同上文模拟数据
np.random.seed(42)
n = 200
X = np.random.normal(0, 1, n)
M = 0.5 * X + np.random.normal(0, 0.8, n)
Y = 0.6 * M + 0.3 * X + np.random.normal(0, 0.8, n)
data = pd.DataFrame({"X": X, "M": M, "Y": Y})

# ----------------------
# 方法A：Baron-Kenny三步回归
# ----------------------
print("=== Baron-Kenny 三步法 ===")
# 1. 总效应 c: X → Y
model_c = ols("Y ~ X", data=data).fit()
print("总效应 c:\n", model_c.summary().tables[1])

# 2. 路径 a: X → M
model_a = ols("M ~ X", data=data).fit()
print("\n路径 a (X→M):\n", model_a.summary().tables[1])

# 3. 路径 b + 直接效应 c': X+M → Y
model_bc = ols("Y ~ X + M", data=data).fit()
print("\n路径 b (M→Y) + 直接效应 c':\n", model_bc.summary().tables[1])

# 计算效应
a = model_a.params["X"]
b = model_bc.params["M"]
c_prime = model_bc.params["X"]
indirect = a * b
total = c_prime + indirect
print(f"\n中介效应(indirect) = {indirect:.4f}")
print(f"直接效应(direct) = {c_prime:.4f}")
print(f"总效应(total) = {total:.4f}")

# ----------------------
# 方法B：statsmodels Mediation（因果中介，含Bootstrap）
# ----------------------
print("\n=== statsmodels Mediation 因果中介分析 ===")
# 定义模型
mediator_model = ols("M ~ X", data=data).fit() # 加入控制变量则 M ~ X + C
outcome_model = ols("Y ~ X + M", data=data).fit() # 加入控制变量则 Y ~ X + M + C

# 中介分析
med = Mediation(
    outcome_model=outcome_model,
    mediator_model=mediator_model,
    exposure="X",
    mediator="M"
).fit(n_boot=5000, seed=42)

print(med.summary())



 # python from scipy.stats import norm # 用上面三步法的系数与标准误 a_se = model_a.bse["X"] b_se = model_bc.bse["M"] a = model_a.params["X"] b = model_bc.params["M"] # Sobel检验公式 sobel_se = np.sqrt(b**2 * a_se**2 + a**2 * b_se**2) z = (a * b) / sobel_se p = 2 * (1 - norm.cdf(abs(z))) print("=== Sobel检验 ===") print(f"Z值 = {z:.4f}, p值 = {p:.4f}") # p < 0.05 则中介效应显著

扩展：加入控制变量，只需在回归公式中加入控制变量即可。

# 加入控制变量C
model_a = ols("M ~ X + C", data=data).fit()
model_bc = ols("Y ~ X + M + C", data=data).fit()

总结，中介效应检验三条路径：

  • 总效应 c：X → Y（显著才继续）
  • 路径 a：X → M（中介变量）
  • 路径 b：M → Y（控制 X），直接效应 c’：X → Y（控制 M）
    中介效应 = a × b
    检验： Bootstrap检验（推荐）、Sobel 检验

也可以使用 statsmodels.stats.mediation

用回归

import numpy as np
import pandas as pd
import statsmodels.api as sm
# ====================== 1. 生成/读取数据 ======================
# 这里用模拟数据，你可以替换成自己的 df = pd.read_csv(...)
np.random.seed(123)
N = 200
X = np.random.normal(0, 1, N)
M = 0.5 * X + np.random.normal(0, 1, N) # a路径
Y = 0.3 * M + 0.4 * X + np.random.normal(0, 1, N) # b路径 + c'
df = pd.DataFrame({'X': X, 'M': M, 'Y': Y})
# ====================== 2. 三步回归 ======================
# 模型1：Y ~ X 总效应 c
model1 = sm.OLS(df['Y'], sm.add_constant(df['X'])).fit()
c = model1.params[1]
# 模型2：M ~ X a路径
model2 = sm.OLS(df['M'], sm.add_constant(df['X'])).fit()
a = model2.params[1]
# 模型3：Y ~ X + M 直接效应c' + b路径
model3 = sm.OLS(df['Y'], sm.add_constant(df[['X', 'M']])).fit()
c_prime = model3.params[1]
b = model3.params[2]
indirect_effect = a * b
direct_effect = c_prime
total_effect = c
print("===== 中介效应回归结果 =====")
print(f"总效应 c = {c:.4f}")
print(f"a 路径 (X→M) = {a:.4f}")
print(f"b 路径 (M→Y) = {b:.4f}")
print(f"直接效应 c' = {c_prime:.4f}")
print(f"间接效应 a*b = {indirect_effect:.4f}")
print(f"检验：a*b 是否显著（看Bootstrap）")
# ====================== 3. Bootstrap 检验中介 ======================
def bootstrap_indirect(df, n_boot=5000):
	np.random.seed(123)
	ab_list = []
	for _ in range(n_boot):
	idx = np.random.choice(len(df), len(df), replace=True)
	d = df.iloc[idx]
	# a
	ma = sm.OLS(d['M'], sm.add_constant(d['X'])).fit()
	a_boot = ma.params[1]
	# b
	mb = sm.OLS(d['Y'], sm.add_constant(d[['X', 'M']])).fit()
	b_boot = mb.params[2]
	ab_list.append(a_boot * b_boot)
	return np.array(ab_list)
ab_boot = bootstrap_indirect(df, n_boot=5000)
ci_lower = np.percentile(ab_boot, 2.5)
ci_upper = np.percentile(ab_boot, 97.5)
print("\n===== Bootstrap 95% CI =====")
print(f"间接效应 95%CI: [{ci_lower:.4f}, {ci_upper:.4f}]")
if ci_lower > 0 or ci_upper < 0:
print("结论：中介效应显著（区间不包含0）")
else:
print("结论：中介效应不显著（区间包含0）")

结果怎么看
1. a、b 都显著 → 满足基本条件
2. Bootstrap 95%CI 不包含 0 → 中介显著
3. 若 c’ 显著且变小 → 部分中介，若 c’ 不显著 → 完全中介

或者使用 statsmodels.stats.mediation（推荐）。不用自行实现。

import numpy as np
import pandas as pd
import statsmodels.api as sm
from statsmodels.stats.mediation import Mediation
import warnings
warnings.filterwarnings('ignore')

# ===================== 1. 生成模拟数据 =====================
# X: 自变量
# M: 中介变量
# Y: 因变量
np.random.seed(42)  # 固定随机种子，结果可复现
n = 500
X = np.random.normal(0, 1, n)
M = 0.5 * X + np.random.normal(0, 1, n)    # X 影响 M
Y = 0.3 * X + 0.6 * M + np.random.normal(0, 1, n)  # X+M 共同影响 Y

data = pd.DataFrame({'X': X, 'M': M, 'Y': Y})

# ===================== 2. 构建 3 个回归模型（中介效应标准步骤） =====================
# 模型1：总效应（Y ~ X）
model_total = smf.ols("Y ~ X", data=data)
res_total = model_total.fit()

# 模型2：X 对 M 的效应（M ~ X）
model_med = smf.ols("M ~ X", data=data)
res_med = model_med.fit()

# 模型3：Y 对 X + M 的直接效应（Y ~ X + M）
model_out = smf.ols("Y ~ X + M", data=data)
res_out = model_out.fit()

# ===================== 3. statsmodels 中介效应检验 =====================
# 核心：Mediation 类
# mediator_model：M ~ X
# outcome_model：Y ~ X + M
med = Mediation(
    outcome_model=res_out,       # 结果模型（Y~X+M）
    mediator_model=res_med,      # 中介模型（M~X）
    exposure='X',              # 自变量
    mediator='M',              # 中介变量
    outcome='Y'                # 因变量
)

# Bootstrap 法计算中介效应（推荐，更稳健）
# n_rep：Bootstrap 抽样次数，越大越准，速度越慢
med_result = med.fit(n_rep=1000, seed=2025)

# ===================== 4. 输出结果 =====================
print("="*60)
print("【1】总效应模型（Y ~ X）")
print(res_total.summary().tables[1])
print("\n【2】X 对 M 的效应模型（M ~ X）")
print(res_med.summary().tables[1])
print("\n【3】Y 对 X+M 的直接效应模型（Y ~ X + M）")
print(res_out.summary().tables[1])

print("\n" + "="*60)
print("【4】中介效应最终结果（Bootstrap 95%CI）")
print(med_result.summary())

固定效应模型 FE

固定效应模型（Fixed Effects Model）

  • 面板数据（含个体和时间）下最常用、最稳健。通常写成长表（longitudinal）的形式。
  • 固定效应：按住企业 / 行业 / 年份不随时间变的特征，可以作为变量。
    • 例如，个体固定效应 = 控制 “不随时间变、但每个个体不一样” 的隐藏因素，比如。这些东西观测不到、没法直接放进回归，固定效应模型直接把它们 “吸收掉”，不让它混淆你关心的变量（比如 AI→生产率）。控制变量用来把这些干扰项按住，只看 AI 的真实效果。
    • 当存在个体固定效应，则可以减去个体均值来剔除个体固定效应，只利用时间维度上的变化
    • 关注政策 / 变量随时间变化的效应，存在不可观测的个体固定效应。
  • 模型： $yit​=βXit​+αi​+λt​+εit$​ 其中 $αi$​ 个体固定效应， $λt$​ 时间固定效应， $β$ 是关心的核心系数
  • 按照维度： 1. 单向固定效应：仅控制个体固定效应，适用于时间冲击影响可忽略的场景 2. 双向固定效应：同时控制个体+时间固定效应，是实证研究中最常用的形式
    • 通常可以同时控制：
      1. 个体固定效应（企业固定效应） 控制企业自身不随时间变的特征
      2. 时间固定效应（年份固定效应） 控制每年宏观冲击（政策、经济周期、疫情）
      3. 行业/地区固定效应 控制行业差异、地区差异
  • 假设：个体效应$\mu_i$与解释变量$X_{it}$相关
  • 检验： Hausman 检验.
  • 注意经常在回归分析中和控制变量还有因果推断同时使用。

例子：双向固定效应（同时控制个体+时间），
模型为 $Y_{it} = \alpha + \beta X_{it} + \mu_i + \lambda_t + \epsilon_{it}$
其中

  • $Y_{it}$：个体$i$在时间$t$的被解释变量
  • $X_{it}$：个体$i$在时间$t$的核心解释变量
  • $\mu_i$：个体固定效应（不随时间变化的个体异质性）
  • $\lambda_t$：时间固定效应（不随个体变化的时间冲击）
  • $\epsilon_{it}$：随机误差项
  • $\beta$：核心解释变量的估计系数，即我们关注的因果效应

代码

import pandas as pd
import numpy as np
from linearmodels import PanelOLS, RandomEffects, compare
import linearmodels as plm

# 构造/加载面板数据（entity + time 双索引）
np.random.seed(42)
n = 50  # 个体数
T = 10  # 时间期数
entity = np.repeat(np.arange(n), T)
time = np.tile(np.arange(2010, 2020), n)
x = np.random.normal(0, 1, n*T)  # 自变量
alpha = np.random.normal(5, 2, n)  # 个体固定效应
y = alpha[entity] + 2*x + np.random.normal(0, 1, n*T)  # 因变量

df = pd.DataFrame({
    "entity": entity,
    "time": time,
    "y": y,
    "x": x
}).set_index(["entity", "time"])  # 必须设双索引

# 拟合固定效应模型（双向：entity + time）
fe_model = PanelOLS.from_formula(
    "y ~ x + EntityEffects + TimeEffects",  # EntityEffects/TimeEffects 固定效应
    data=df
)
fe_results = fe_model.fit(cov_type="clustered", cluster_entity=True)  # 聚类标准误
print(fe_results)

Hausman 检验（FE vs RE 选择）

# 随机效应模型
re_model = RandomEffects.from_formula("y ~ x", data=df)
re_results = re_model.fit()

# Hausman检验（p<0.05选FE）
hausman_test = compare({
    "FE": fe_results,
    "RE": re_results
}, "hausman")
print(hausman_test)

个体固定效应

import pandas as pd
import statsmodels.api as sm
from statsmodels.formula.api import ols

# 数据同上（不设双索引）
df = df.reset_index()

# 个体固定效应：加入 entity 虚拟变量
model = ols("y ~ x + C(entity)", data=df)
results = model.fit()
print(results.summary())

双向固定效应

import pandas as pd
import statsmodels.api as sm
from statsmodels.formula.api import ols

# 构造示例面板数据
df = pd.DataFrame({
    'id':      [1,1,1, 2,2,2, 3,3,3],   # 个体id
    'year':    [2018,2019,2020]*3,      # 年份
    'y':       [1,2,3, 2,3,4, 3,4,5],
    'x1':      [0.1,0.2,0.3, 0.2,0.3,0.4, 0.3,0.4,0.5]
})

# 双向固定效应：C(id) 个体FE + C(year) 时间FE
model = ols('y ~ x1 + C(id) + C(year)', data=df)
# model = ols('y ~ x1 + C(id) + C(year) - 1', data=df) # 有时会去掉截距，不加常数项
res = model.fit()
print(res.summary())

更推荐：python运行 from linearmodels import PanelOLS

因为它 直接支持 EntityEffects + TimeEffects, 自动输出组内 R2, 支持聚类标准误, 支持 Hausman 检验.

时间序列分析

移动平均 (MA)：消除短期波动

自回归模型

  • AR (p) 自回归模型：当前值 = 历史值线性组合 + 噪声
  • MA (q) 移动平均模型：当前值 = 历史误差项线性组合 + 均值
  • ARMA(p,q)：AR + MA，适用于平稳序列
  • ARIMA(p,d,q)：ARMA 扩展，通过 d 阶差分 将非平稳转为平稳，是最经典通用模型
  • from statsmodels.tsa.arima.model import ARIMA

典型应用领域

  • 金融：股价预测、风险价值（VaR）、高频交易、宏观经济指标（GDP / 通胀）
  • 零售电商：销量预测、库存优化、节假日需求分析
  • 能源：电力负荷预测、风电 / 光伏出力预测
  • 工业物联网：设备故障预警、传感器信号监控、寿命预测
  • 交通：车流量预测、拥堵分析、网约车调度
  • 医疗：心率 / 血压监测、流行病趋势预测
  • 气象：温度、降水、空气质量预报
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.tsa.arima.model import ARIMA
from statsmodels.tsa.stattools import adfuller
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
import warnings
warnings.filterwarnings("ignore")

# ====================== 1. 构造模拟时间序列数据 ======================
np.random.seed(42)
n = 120
t = np.arange(n)
# 趋势 + 季节 + 噪声
trend = 0.1 * t
season = 3 * np.sin(2 * np.pi * t / 12)
noise = np.random.normal(0, 1, n)
ts = trend + season + noise

# 转成时间索引
dates = pd.date_range(start='2014-01-01', periods=n, freq='M')
ts = pd.Series(ts, index=dates)

# ====================== 2. 查看数据 ======================
plt.figure(figsize=(12,4))
plt.plot(ts)
plt.title('时间序列')
plt.show()

# ====================== 3. ADF平稳性检验 ======================
def adf_test(series):
    res = adfuller(series)
    print('ADF统计量:', res[0])
    print('p值:', res[1])
    print('临界值:', res[4])
    if res[1] < 0.05:
        print('结论：平稳')
    else:
        print('结论：非平稳')

print("=== ADF平稳性检验 ===")
adf_test(ts)

# ====================== 4. 画ACF/PACF定阶 ======================
fig, (ax1, ax2) = plt.subplots(1,2,figsize=(12,4))
plot_acf(ts, lags=20, ax=ax1)
plot_pacf(ts, lags=20, ax=ax2)
plt.show()

# ====================== 5. 拟合 ARIMA 模型 ======================
# 这里手动选 p=2, d=1, q=2
model = ARIMA(ts, order=(2,1,2))
result = model.fit()
print(result.summary())

# ====================== 6. 预测未来12个月 ======================
forecast_steps = 12
forecast = result.get_forecast(steps=forecast_steps)
forecast_mean = forecast.predicted_mean
forecast_ci = forecast.conf_int()

# ====================== 7. 绘图展示 ======================
plt.figure(figsize=(12,5))
plt.plot(ts, label='历史')
plt.plot(forecast_mean.index, forecast_mean, label='预测', color='red')
plt.fill_between(forecast_ci.index,
                 forecast_ci.iloc[:,0],
                 forecast_ci.iloc[:,1],
                 color='pink', alpha=0.3)
plt.legend()
plt.title('ARIMA 预测')
plt.show()

面板数据分析

用固定效应模型

面板数据（同时包含横截面与时间序列维度的数据）

核心要点

  • 优先选固定效应（控制不可观测异质性最稳妥）
  • 动态面板必用 GMM（解决滞后因变量内生性）
  • 因果识别：DID、合成控制、工具变量是主流
  • 必须报告：模型选择依据、稳健标准误、内生性处理

分析流程

  • 数据准备：长格式、缺失值、描述统计
  • 平稳性检验（单位根）→ 避免伪回归
  • 模型选择：F / LM / Hausman 检验
  • 回归估计（FE/RE/GMM 等）
  • 诊断：异方差、自相关、截面相关
  • 稳健标准误（聚类标准误）
  • 结果解读与机制分析

经济学

经济学规律

规律和指标。

微观经济学
研究个体决策：消费者、厂商、市场价格与资源配置。

  • 需求法则（Law of Demand）：其他条件不变时，价格↑ → 需求量↓；价格↓ → 需求量↑。 - 表现为向右下方倾斜的需求曲线。
  • 供给法则（Law of Supply）：其他条件不变时，价格↑ → 供给量↑；价格↓ → 供给量↓。 - 表现为向右上方倾斜的供给曲线。
  • 供求均衡法则：价格高于均衡价：供给 > 需求 → 过剩 → 价格下跌 - 价格低于均衡价：需求 > 供给 → 短缺 → 价格上涨 - 市场最终趋向均衡价格与均衡数量。
  • 边际效用递减法则：连续消费同一种商品，每新增一单位带来的满足感越来越小。 - 解释：为什么多吃一碗饭没第一碗香。
  • 边际产量/边际收益递减法则：短期内，其他要素不变，持续增加某一投入，边际产量最终会下降。 - 是成本曲线向上弯曲的根本原因。
  • 等边际原则（消费者最优）：消费者效用最大时 $\frac{MU_1}{P_1}=\frac{MU_2}{P_2}=\dots$ - 花在每种商品上的最后一块钱带来的边际效用相等。
  • 利润最大化原则（厂商最优）： 任何市场结构下，最优产量满足： $MR = MC$ 边际收益 = 边际成本。
  • 成本法则：短期：平均总成本先降后升（U型），源于边际收益递减。 - 长期：规模报酬先递增、后不变、最终递减。
  • 市场效率法则：完全竞争市场长期均衡： $P = MC = AC_{\text{min}}$，资源配置最有效率。 - 垄断、外部性、公共品、信息不对称 → 市场失灵。
  • 比较优势法则（贸易基础）：即使一国所有产品生产率都更低，仍可通过机会成本更低的产品进行贸易，双方都获益。

宏观经济学
研究整体经济：GDP、通胀、失业、增长、货币、财政政策。

  • 国民收入核算恒等式： $Y = C + I + G + NX$ ，总产出 = 消费 + 投资 + 政府购买 + 净出口。
  • 边际消费倾向递减法则：收入越高，新增收入中用于消费的比例越小，储蓄比例越高。 - 是凯恩斯总需求不足理论的基础。
  • 乘数效应：投资、政府支出增加 → 国民收入成倍增加。 $\text{乘数} = \frac{1}{1-MPC}$
  • 菲利普斯曲线法则（短期）： 短期内：通胀↑ → 失业↓；通胀↓ → 失业↑。 - 长期：失业回到自然失业率，通胀与失业无关。
  • 货币数量论 $MV = PY$ - 货币供应量增速长期决定通胀率。 - 货币中性：长期货币只影响价格，不影响实际产出。
  • 奥肯定律：产出缺口与失业率负相关 实际GDP每低于潜在GDP 2%左右，失业率大约上升1%。
  • 总需求-总供给法则：总需求AD右移 产出↑、价格↑ → 经济扩张、通胀上升，总供给AS左移 产出↓、价格↑ → 滞胀 - 长期AS垂直：需求只影响价格，不影响产出。
  • 挤出效应：政府大幅增加支出 → 利率上升 → 私人投资减少。 - 财政政策效果被部分抵消。
  • 购买力平价（汇率长期法则）： 长期汇率由两国物价水平之比决定，同一商品在各国价格应趋同。
  • 经济增长核心法则：长期增长主要来自： 1) 资本积累 2) 劳动投入 3) 技术进步（最重要） - 无技术进步时，增长最终收敛到稳态。

最精简总结

  • 微观：价格调节供求、边际递减、最优决策、市场效率与失灵。
  • 宏观：总需求决定短期产出、货币决定通胀、失业与增长交替、长期靠技术。

经济学指标

微观经济学指标：

  • 价格与数量类 - 价格 P：商品或要素的单位价格 - 需求量 Qd：消费者愿意且能购买的数量 - 供给量 Qs：厂商愿意且能出售的数量 - 均衡价格 Pe、均衡数量 Qe：市场出清时的价格与数量
  • 弹性指标 - 需求价格弹性 (E_d)：需求量对价格变动的反应 - 供给价格弹性 (E_s)：供给量对价格变动的反应 - 需求收入弹性 (E_I)：判断正常品/劣等品/奢侈品 - 需求交叉价格弹性 (E_{xy})：判断替代品/互补品
  • 消费者行为 - 总效用 TU、边际效用 MU - 消费者剩余 CS：消费者愿意支付与实际支付的差额 - 无差异曲线、边际替代率 MRS - 预算线、最优消费组合
  • 生产者行为 - 总产量 TP、平均产量 AP、边际产量 MP - 边际技术替代率 MRTS - 总成本 TC、固定成本 FC、可变成本 VC - 平均成本 AC、AVC、AFC - 边际成本 MC - 总收益 TR、平均收益 AR、边际收益 MR - 利润 π = TR − TC
  • 市场结构与效率 - 生产者剩余 PS - 总剩余 TS = CS + PS（衡量社会福利） - 无谓损失 DWL：垄断/税收/管制造成的效率损失 - 勒纳指数：衡量垄断势力 (L=(P-MC)/P) - 集中度指标：CR4、HHI（市场垄断程度）
  • 要素市场 - 工资 w、利率 r、地租 R - 边际产品价值 VMP、边际收益产品 MRP - 要素需求、要素供给
  • 福利与一般均衡 - 帕累托最优（无改进空间） - 基尼系数（微观+宏观通用，衡量收入分配） - 外部性、公共物品、搭便车相关效率指标

宏观经济学指标：

  • 国民经济总量： - GDP（国内生产总值）一定时期内一国境内生产的全部最终产品和服务的市场价值。 - GNP（国民生产总值）一国国民在国内外生产的最终产品和服务价值。 - GNI（国民总收入）居民、企业、政府在国内外获得的全部收入总和。
  • 增长与产出： - 经济增长率：GDP 同比 / 环比增速 - 潜在 GDP：资源充分利用时的最大产出 - 产出缺口：实际 GDP − 潜在 GDP
  • 就业与失业： - 失业率：失业人口 / 劳动力人口 - 劳动参与率：劳动力 / 劳动年龄人口 - 非农就业、就业人数
  • 物价与通胀： - CPI（居民消费价格指数） - PPI（工业生产者价格指数） - GDP 平减指数：衡量整体通胀 - 核心通胀：剔除食品、能源后的通胀
  • 货币与金融： - M0、M1、M2（货币供应量） - 利率：政策利率、市场利率 - 汇率：本币对外币比价 - 外汇储备
  • 财政与债务： - 财政收入 / 支出 - 财政赤字 / 盈余 - 政府债务率：政府债务 / GDP - 税收收入
  • 国际收支： - 经常账户：贸易 + 服务 + 收益 + 转移 - 资本与金融账户 - 贸易顺差 / 逆差 - 国际收支总差额
  • 居民与消费： - 人均可支配收入 - 社会消费品零售总额 - 储蓄率、消费率
  • 投资与产业： - 固定资产投资 - 工业增加值 - PMI（采购经理人指数）：景气先行指标

参考

工具： statsmodels、linearmodels
论文，期刊：

经典论文

  • 《》
  • 《》

教材

  • 《基本无害的计量经济学》

中文期刊

  1. 数字技术创新与中国企业高质量发展——来自企业数字专利的证据 黄勃 等 2023 经济研究
  2. 数智化如何影响双循环参与度与收入差距——基于省级—行业层面… 张云和柏培文 2023 管理世界
  3. 人工智能如何提升企业生产效率？——基于劳动力技能结构调整的… 姚加权 等 2024 管理世界
  4. 人工智能技术创新与中国劳动力市场势力：从“二元增长分化”到劳…

回归分析

import numpy as np
import statsmodels.api as sm
import statsmodels.formula.api as smf
# Load data
dat = sm.datasets.get_rdataset("Guerry", "HistData").data
# Fit regression model (using the natural log of one of the regressors)
results = smf.ols('Lottery ~ Literacy + np.log(Pop1831)', data=dat).fit()
# Inspect the results
print(results.summary())

返回结果

  • Dep. Variable: Lottery 因变量：你要预测的东西 → 彩票参与率
  • R-squared: 模型解释力： Adj. R-squared: 调整后解释力：
  • F-statistic 整体模型是否显著 Prob (F) 模型整体非常显著（p<0.05）
  • No. Observations 一共个样本
  • 逐个变量的 P 值，越小越显著。模型是否健康 所有 p值都 >0.05

面板数据

from linearmodels import PanelOLS
mod = PanelOLS.from_formula('invest ~ value + capital + EntityEffects', data)
res = mod.fit(cov_type='clustered', cluster_entity=True)

工具变量法 Instrumental Variable Models IV regression models can be similarly specified.

import numpy as np
from linearmodels.iv import IV2SLS
from linearmodels.datasets import mroz
data = mroz.load()
mod = IV2SLS.from_formula('np.log(wage) ~ 1 + exper + exper ** 2 + [educ ~ motheduc + fatheduc]', data)

实际输出：

1. OLS（Guerry 数据）

OLS Regression Results                            
==============================================================================
Dep. Variable:                Lottery   R-squared:                       0.348
Model:                            OLS   Adj. R-squared:                  0.333
Method:                 Least Squares   F-statistic:                     22.20
Date:                Sun, 17 May 2026   Prob (F-statistic):           1.90e-08
No. Observations:                  86   AIC:                             765.6
===================================================================================
                      coef    std err          t      P>|t|
-----------------------------------------------------------------------------------
Intercept         246.4341     35.233      6.995      0.000
Literacy           -0.4889      0.128     -3.832      0.000
np.log(Pop1831)   -31.3114      5.977     -5.239      0.000

三个变量全部显著（p < 0.001），R² = 0.348。

2. PanelOLS（Munnell 基础设施数据）

                          PanelOLS Estimation Summary                           
================================================================================
Dep. Variable:                    GSP   R-squared:                        0.8406
Estimator:                   PanelOLS   R-squared (Within):              0.8809
No. Observations:                 816   F-statistic (robust):             19.959
Entities:                          48   P-value                           0.0000
Time periods:                      17   
======================================================================
             Parameter  Std. Err.     T-stat    P-value
Intercept   2.787e+04     8279.2     3.3665     0.0008
PC             0.8306     0.0846     9.8170     0.0000
UNEMP         -1193.3     246.71    -4.8367     0.0000

双向固定效应（Entity + Time），PC 和 UNEMP 显著，P_CAP/HWY/WATER/UTIL 不显著。

3. IV-2SLS（Mroz 工资数据）

                          IV-2SLS Estimation Summary                          
==============================================================================
Dep. Variable:           np.log(wage)   R-squared:                      0.1357
Estimator:                    IV-2SLS   Adj. R-squared:                 0.1296
No. Observations:                 428   F-statistic:                    24.653
Cov. Estimator:            unadjusted
====================================================================
                  Parameter  Std. Err.     T-stat    P-value
Intercept            0.0481     0.3985     0.1207     0.9039
exper                0.0442     0.0134     3.3038     0.0010
np.square(exper)    -0.0009     0.0004    -2.2485     0.0245
educ                 0.0614     0.0313     1.9622     0.0497
====================================================================
Endogenous: educ
Instruments: motheduc, fatheduc

用父母教育作为工具变量估计教育回报率，educ 在 5% 水平显著。

为什么要用 IV：教育（educ）与能力等不可观测因素相关，OLS 估计的”教育回报率”有偏。父母教育可以作为工具变量——影响个人教育但影响无法直接影响工资（通过能力以外的途径）。

相关工具

初步数据的整理通常用 Excel，有时候包括 Power Query。

Stata

Stata 因果推断核心代码（最常用方法）
我给你整理了实证研究中最主流、最常用的因果推断方法，全部是可直接运行、注释完整的 Stata 代码，包含：

  • 双重差分 DID（最常用）
  • 工具变量 IV
  • 匹配法 PSM
  • 断点回归 RDD
  • 固定效应模型 FE

所有代码都用标准公开数据格式，你直接替换变量名就能用。

你告诉我3个信息，我可以直接给你定制可运行代码

  1. 你的数据是截面/面板？
  2. 你要做的因果问题（比如：政策效果、某变量影响）
  3. 变量名：被解释变量、核心自变量、控制变量

我可以直接按你的数据写完整可复现的实证代码+结果解读。

双重差分法 DID（政策评估首选）
适用场景：政策实施前后、处理组 vs 对照组

* 1. 读取数据（自行替换你的数据）
use "your_data.dta", clear

* 2. 生成核心变量
* treat = 1 处理组，0 对照组
* time = 1 政策后，0 政策前
gen did = treat * time  // 交互项（核心系数）

* 3. 基准DID回归
reg y treat time did, robust  // 稳健标准误

* 4. 加入控制变量 + 个体/时间固定效应（推荐）
xtset id year  // 面板数据设定
xtreg y did x1 x2 i.year, fe robust  // 双向固定效应DID

* 5. 平行趋势检验（事件研究法）
gen period = year - policy_year  // 相对时间
eventstudy y, period(period) treat(treat) controls(x1 x2) fe robust

工具变量法 IV（解决内生性）
适用场景：存在遗漏变量、双向因果、测量误差

* 1. 两阶段最小二乘法 2SLS
ivregress 2sls y x1 x2 (endog_var = instrument1), robust

* 2. 检验工具变量有效性（必须做）
estat firststage  // 一阶段F值（>10才算强工具变量）
estat overid  // 过度识别检验（多个工具变量时用）

倾向得分匹配 PSM（观测数据因果识别）
适用场景：非随机分配的处理效应

* 1. 估计倾向得分
logit treat x1 x2 x3  // 用协变量预测处理状态
predict ps, pr  // 生成倾向得分

* 2. 最近邻匹配（1:1）
psmatch2 treat, pscore(ps) neighbor(1) ate robust

* 3. 匹配平衡性检验（必须做）
pstest x1 x2 x3, both graph

断点回归 RDD（准自然实验最强）
适用场景：以某个分数/阈值为界产生处理

* 1. 清晰断点回归
rdrobust y running_var, c(60)  // c=阈值，如60分及格

* 2. 画图检验断点
rdplot y running_var, c(60)

* 3. 带宽选择（最优带宽）
rdbwselect y running_var, c(60)

固定效应模型（面板数据基础）

* 个体固定效应
xtreg y x1 x2, fe robust

* 双向固定效应（个体+时间）
xtreg y x1 x2 i.year, fe robust

* 随机效应
xtreg y x1 x2, re robust

* Hausman检验（选FE还是RE）
hausman fe re

合成控制法 Synthetic Control（单案例政策）

* 安装包
ssc install synth

* 运行合成控制
synth y x1 x2, trunit(10) trperiod(2015)

R

核心概念与R生态
因果推断的核心是估计干预（处理）对结果的真实效应（ATE/ATT），解决观测数据中的混杂偏倚；R有成熟包生态，核心包如下：

  • 倾向得分/匹配/加权：MatchIt、WeightIt、CBPS
  • 双重差分（DID）：did、DRDID、fixest
  • 工具变量（IV）：ivreg、estimatr、DoubleML
  • 合成控制：Synth、gsynth、tidysynth
  • 因果图/可识别性：dagitty、causaleffect
  • 双重稳健/机器学习因果：DoubleML、grf

关键步骤总结

  1. 明确因果问题：定义处理(T)、结果(Y)、混杂/工具变量
  2. 画DAG，验证效应可识别、确定调整集
  3. 选方法：观测截面→PSM/IPTW/DR；面板→DID；内生→IV
  4. 平衡检验（匹配/加权后）、稳健性检验
  5. 估计ATE/ATT，报告标准误与置信区间

经典方法+完整R代码（基于lalonde数据集）

关键指标解读（直接看结果）

  1. ATT：仅看「实际接受干预人群」的平均效应
  2. ATE：全人群 假设干预的平均效应
  3. PSM/IPTW：看 treat 系数 = 因果效应值
  4. DID：显著为正 / 负 代表政策显著有效
  5. IV：规避双向因果、遗漏变量偏倚

R语言 因果推断

倾向得分匹配（PSM，估计ATT）
最常用，平衡处理/对照组协变量，消除混杂

# 安装并加载包
install.packages(c("MatchIt", "dplyr", "tableone"))
library(MatchIt)
library(dplyr)
library(tableone)

# 数据：lalonde（职业培训对收入的影响，treat=1为处理组，re78为结果）
data("lalonde", package = "MatchIt")
head(lalonde)

# 1) 构建倾向得分模型+最近邻1:1匹配
ps_match <- matchit(
  treat ~ age + educ + black + hisp + married + nodegr + re74 + re75,
  data = lalonde,
  method = "nearest",  # 最近邻
  ratio = 1,          # 1:1匹配
  caliper = 0.2        # 卡尺，控制匹配质量
)

# 2) 检验协变量平衡（关键！标准化差异<0.1为平衡）
summary(ps_match, standardize = TRUE)
# 可视化平衡
plot(ps_match, type = "jitter")

# 3) 提取匹配数据，估计ATT
matched_data <- match.data(ps_match)
att_model <- lm(re78 ~ treat, data = matched_data)
summary(att_model)
# ATT：培训使1978年收入平均增加XX元

倾向得分匹配 PSM（最常用，观测数据）

  • 最常用，平衡处理 / 对照组协变量，消除混杂
  • 使用 R 内置经典因果数据集 lalonde
  • 看 treat 系数 = 因果效应值
  • 估计 ATT
    ```r

    数据：treat=1 参与培训(处理组)，re78=1978年收入(结局)

    data(lalonde)
    head(lalonde)

1:1最近邻匹配 + 卡尺约束

psm_fit <- matchit(
treat ~ age + educ + black + hisp + married + re74 + re75,
data = lalonde,
method = “nearest”,
ratio = 1,
caliper = 0.2
)

协变量平衡检验（标准化均值差<0.1 平衡合格）

summary(psm_fit, standardize = TRUE)

提取匹配后数据

match_data <- match.data(psm_fit)

估计ATT 平均处理效应

psm_reg <- lm(re78 ~ treat, data = match_data)
summary(psm_reg)


---

逆概率加权（IPTW，估计ATE/ATT）
用倾向得分加权，无需删样本，适合大样本
```r
library(WeightIt)

# 估计倾向得分权重
ipw <- weightit(
  treat ~ age + educ + black + hisp + married + nodegr + re74 + re75,
  data = lalonde,
  method = "ps",  # 倾向得分加权
  estimand = "ATE" # 目标：平均处理效应ATE
)

# 检查权重平衡
summary(ipw)

# 加权回归估计ATE
ate_model <- lm(re78 ~ treat, data = lalonde, weights = ipw$weights)
summary(ate_model)

IPTW 逆概率加权（无需删样本）

  1. PSM/IPTW：看 treat 系数 = 因果效应值
    ```r

    计算倾向得分权重

    ipw_fit <- weightit(
    treat ~ age + educ + black + hisp + married + re74 + re75,
    data = lalonde,
    method = “ps”,
    estimand = “ATE” # 估计总体平均处理效应
    )

权重平衡诊断

summary(ipw_fit)

加权回归 估计因果效应

ipw_reg <- lm(re78 ~ treat, data = lalonde, weights = ipw_fit$weights)
summary(ipw_reg)


---

双重差分（DID，面板/重复截面，估计政策效应）
适合**干预前后+处理/对照组**的面板数据，控制时间趋势与个体固定效应
```r
install.packages(c("did", "fixest"))
library(did)
library(fixest)

# 构造模拟面板数据（或用真实面板）
# 核心：att_gt（Callaway & Sant'Anna，交错DID，处理时间不同）
# 示例：处理组(treat=1)在time=2后接受干预，outcome为结果
# 数据需含：id(个体)、time(时期)、treat(是否处理)、group(处理时间)

# 估计DID-ATT
did_model <- att_gt(
  yname = "outcome",
  tname = "time",
  idname = "id",
  gname = "group",  # 处理开始时间
  xformla = ~ age + educ, # 协变量
  data = panel_data,
  est_method = "dr" # 双重稳健，更稳健
)

# 汇总结果
summary(did_model)
# 动态效应+平均效应
aggdid <- aggte(did_model, type = "dynamic")
summary(aggdid)

双重差分 DID（政策评估/面板数据）

  • 模拟交错政策冲击面板数据，经典政策因果评估
    1. DID：显著为正/负 代表政策显著有效
      ```r

      构造面板数据

      n_id <- 200
      time_len <- 6
      panel_data <- expand.grid(id = 1:n_id, time = 1:time_len) %>%
      mutate(
      # 分组：g=3 代表第3期开始实施政策
      group = sample(c(Inf,3,4), nrow(.), replace = TRUE),
      treat = ifelse(time >= group & group != Inf, 1, 0),
      x1 = rnorm(nrow(.)),
      # 生成结果：政策效应+时间趋势+个体效应
      Y = 2treat + 0.3time + 0.5*x1 + rnorm(nrow(.))
      )

Callaway & Sant’Anna 交错DID

did_fit <- att_gt(
yname = “Y”,
tname = “time”,
idname = “id”,
gname = “group”,
xformla = ~ x1,
data = panel_data,
control_group = “never_treated”
)

总体平均政策效应

did_avg <- aggte(did_fit, type = “simple”)
summary(did_avg)

动态效应（事件研究）

did_dyn <- aggte(did_fit, type = “dynamic”)
summary(did_dyn)


---

 


工具变量（IV，解决内生性，如遗漏变量/双向因果）
用**仅影响处理、不直接影响结果**的工具变量，2SLS估计
```r
install.packages("estimatr")
library(estimatr)

# 示例：教育(treat)对收入(outcome)的影响，用"是否出生在大学附近"作IV
# 公式：outcome ~ treat | IV + 协变量
iv_model <- iv_robust(
  re78 ~ treat + age + educ | instrument + age + educ, # 第二阶段|第一阶段
  data = lalonde_iv,
  se_type = "HC1" # 稳健标准误
)

summary(iv_model)
# 第一阶段F值>10，说明IV不是弱工具

工具变量 IV（解决内生性/双向因果）

  • 模拟数据： - T = 受教育程度(内生处理) - Y = 工资收入 - IV = 出生地是否靠近大学（工具变量）
    1. IV：规避双向因果、遗漏变量偏倚
      ```r

      模拟数据

      n <- 1000
      iv_data <- tibble(
      IV = rbinom(n,1,0.5),
      C = rnorm(n),
      T = 0.6IV + 0.4C + rnorm(n), # 第一阶段
      Y = 1.2T + 0.5C + rnorm(n) # 结果方程
      )

2SLS 工具变量回归

iv_fit <- iv_robust(
Y ~ T + C | IV + C, # 结果方程 | 工具变量+外生协变量
data = iv_data,
se_type = “HC2”
)
summary(iv_fit)


---

双重稳健估计（DR，PSM+回归，更稳健）
结合加权/匹配与回归，任一模型正确即一致估计
```r
library(DRDID)
# 适用于面板DID的双重稳健
dr_did <- drdid(
  yname = "re78",
  tname = "time",
  idname = "id",
  dname = "treat",
  xformla = ~ age + educ + re74,
  data = panel_data,
  boot = TRUE, boot_reps = 200
)
summary(dr_did)

双重稳健 DR 估计（抗混杂、模型容错）

# 双重稳健：倾向得分加权 + 结果回归
dr_fit <- weightit(
  treat ~ age + educ + re74 + re75,
  data = lalonde,
  method = "ps"
)

# 加权+回归调整，双重稳健
dr_reg <- lm(re78 ~ treat + age + educ + re74 + re75, 
             data = lalonde, 
             weights = dr_fit$weights)
summary(dr_reg)

因果图与可识别性（dagitty）
先画DAG，识别混杂/工具变量/后门路径，确保效应可识别

install.packages("dagitty")
library(dagitty)

# 构建DAG：处理(T)→结果(Y)，混杂(C)影响T和Y
dag <- dagitty("dag {
  C -> T
  C -> Y
  T -> Y
}")

# 识别后门调整集（需要控制的混杂）
adjustmentSets(dag, exposure = "T", outcome = "Y")
# 可视化DAG
plot(dag)

因果DAG图（先识别混杂，因果分析前提）

# 定义因果有向无环图
dag <- dagitty('dag {
  混杂变量C -> 处理T
  混杂变量C -> 结果Y
  处理T -> 结果Y
}')

# 寻找需要控制的混杂变量（后门准则）
adjustmentSets(dag, exposure = "处理T", outcome = "结果Y")

# 绘图
plot(dag)

回归分析

固定效应模型

Python

DID

常用 statsmodels、linearmodels（面板数据）或 CausalPy 库。

  • https://causalpy.readthedocs.io/en/stable/

Models can also be specified using the formula interface.

from linearmodels import PanelOLS
mod = PanelOLS.from_formula('invest ~ value + capital + EntityEffects', data)
res = mod.fit(cov_type='clustered', cluster_entity=True)

The formula interface for PanelOLS supports the special values EntityEffects and TimeEffects which add entity (fixed) and time effects, respectively.

Formula support comes from the formulaic package which is a replacement for patsy.

调用 OLS 有矩阵写法（类似 MatLab）和公式写法（类似R），个人偏向后者。

可以复现一篇完整的论文来得到具体的代码，如果作者有提供最好，是在不行也可以借助 ChatGPT 的帮忙。

数据预处理
对齐论文样本筛选、变量定义、1%缩尾、面板索引设定：

def data_preprocess(df):
    """
    数据预处理（严格对齐论文3.1研究设计）
    1. 剔除金融/IT/科研行业、ST/*ST、缺失值
    2. 连续变量1%缩尾
    3. 设置面板索引（公司-年份）
    """
    # 样本筛选
    exclude_ind = ['金融', '信息传输、软件和信息技术服务业', '科学研究和技术服务业']
    df = df[~df['industry'].isin(exclude_ind)]
    df = df[~df['st_status'].isin(['ST', '*ST'])]
    
    # 剔除核心变量缺失值
    core_vars = ['TFP', 'Lnwords', 'Size', 'Age', 'Leverage', 'Growth', 'BoardSize', 'Dual', 'Top1', 'Lnallpats']
    df = df.dropna(subset=core_vars)
    
    # 连续变量1%缩尾
    cont_vars = ['TFP', 'Lnwords', 'Lnwords_MD&A', 'Lnpatents', 'Routine', 'Non_routine', 'Tobinq'] + core_vars[2:]
    for var in cont_vars:
        df[var] = df[var].clip(lower=df[var].quantile(0.01), upper=df[var].quantile(0.99))
    
    # 面板数据索引
    df = df.set_index(['code', 'year'])
    return df

测算全要素生产率TFP

# 补充：OP法测算全要素生产率TFP（论文3.3核心方法）
def calculate_tfp(df):
    """OP法计算企业TFP（Y=销售额，L=员工数，K=固定资产，I=投资）"""
    df['lnY'] = np.log(df['sales'] + 1)
    df['lnL'] = np.log(df['employee'] + 1)
    df['lnK'] = np.log(df['fixed_assets'] + 1)
    df['lnI'] = np.log(df['investment'] + 1)
    
    # 面板回归估算TFP
    model = PanelOLS.from_formula('lnY ~ lnL + lnK + EntityEffects + TimeEffects', data=df)
    res = model.fit(cov_type='clustered', cluster_entity=True)
    df['TFP'] = res.resids
    return df

描述性统计
复现论文核心变量描述、年度AI披露/专利统计：

def descriptive_stats(df):
    """核心变量描述性统计+年度AI应用情况"""
    # 核心变量统计
    key_vars = ['TFP', 'Lnwords', 'Lnwords_MD&A', 'Lnpatents', 'Routine', 'Non_routine', 'Tobinq', 'Size', 'Leverage']
    desc = df[key_vars].describe().T.round(4)
    
    # 年度AI披露/专利企业数
    ai_disclosure = (df['Lnwords'] > 0).groupby('year').sum()
    ai_patent = (df['Lnpatents'] > 0).groupby('year').sum()
    
    print("="*50)
    print("核心变量描述性统计（论文表2 Panel B）")
    print("="*50)
    print(desc)
    
    print("\n" + "="*50)
    print("年度AI披露与专利申请（论文表2 Panel A）")
    print("="*50)
    print(pd.DataFrame({'AI披露企业数': ai_disclosure, 'AI专利企业数': ai_patent}))
    return desc

基准回归，固定效应模型
OLS+行业聚类标准误+固定效应，对齐论文基准模型：

def baseline_regression(df):
    """
    基准回归：TFP = α + βAI + γControls + year+ind+pro + ε
    标准误：行业层面聚类
    """
    # 被解释变量/解释变量/控制变量
    y = df['TFP']
    X = df[['Lnwords', 'Size', 'Age', 'Leverage', 'Growth', 'BoardSize', 'Dual', 'Top1', 'Lnallpats']]
    X = sm.add_constant(X)
    
    # 固定效应（年份、行业、省份）
    fe = pd.get_dummies(df[['year', 'industry', 'province']], drop_first=True)
    X = pd.concat([X, fe], axis=1)
    
    # OLS回归+行业聚类
    model = sm.OLS(y, X).fit(cov_type='cluster', cov_kwds={'groups': df['industry']})
    
    print("\n" + "="*50)
    print("基准回归结果（论文表3）")
    print("="*50)
    print(model.summary().tables[1])
    return model

用于面板数据，固定效应，如实体，时间，各有一个基数。

有时候可以先回归再因果，把因果作为稳健性检验。
倾向得分匹配 PSM

def psm_test(df):
    """1:1近邻PSM匹配+回归（缓解自选择偏差）"""
    # 定义处理组：应用AI企业（Lnwords>0）
    df['treat'] = (df['Lnwords'] > 0).astype(int)
    covars = ['Size', 'Age', 'Leverage', 'Growth', 'BoardSize', 'Dual', 'Top1', 'Lnallpats']
    
    # 倾向得分Logit估计
    logit_X = sm.add_constant(df[covars])
    logit = sm.Logit(df['treat'], logit_X).fit(disp=0)
    df['pscore'] = logit.predict(logit_X)
    
    # 1:1近邻匹配
    treat = df[df['treat']==1]
    control = df[df['treat']==0]
    nn = NearestNeighbors(n_neighbors=1).fit(control[['pscore']])
    _, idx = nn.kneighbors(treat[['pscore']])
    matched_df = pd.concat([treat, control.iloc[idx.flatten()]])
    
    # 匹配后回归
    y_m = matched_df['TFP']
    X_m = sm.add_constant(matched_df[['Lnwords'] + covars])
    fe_m = pd.get_dummies(matched_df[['year', 'industry', 'province']], drop_first=True)
    X_m = pd.concat([X_m, fe_m], axis=1)
    psm_model = sm.OLS(y_m, X_m).fit(cov_type='cluster', cov_kwds={'groups': matched_df['industry']})
    
    print("\n" + "="*50)
    print("PSM匹配后回归（论文表3列3）")
    print("="*50)
    print(psm_model.summary().tables[1])
    return matched_df, psm_model

先匹配出处理组和对照组从而平衡样本，然后再回归。
以缓解直接回归时的偏差。

2SLS工具变量法

def iv_2sls_test(df):
    """
    工具变量：SEA_PORT×LnAI（清朝通商口岸×滞后1期全球AI专利）
    2SLS面板回归（论文表4）
    """
    # 第一阶段：Lnwords ~ IV + Controls + FE
    # 第二阶段：TFP ~ 拟合Lnwords + Controls + FE
    iv_model = IV2SLS.from_formula(
        'TFP ~ 1 + Size + Age + Leverage + Growth + BoardSize + Dual + Top1 + Lnallpats + '
        'C(year) + C(industry) + C(province) ~ [Lnwords ~ IV]',
        data=df.reset_index()
    ).fit(cov_type='clustered', cluster_entity=True)
    
    print("\n" + "="*50)
    print("2SLS工具变量回归（论文表4）")
    print("="*50)
    print(iv_model)
    return iv_model

工具变量->结果
双重差分法 DID

TODO

中介效应检验
劳动力技能结构中介：AI→替代低技能劳动力/互补高技能劳动力→TFP

def mediation_test(df):
    """逐步回归法检验中介效应（论文表7）"""
    controls = 'Size + Age + Leverage + Growth + BoardSize + Dual + Top1 + Lnallpats'
    fe = 'C(year) + C(industry) + C(province)'
    
    # 1. 总效应：TFP ~ AI
    model1 = PanelOLS.from_formula(f'TFP ~ Lnwords + {controls} + {fe}', data=df).fit()
    
    # 2. 中介1：常规低技能劳动力（Routine）
    model2_r = PanelOLS.from_formula(f'Routine ~ Lnwords + {controls} + {fe}', data=df).fit()
    model3_r = PanelOLS.from_formula(f'TFP ~ Lnwords + Routine + {controls} + {fe}', data=df).fit()
    
    # 3. 中介2：非常规高技能劳动力（Non_routine）
    model2_n = PanelOLS.from_formula(f'Non_routine ~ Lnwords + {controls} + {fe}', data=df).fit()
    model3_n = PanelOLS.from_formula(f'TFP ~ Lnwords + Non_routine + {controls} + {fe}', data=df).fit()
    
    print("\n" + "="*50)
    print("中介效应检验（论文表7）")
    print("="*50)
    print(f"总效应（AI→TFP）：{model1.params['Lnwords'].round(4)}")
    print(f"AI→低技能劳动力：{model2_r.params['Lnwords'].round(4)}")
    print(f"低技能中介效应：{model3_r.params['Routine'].round(4)}")
    print(f"AI→高技能劳动力：{model2_n.params['Lnwords'].round(4)}")
    print(f"高技能中介效应：{model3_n.params['Non_routine'].round(4)}")
    return model1, model2_r, model3_r, model2_n, model3_n

比较系数的大小，以及看是否显著。

异质性分析

def heterogeneity_test(df):
    """分组回归：产权/行业/地区异质性（论文表8-10）"""
    # 1. 产权异质性
    soe = df[df['soe']==1]
    non_soe = df[df['soe']==0]
    m_soe = PanelOLS.from_formula(f'TFP ~ Lnwords + Size + Age + Leverage + Growth + BoardSize + Dual + Top1 + Lnallpats + C(year)+C(industry)+C(province)', data=soe).fit()
    m_non = PanelOLS.from_formula(f'TFP ~ Lnwords + Size + Age + Leverage + Growth + BoardSize + Dual + Top1 + Lnallpats + C(year)+C(industry)+C(province)', data=non_soe).fit()
    
    # 2. 行业异质性：高技术vs低技术
    high_tech = df[df['high_tech']==1]
    m_high = PanelOLS.from_formula(f'TFP ~ Lnwords + Controls + FE', data=high_tech).fit()
    
    print("\n" + "="*50)
    print("异质性分析（论文表8-10）")
    print("="*50)
    print(f"国企AI系数：{m_soe.params['Lnwords'].round(4)}")
    print(f"非国企AI系数：{m_non.params['Lnwords'].round(4)}")
    print(f"高技术行业AI系数：{m_high.params['Lnwords'].round(4)}")
    print(f"低技术行业AI系数：{m_low.params['Lnwords'].round(4)}")

比较数据分组后，每组系数大小。

企业价值回归

def firm_value_test(df):
    """AI对企业价值（Tobinq）的影响（论文表11）"""
    model = PanelOLS.from_formula(
        'Tobinq ~ Lnwords + Size + Age + Leverage + Growth + BoardSize + Dual + Top1 + Lnallpats + '
        'C(year) + C(industry) + C(province)',
        data=df
    ).fit(cov_type='clustered', cluster_entity=True)
    
    print("\n" + "="*50)
    print("人工智能与企业价值（论文表11）")
    print("="*50)
    print(model.summary())
    return model

C (…) = Categorical 的缩写 → 告诉模型：这一列是「分类变量」，要自动生成「虚拟变量（哑变量）」纳入回归。 简单说： 把数字 / 文本的分类变量，转成回归能用的 0/1 虚拟变量，不用你手动生成。

casualpy

注意结果的可视化解读，比如说差异体现在什么地方。

import causalpy as cp  # 导入因果推断库 CausalPy

# ==================== 1. Difference-in-Differences ====================
# DiD 通过比较处理组和对照组在干预前后的变化差异来估计因果效应。
# 核心思想：用对照组的变化作为处理组的反事实。
# group:post_treatment 交互项系数即为因果效应。
# 运行结果：Causal impact = 0.50, 94% CI [0.41, 0.60]
print("=" * 60)
print("1. Difference-in-Differences (DiD)")
print("Compares treatment vs control over pre/post periods")
print("=" * 60)

df = cp.load_data("did")                          # 加载内置 DiD 数据集
result = cp.DifferenceInDifferences(                # 创建 DiD 模型
	df,                                             # 输入数据
	formula="y ~ 1 + group + t + group:post_treatment",  # 回归公式，交互项 group:post_treatment 的系数 = 因果效应
	time_variable_name="t",                         # 时间变量名（0=干预前, 1=干预后）
	group_variable_name="group",                    # 分组变量名（0=对照组, 1=处理组）
)
fig, ax = result.plot()                             # 绘制 DiD 结果图（含反事实和因果效应）
fig.savefig("1_did.png")                            # 保存为图片
print(result.summary())                             # 打印模型摘要（包含系数和因果效应）

# ==================== 2. Regression Discontinuity ====================
# RD 在驱动变量的阈值处附近比较结果变量的跳跃，用于估计因果效应。
# 适用于"临界值决定处理"的场景（如分数线、年龄门槛）。
# treated 是是否超过阈值（x > 0.5）的指示变量。
# 运行结果：Discontinuity at threshold = 0.19, 94% CI [-0.07, 0.45]
print("\n" + "=" * 60)
print("2. Regression Discontinuity (RD)")
print("Causal effect at a cutoff threshold on running variable")
print("=" * 60)

df = cp.load_data("rd")                             # 加载内置 RD 数据集
result = cp.RegressionDiscontinuity(                 # 创建 RD 模型
	df,                                             # 输入数据
	formula="y ~ 1 + x + treated",                  # y 对 x 和 treated 回归，treated 系数 = 断点处的跳跃
	running_variable_name="x",                      # 驱动变量名
	treatment_threshold=0.5,                        # 处理阈值，x > 0.5 时 treated = 1
)
fig, ax = result.plot()                             # 绘制 RD 结果图（含断点两侧的拟合线）
fig.savefig("2_rd.png")                             # 保存为图片
print(result.summary())                             # 打印模型摘要（含断点效应估计）

# ==================== 3. Interrupted Time Series ====================
# ITS 分析单个时间序列在已知干预时刻前后的变化。
# 通过比较干预前后的斜率和截距变化来估计影响。
# t 为时间索引，treatment_time=100 处为干预点。
# 运行结果：Intercept = 33, t (斜率) = 0.24
print("\n" + "=" * 60)
print("3. Interrupted Time Series (ITS)")
print("Impact of intervention at known time in a single time series")
print("=" * 60)

df = cp.load_data("its")                            # 加载内置 ITS 数据集
result = cp.InterruptedTimeSeries(                   # 创建 ITS 模型
	df,                                             # 输入数据
	formula="y ~ 1 + t",                            # y 对时间 t 回归，比较干预前后的参数变化
	time_variable_name="t",                         # 时间变量名
	treatment_time=100,                             # 干预发生的时间点
)
fig, ax = result.plot()                             # 绘制 ITS 结果图（含干预前后的拟合线 + 反事实）
fig.savefig("3_its.png")                            # 保存为图片
result.print_coefficients()                         # 打印模型系数

# ==================== 4. Instrumental Variable ====================
# IV 通过工具变量解决内生性问题。这里用 nearcollege (幼年是否靠近大学)
# 作为 college (是否上大学) 的工具变量，估计大学教育对工资的因果效应。
# 第一阶段：college ~ nearcollege_bin；第二阶段：wage ~ college。
# beta_t 是 college 的因果效应系数。
# 运行结果：Causal effect = 0.1, 95% CI [0.1, 0.2]
print("\n" + "=" * 60)
print("4. Instrumental Variable (IV)")
print("Effect of college on wage, instrumented by nearcollege")
print("=" * 60)

df = cp.load_data("schoolReturns")                  # 加载内置 IV 数据集（Card, 1995）
df["nearcollege_bin"] = (df["nearcollege"] == "yes").astype(int)  # 将字符串工具变量转为 0/1
df["college"] = (df["education"] >= 16).astype(int) # 创建二值处理变量：是否上大学（education >= 16）

result = cp.InstrumentalVariable(                    # 创建 IV 模型（两阶段最小二乘）
	instruments_data=df,                            # 工具变量所在数据
	data=df,                                        # 结果变量和自变量所在数据
	instruments_formula="college ~ nearcollege_bin", # 第一阶段：用工具变量预测 college
	formula="wage ~ 1 + college",                   # 第二阶段：用 college 预测 wage（college 的系数 = 因果效应）
)
beta_t = result.idata.posterior["beta_t"]            # 从后验分布中提取因果效应系数
effect = float(beta_t.mean())                        # 后验均值
ci = [float(beta_t.quantile(q)) for q in [0.025, 0.975]]  # 95% 可信区间
print(f"Causal effect (college on wage): {effect:.1f}")
print(f"95% CI: [{ci[0]:.1f}, {ci[1]:.1f}]")

# ==================== 5. Inverse Propensity Weighting ====================
# IPW 通过倾向得分加权来模拟随机化，减少选择偏差。
# 这里估计 NSW 就业培训项目对 1978 年收入的因果效应。
# 用协变量预测接受处理的概率（倾向得分），再以 IPW 方案加权。
# 运行结果：Raw diff = -$635, Propensity score mean = 0.301
print("\n" + "=" * 60)
print("5. Inverse Propensity Weighting (IPW)")
print("ATE of NSW job training on 1978 earnings")
print("=" * 60)

df = cp.load_data("lalonde")                        # 加载 LaLonde 数据集（经典因果推断数据）
df["treat"] = df["treat"].astype(int)               # 确保处理变量为整型

result = cp.InversePropensityWeighting(               # 创建 IPW 模型
	df,                                             # 输入数据
	formula="treat ~ 1 + age + educ + married + nodegree + re74 + re75",  # 用协变量估计倾向得分（处理概率）
	outcome_variable="re78",                        # 结果变量：1978 年收入
	weighting_scheme="ipw",                         # 加权方案：逆概率加权
)
treated = df[df["treat"] == 1]                      # 筛选处理组
control = df[df["treat"] == 0]                      # 筛选对照组
raw_diff = treated["re78"].mean() - control["re78"].mean()  # 未加权的原始均值差（有偏估计）
print(f"Raw difference in means: ${raw_diff:.0f}")  # 打印未调整的简单差异
print(f"Propensity score mean:   {result.t.mean():.3f}")   # 倾向得分均值
print(f"Average propensity score model coefficient: {float(result.idata.posterior['b'].mean()):.3f}")  # 倾向得分模型系数的后验均值

print("\n" + "=" * 60)
print("All results saved as PNG plots")
print("=" * 60)

SPSS

SPSS 通常用于调查问卷分析。开源的有 PSPP。

附录

//@狸角兽:我觉得主要靠的是“作者对于真理的追求”和“愿意为真相付出的代价”。若是要求不高的话，把相关关系当因果关系的文章，大把；更广泛的，跑个IV跑个DID 因果关系就搭好了，也能上顶刊。要求高时，哪怕穷尽观察与测量，也不过是逼近真相的序章。 所谓推断，根本上是倔强：我们终究无法见证所 …展开
@包特_ExpEcon 我觉得“因果推断”这个概念可能是有一些问题的，因为更加好的识别因果一般都是靠观察/测量，而不是推断～ ​​​

//@狸角兽:我觉得主要靠的是“作者对于真理的追求”和“愿意为真相付出的代价”。若是要求不高的话，把相关关系当因果关系的文章，大把；更广泛的，跑个IV跑个DID 因果关系就搭好了，也能上顶刊。要求高时，哪怕穷尽观察与测量，也不过是逼近真相的序章。所谓推断，根本上是倔强：我们终究无法见证所有反事实，但仍选择不轻易妥协。//@包特_ExpEcon:我觉得“因果推断”这个概念可能是有一些问题的，因为更加好的识别因果一般都是靠观察/测量，而不是推断～ ​​​

补充1

读物

  • 《Mostly Harmless Econometrics: An Empiricist’s Companion》

参考论文

中文论文

  • 管理世界 《人工智能如何提升企业生产效率》

稳健性检验

目前有点罗列资料了，还缺乏重点和故事（每一步的动机）。

本科常用流程：描述性统计 → 相关性分析 → 基准回归 → 稳健性检验 → 异质性检验
硕士常用流程：在本科基础上，增加：内生性处理、机制分析

这是本科论文最标准的“基础款”流程，核心目标是验证「X和Y之间有没有稳定的相关关系」。

  1. 描述性统计：先给读者看数据的基本情况（均值、标准差、极值等），让大家对变量有个直观印象。
  2. 相关性分析：用相关系数矩阵，初步判断变量之间的方向和强度，也能顺便排查多重共线性问题。
  3. 基准回归：论文的核心，用OLS等方法估计核心解释变量对被解释变量的影响系数和显著性，回答“X对Y有没有影响、影响多大”。
  4. 稳健性检验：换模型、换变量定义、换子样本，看核心结果会不会变，证明你的结论不是偶然得到的。
  5. 异质性检验：把样本分组（比如分地区、分企业规模），看X对Y的影响在不同群体里有没有差异，让论文更丰富。

硕士论文的核心要求，是从「证明相关关系」升级到「论证因果关系」，所以在本科流程的基础上，多了两个关键步骤：

  1. 内生性处理：解决“X和Y的关系可能是假的”问题（比如双向因果、遗漏变量偏差）。常用方法就是工具变量、双重差分、倾向得分匹配等，目的是让你的估计结果更接近真实的因果效应。
  2. 机制分析：不仅要证明“X会影响Y”，还要说清楚「X是通过什么路径/渠道影响Y的」。比如“数字经济（X）提升企业绩效（Y），是通过技术创新（M）这个中介实现的”，让论文的逻辑链条更完整，研究深度大大提升。

一句话总结：本科流程是“把故事讲圆，证明相关关系”；硕士流程是“把因果坐实，讲清影响路径”。

博士可以向理论深度、方法创新和问题边界全面升级，核心是从 “验证别人的理论” 走向 “构建 / 修正理论”，并做出具有原创性的学术贡献。

政策评估

政策评估的通用步骤（计量实操）

  1. 明确研究问题：定义政策（处理）、结果变量、处理组 / 对照组、政策时点。
  2. 数据准备：面板数据优先（个体 × 时间），整理协变量（人口、经济、制度等）。
  3. 基准回归：OLS / 固定效应，得到初步政策系数。
  4. 内生性处理：选择 DID/RDD/PSM/IV 等准实验方法，构建反事实。
  5. 假设检验：平行趋势（DID）、连续性（RDD）、平衡性（PSM）、工具有效性（IV）。
  6. 稳健性检验：替换变量、调整样本、安慰剂检验（随机分配政策）、动态效应分析。
  7. 异质性分析：分地区、分群体、分时间检验政策效果差异。
  8. 机制分析：检验政策通过何种渠道（如投资、消费、就业）影响结果。

前沿拓展（近年主流）

  • 多期 DID（ staggered DID）：解决政策分批实施的估计偏误（如 TWFE 偏误），采用事件研究 + 分组固定效应或Callaway-Sant’Anna（CS）估计。
  • 机器学习 + 政策评估：用随机森林、LASSO 选择协变量，提升 PSM/RDD 的匹配与拟合精度。
  • 空间计量 + 政策评估：考虑政策空间溢出效应（如邻区政策影响），采用空间 DID、空间滞后模型。
  • 中介效应 / 调节效应：分析政策的作用机制（中介）与边界条件（调节，如制度环境、市场化程度）。

总结

政策评估计量方法的核心是 “找反事实、控内生性、证因果”。DID是通用首选，RDD是高可信度首选，PSM用于匹配校正，SCM用于单一主体，IV用于内生性难题，事件研究用于动态分析。实际研究中，常以**一种主方法 + 多种稳健性检验 + 异质性 / 机制分析 ** 组合，确保结论可靠。

需要我帮你整理一份政策评估计量方法的实操模板（含 DID/RDD/PSM 的 Stata/R 代码框架、假设检验命令与稳健性检验步骤）吗？

量化交易

股票三大经典技术指标实战代码（Python）
我给你写一套可直接运行、开箱即用的Python代码，包含：

  1. 双均线策略（金叉/死叉）
  2. MACD + 顶底背离
  3. 布林带突破策略

使用最常用的库：pandas + ta（技术指标库），数据用模拟K线，你替换成真实数据即可。

先安装依赖

pip install pandas ta numpy

完整代码（直接复制运行）

import pandas as pd
import ta
import numpy as np

# ====================== 1. 生成模拟K线数据（可替换为你的真实数据） ======================
def generate_sample_data():
    # 生成200根K线：时间、开盘价、收盘价、最高价、最低价
    date_range = pd.date_range(start="2025-01-01", periods=200, freq="D")
    close_prices = np.cumsum(np.random.randn(200)) + 100  # 随机收盘价
    high = close_prices + np.random.uniform(0, 2, 200)
    low = close_prices - np.random.uniform(0, 2, 200)
    open_ = close_prices + np.random.uniform(-1, 1, 200)

    df = pd.DataFrame({
        "date": date_range,
        "open": open_,
        "high": high,
        "low": low,
        "close": close_prices
    })
    return df

# ====================== 2. 双均线指标（金叉 / 死叉） ======================
def add_ma_signal(df):
    # 计算5日、20日均线（常用短长周期）
    df["ma_short"] = ta.trend.sma_indicator(df["close"], window=5)
    df["ma_long"] = ta.trend.sma_indicator(df["close"], window=20)

    # 金叉：短均线上穿长均线 → 买入信号
    df["ma_gold_cross"] = (df["ma_short"] > df["ma_long"]) & (df["ma_short"].shift(1) <= df["ma_long"].shift(1))
    # 死叉：短均线下穿长均线 → 卖出信号
    df["ma_dead_cross"] = (df["ma_short"] < df["ma_long"]) & (df["ma_short"].shift(1) >= df["ma_long"].shift(1))
    return df

# ====================== 3. MACD + 顶背离 / 底背离 ======================
def add_macd_and_divergence(df):
    # 计算MACD指标
    macd = ta.trend.MACD(df["close"])
    df["macd"] = macd.macd()
    df["macd_signal"] = macd.macd_signal()
    df["macd_diff"] = macd.macd_diff()  # MACD柱

    # === 顶背离：价格创新高，MACD不创新高 → 看跌
    df["price_high"] = df["close"].rolling(window=10).max()
    df["macd_high"] = df["macd"].rolling(window=10).max()
    df["bearish_divergence"] = (df["close"] == df["price_high"]) & (df["macd"] < df["macd_high"].shift(1))

    # === 底背离：价格创新低，MACD不创新低 → 看涨
    df["price_low"] = df["close"].rolling(window=10).min()
    df["macd_low"] = df["macd"].rolling(window=10).min()
    df["bullish_divergence"] = (df["close"] == df["price_low"]) & (df["macd"] > df["macd_low"].shift(1))
    return df

# ====================== 4. 布林带突破策略 ======================
def add_bollinger_signal(df):
    # 计算布林带（默认20日均值，2倍标准差）
    bollinger = ta.volatility.BollingerBands(df["close"])
    df["boll_mid"] = bollinger.bollinger_mavg()
    df["boll_up"] = bollinger.bollinger_hband()
    df["boll_low"] = bollinger.bollinger_lband()

    # 突破上轨 → 强势上涨信号
    df["boll_break_up"] = df["close"] > df["boll_up"]
    # 跌破下轨 → 弱势下跌信号
    df["boll_break_down"] = df["close"] < df["boll_low"]
    return df

# ====================== 主函数：一键计算所有指标 ======================
if __name__ == "__main__":
    # 1. 获取数据
    df = generate_sample_data()

    # 2. 计算所有指标
    df = add_ma_signal(df)          # 双均线
    df = add_macd_and_divergence(df)# MACD + 背离
    df = add_bollinger_signal(df)   # 布林带

    # 3. 输出结果（只看关键列）
    result = df[[
        "date", "close",
        "ma_short", "ma_long", "ma_gold_cross", "ma_dead_cross",
        "macd", "bullish_divergence", "bearish_divergence",
        "boll_up", "boll_low", "boll_break_up", "boll_break_down"
    ]].tail(30)  # 显示最后30行

    print("📊 技术指标计算完成（最后30根K线）")
    print(result)

代码说明（一看就懂）

  1. 双均线
    • ma_short：5日均线（快线）
    • ma_long：20日均线（慢线）
    • 金叉：ma_gold_cross = True → 买入信号
    • 死叉：ma_dead_cross = True → 卖出信号
  2. MACD + 背离
    • macd：MACD线
    • macd_signal：信号线
    • 底背离：bullish_divergence = True → 看涨
    • 顶背离：bearish_divergence = True → 看跌
  3. 布林带突破
    • boll_up：上轨
    • boll_low：下轨
    • 突破上轨：boll_break_up = True → 多头强势
    • 跌破下轨：boll_break_down = True → 空头强势

如何替换成你的真实股票数据？
把 generate_sample_data() 替换成：

df = pd.read_csv("你的K线文件.csv")
# 确保包含字段：open, high, low, close

总结
这套代码包含了你要的双均线、MACD背离、布林带突破三大策略，直接运行就能出信号，适合量化交易、策略回测、看盘分析。
还可以：

  • 改成可视化K线图
  • 做成自动交易信号提示
  • 适配A股/期货/加密货币

演示用数据集

R自带的？