## 1. 为什么金融实证研究要从STATA转向Python？ 我刚开始做公司金融实证研究那会儿，身边几乎所有人都在用STATA。导师给的代码是STATA的，师兄师姐分享的经验是STATA的，甚至很多顶级期刊论文的复现资料也是STATA的。STATA确实是个好工具，点几下菜单、敲几行命令，结果就出来了，特别适合快速出结果。但用久了，我总觉得哪里不对劲——我好像变成了一个“命令操作员”，而不是一个真正理解模型的研究者。 后来因为一个偶然的项目需要处理超大规模的面板数据，STATA有点力不从心，我才硬着头皮开始用Python。这一用，就再也回不去了。我最大的感受是，Python让研究的“黑箱”变得透明了。在STATA里，你输入 `reg y x1 x2 i.province`，它直接给你吐出一张漂亮的回归表。但中间的数据是怎么处理的？虚拟变量是怎么生成的？标准误是怎么计算的？如果你不去深究，这些细节就像被封装在一个魔法盒子里。 而Python，尤其是 `statsmodels` 这个库，逼着你去理解每一个步骤。你得自己用 `pandas` 清洗数据，用 `sm.categorical` 或者 `pd.get_dummies` 去创建虚拟变量，亲手用 `sm.add_constant` 添加截距项，最后调用 `sm.OLS().fit()` 得到结果。这个过程听起来麻烦，但实际上，它让你对模型构建的逻辑有了前所未有的掌控感。你的代码就是你的研究笔记，每一步都清晰可循，可重复性极高。 更重要的是，Python的生态太强大了。数据预处理有 `pandas` 和 `numpy`，复杂到让人头疼的缺失值处理、变量转换，几行代码就能搞定。做完回归，你想画一张更精美的趋势图来展示核心变量关系？`matplotlib` 和 `seaborn` 让你随心所欲。你想把多个回归结果输出成论文要求的格式？`summary_col` 函数可以帮你生成媲美STATA `outreg2` 的表格。这些工作流上的无缝衔接，是STATA难以比拟的。所以，如果你已经熟悉了STATA的基本计量思想，那么用Python的 `statsmodels` 上手会非常快，它就像是STATA的一个“开源透明版”，既能让你沿用熟悉的计量逻辑，又能打开新世界的大门。 ## 2. 快速上手：statsmodels核心模块与OLS初体验 咱们先别被 `statsmodels` 官网那一长串的模块列表吓到。对于从STATA迁移过来的金融实证研究，你最开始只需要盯紧两个核心：**`statsmodels.api`** (通常简写为 `sm`) 和 **`statsmodels.formula.api`** (通常简写为 `smf`)。`smf` 允许你使用类似R语言的公式语法，比如 `y ~ x1 + x2`，对于习惯STATA命令式思维的同学可能更亲切。但我个人更推荐从 `sm` 的数组风格入手，因为它能让你更清楚地看到数据是如何被组织进模型的，这和STATA底层处理数据的逻辑更接近。 安装就是一行命令的事：`pip install statsmodels`。如果遇到网络问题，加上清华镜像源 `-i https://pypi.tuna.tsinghua.edu.cn/simple` 会快很多。 接下来，我们直接看一个最小二乘（OLS）的例子，这是所有实证研究的基石。假设我们想研究公司规模（`size`）和资产负债率（`lev`）对企业研发投入（`rd`）的影响，同时控制年份效应。 ```python import pandas as pd import numpy as np import statsmodels.api as sm # 1. 模拟一份简单的公司年度面板数据 np.random.seed(2023) n_firms = 100 n_years = 5 firm_ids = np.repeat(range(n_firms), n_years) years = np.tile(range(2018, 2023), n_firms) # 生成变量 size = np.random.lognormal(mean=10, sigma=1.0, size=n_firms*n_years) lev = np.random.uniform(0.1, 0.8, size=n_firms*n_years) # rd 受到 size 和 lev 的影响，并加入年份固定效应和随机噪声 year_dummy = (years - 2018) * 0.05 # 模拟一个随时间增长的趋势 rd = 0.3 * size - 0.5 * lev + year_dummy + np.random.normal(0, 0.5, size=n_firms*n_years) df = pd.DataFrame({'firm_id': firm_ids, 'year': years, 'rd': rd, 'size': size, 'lev': lev}) print(df.head()) ``` 数据有了，现在进行OLS回归。在STATA里，我们可能直接 `reg rd size lev i.year`。在Python里，我们需要显式地处理虚拟变量： ```python # 2. 准备因变量和自变量 y = df['rd'] # 因变量：研发投入 # 自变量：连续变量 size 和 lev X_continuous = df[['size', 'lev']] # 为年份生成虚拟变量（注意：避免虚拟变量陷阱，我们丢弃第一年作为基准组） year_dummies = pd.get_dummies(df['year'], prefix='year', drop_first=True) # 将所有自变量水平拼接（concat） X = pd.concat([X_continuous, year_dummies], axis=1) # 3. 添加常数项（截距）。在STATA的reg命令中，截距是默认包含的。 X = sm.add_constant(X) # 4. 建立并拟合OLS模型 model = sm.OLS(y, X) results = model.fit() # 5. 查看完整的回归结果摘要 print(results.summary()) ``` 运行这段代码，你会得到一张非常熟悉的回归结果表，里面包含了系数估计、标准误、t值、p值、R方等所有关键统计量。这张表和STATA输出的核心部分几乎一模一样。第一次在Python里看到这个，我当时的感受是：“稳了！” 这意味着，我可以用Python完整地复现STATA的基准分析。`results.summary()` 这个函数是 `statsmodels` 的神器之一，它把计量结果格式化得清清楚楚，直接就能贴到论文初稿里。 ## 3. 攻克核心难点：虚拟变量的生成与固定效应控制 虚拟变量，或者说固定效应控制，是面板数据实证研究的灵魂。在STATA里，用 `i.year` 或 `i.province` 就能轻松搞定，后台自动帮你处理好了虚拟变量的生成和共线性问题。迁移到Python，我们需要自己动手，但这恰恰是理解固定效应本质的好机会。 **方法一：使用 `pandas.get_dummies`** 这是最直观、最常用的方法，特别适合处理单个或多个分类变量。 ```python # 假设我们的df中还有‘province’省份列 # 生成省份虚拟变量，并丢弃第一个省份作为基准组（避免完全共线性） province_dummies = pd.get_dummies(df['province'], prefix='prov', drop_first=True) # 如果你要同时控制年份和省份，就分别生成再合并 year_dummies = pd.get_dummies(df['year'], prefix='year', drop_first=True) # 合并到自变量矩阵X中 X = pd.concat([df[['size', 'lev']], year_dummies, province_dummies], axis=1) X = sm.add_constant(X) ``` **方法二：使用 `statsmodels` 自带的 `categorical` 函数** 这个方法更接近计量软件底层的处理方式，能更好地处理复杂的面板数据结构。在原始文章的论文复现中，作者就用了这个方法。 ```python # 使用 sm.categorical 生成虚拟变量。drop=True 表示丢弃第一组以避免共线性。 # 它返回的是一个数组（或稀疏矩阵），适合直接用于模型拟合。 dummy_year = sm.categorical(df['year'].values, drop=True) dummy_province = sm.categorical(df['province'].values, drop=True) # 将连续变量和虚拟变量堆叠起来 import numpy as np X_cont = df[['size', 'lev']].values X = np.column_stack((X_cont, dummy_year, dummy_province)) # 别忘了添加常数项 X = sm.add_constant(X) # 这时X就是一个纯粹的NumPy数组，可以直接用于sm.OLS model = sm.OLS(df['rd'].values, X) results = model.fit() ``` **这里有个非常重要的坑我踩过：** 当你使用 `pd.get_dummies` 时，如果某些分类在某个样本子集中不存在（比如某省份只有部分年份有数据），在合并或后续样本筛选时，可能会导致虚拟变量列数不一致，从而报错。而 `sm.categorical` 在处理面板数据时，能更稳定地保持虚拟变量结构与样本的对应关系。我个人的经验是，对于简单的截面数据或平衡面板，用 `pandas` 更方便；对于复杂的非平衡面板数据，用 `sm.categorical` 更稳妥。 **什么是“固定效应”？** 其实在上述操作中，我们通过加入年份、省份的虚拟变量，就是在控制“年份固定效应”和“省份固定效应”。这相当于在回归中允许不同年份、不同省份有自己的截距项，从而吸收掉那些不随时间（或个体）变化，但会影响因变量的不可观测因素。这和STATA中 `xtreg, fe` 的思想在本质上是相通的，只不过 `xtreg, fe` 是通过组内离差变换来实现的，而加入虚拟变量是另一种等价形式（LSDV法）。 ## 4. 实战进阶：复现一篇论文的基准回归 光说不练假把式。我们现在就模仿原始文章，尝试复现一篇经典论文的核心部分。假设我们想研究“数字化转型”（`digital`）对企业全要素生产率（`tfp`）的影响，并控制企业规模（`size`）、年龄（`age`）、年份固定效应和行业固定效应。 我们模拟一份更贴近真实研究的数据： ```python # 模拟数据 n = 5000 # 5000个观测值 df = pd.DataFrame({ 'firm_id': np.random.choice(range(500), n), 'year': np.random.choice([2019, 2020, 2021, 2022], n), 'industry': np.random.choice(['制造业', '服务业', '信息技术业', '金融业'], n), 'size': np.random.lognormal(12, 1.2, n), 'age': np.random.randint(1, 30, n), 'digital': np.random.normal(0, 1, n), # 数字化转型指数 }) # 模拟TFP，使其与digital、size、age相关，并加入行业和年份的固定效应 industry_effect = {'制造业': -0.2, '服务业': 0.1, '信息技术业': 0.5, '金融业': 0.3} year_effect = {2019: 0, 2020: -0.1, 2021: 0.2, 2022: 0.4} df['tfp'] = (0.15 * df['digital'] + 0.3 * np.log(df['size']) - 0.01 * df['age'] + df['industry'].map(industry_effect) + df['year'].map(year_effect) + np.random.normal(0, 0.5, n)) ``` 现在，我们构建回归模型。目标：`tfp` 对 `digital`、`ln(size)`、`age` 回归，同时控制年份和行业固定效应。 ```python # 1. 数据预处理：对size取对数，这是金融实证中处理规模变量的常见做法 df['ln_size'] = np.log(df['size']) # 2. 生成固定效应虚拟变量 year_dummies = pd.get_dummies(df['year'], prefix='year', drop_first=True) industry_dummies = pd.get_dummies(df['industry'], prefix='ind', drop_first=True) # 3. 组装自变量矩阵 X = pd.concat([df[['digital', 'ln_size', 'age']], year_dummies, industry_dummies], axis=1) X = sm.add_constant(X) y = df['tfp'] # 4. 执行OLS回归 model = sm.OLS(y, X, missing='drop') # missing='drop'会自动删除含有缺失值的行 results = model.fit() # 5. 输出漂亮的回归结果 print(results.summary()) ``` 查看 `results.summary()`，你会重点关注 `digital` 变量的系数、标准误和p值。如果系数显著为正，就初步支持了“数字化转型提升生产率”的假设。这和你在STATA里跑回归、看结果的过程完全一致。 但Python的强大不止于此。我们可以轻松地扩展分析，比如进行**分组回归**： ```python # 按企业规模中位数分组，检验数字化转型的效果是否存在异质性 median_size = df['size'].median() df_large = df[df['size'] > median_size] df_small = df[df['size'] <= median_size] # 定义一个函数来重复上面的回归流程 def run_ols(data): X = pd.concat([data[['digital', 'ln_size', 'age']], pd.get_dummies(data['year'], prefix='year', drop_first=True), pd.get_dummies(data['industry'], prefix='ind', drop_first=True)], axis=1) X = sm.add_constant(X) y = data['tfp'] model = sm.OLS(y, X, missing='drop') return model.fit() res_large = run_ols(df_large) res_small = run_ols(df_small) print("=== 大企业组回归结果 ===") print(res_large.summary().tables[1]) # 只打印系数表 print("\n=== 小企业组回归结果 ===") print(res_small.summary().tables[1]) ``` 通过对比两组中 `digital` 的系数，你可以分析其影响是否因企业规模而异。在STATA里，你可能需要用 `if` 条件或 `by` 命令，但在Python里，这种灵活的数据切片和函数封装让分析流程更加清晰和自动化。 ## 5. 结果输出与可视化：让实证结果“会说话” 做完回归，得到了一堆系数和p值，接下来就是如何呈现它们。STATA有 `outreg2` 或 `esttab` 这样的神器，可以把多个回归结果输出成LaTeX或Word兼容的格式。Python同样可以做到，而且更加灵活。 **技巧一：使用 `summary_col` 输出专业回归表格** `statsmodels.iolib.summary2` 模块中的 `summary_col` 函数，是模仿 `outreg2` 的利器。 ```python from statsmodels.iolib.summary2 import summary_col # 假设我们跑了三个模型：只放digital；加入控制变量；加入所有固定效应 X1 = sm.add_constant(df[['digital']]) X2 = sm.add_constant(df[['digital', 'ln_size', 'age']]) # X3 是之前包含所有虚拟变量的完整模型 model1 = sm.OLS(y, X1, missing='drop').fit() model2 = sm.OLS(y, X2, missing='drop').fit() model3 = results # 上面已经拟合过的完整模型 # 将三个模型的结果合并输出 table = summary_col([model1, model2, model3], model_names=['模型(1)', '模型(2)', '模型(3)'], stars=True, # 显示显著性星标 float_format='%0.4f', # 系数保留4位小数 info_dict={'N': lambda x: f"{int(x.nobs):,}", # 样本量 'R2': lambda x: f"{x.rsquared:.4f}"}) # R方 print(table) ``` 这段代码会生成一个规整的表格，横向是三个不同的模型设定，纵向是各个变量及其统计量，下方还有样本量和R方的信息。你可以直接把打印出来的文本复制到文本编辑器里，稍作调整就能用于论文。 **技巧二：关键结果的可视化** 统计表格虽然精确，但不够直观。我们可以用图表来辅助说明核心发现。 ```python import matplotlib.pyplot as plt import seaborn as sns # 示例1：绘制核心解释变量与被解释变量的散点图与拟合线 plt.figure(figsize=(8,5)) sns.regplot(x=df['digital'], y=df['tfp'], scatter_kws={'alpha':0.3, 's':10}, line_kws={'color':'red'}) plt.xlabel('数字化转型指数 (digital)') plt.ylabel('全要素生产率 (tfp)') plt.title('数字化转型与企业生产率的关系') plt.grid(True, linestyle='--', alpha=0.5) plt.show() # 示例2：比较分组回归的核心系数（假设我们上面已经得到了res_large和res_small） coef_data = pd.DataFrame({ '组别': ['大企业', '小企业'], '系数': [res_large.params['digital'], res_small.params['digital']], '标准误': [res_large.bse['digital'], res_small.bse['digital']] }) plt.figure(figsize=(6,4)) plt.errorbar(coef_data['组别'], coef_data['系数'], yerr=coef_data['标准误']*1.96, fmt='o', capsize=5, markersize=8) plt.axhline(y=0, color='grey', linestyle='--', linewidth=0.8) plt.ylabel('数字化转型的系数估计') plt.title('数字化转型对企业生产率的影响：按规模分组') plt.grid(axis='y', linestyle='--', alpha=0.5) plt.show() ``` 第一张图可以直观展示 `digital` 和 `tfp` 之间的正相关关系，第二张图则清晰地显示了大企业和小企业中 `digital` 系数估计值及其置信区间的差异。这些图表能极大地增强你论文结果部分的说服力。在STATA中做类似的图可能没那么方便，通常需要将结果导出再用其他软件绘制，而在Python中，从数据分析到结果呈现是一条流畅的流水线。 ## 6. 避坑指南：从STATA到Python的常见问题 迁移过程中肯定会遇到一些“水土不服”的问题，我这里总结几个最常见的坑和解决方案。 **坑1：标准误不一致** 在STATA中，默认的 `reg` 命令使用的是普通标准误。而在金融、经济学实证中，我们经常需要使用稳健标准误（Robust Standard Errors）或聚类标准误（Clustered Standard Errors）来处理异方差或组内相关性。在Python中，我们需要在拟合模型时显式指定。 ```python # 使用稳健标准误（HC3是一种常用的异方差稳健估计） results_robust = model.fit(cov_type='HC3') print(results_robust.summary()) # 此时表格中的标准误和t值就是稳健的 # 使用聚类标准误（例如，在公司层面聚类） # 注意：statsmodels的OLS本身不支持聚类标准误，需要使用`get_robustcov_results`或`linearmodels`库 # 这里展示使用`linearmodels`的简便方法（需安装：pip install linearmodels） from linearmodels.panel import PooledOLS from linearmodels import PanelOLS # 先将数据设置为面板格式 df_panel = df.set_index(['firm_id', 'year']) # 使用PooledOLS并指定聚类稳健协方差矩阵 model_pooled = PooledOLS(df_panel['tfp'], sm.add_constant(df_panel[['digital', 'ln_size', 'age']])) results_pooled_cluster = model_pooled.fit(cov_type='clustered', cluster_entity=True) print(results_pooled_cluster) ``` **坑2：虚拟变量陷阱与多重共线性** 当你手动生成虚拟变量时，如果不丢弃一个基准组，就会导致“虚拟变量陷阱”（完全共线性），模型会无法估计。`pd.get_dummies(..., drop_first=True)` 和 `sm.categorical(..., drop=True)` 中的 `drop` 参数就是为了解决这个问题。如果你在结果摘要中看到某个变量被自动丢弃，或者条件数（Cond. No.）极大，通常就是遇到了严重的多重共线性，需要检查你的虚拟变量设置。 **坑3：缺失值处理** STATA的回归命令通常会默认删除含有缺失值的样本行。在Python的 `pandas` 中，你需要格外小心。`sm.OLS` 默认不会自动删除缺失值，如果 `X` 或 `y` 里有 `NaN`，它会直接报错。因此，要么在数据预处理阶段用 `df.dropna(subset=[...])` 清理，要么在建模时使用 `sm.OLS(y, X, missing='drop')` 参数。 **坑4：内存与速度** 对于超大型数据集（比如百万级观测值），使用 `pandas.get_dummies` 生成大量虚拟变量可能会消耗巨大内存。这时可以考虑： 1. 使用 `sm.categorical` 并设置 `drop=True`，它默认返回一个内存效率更高的稀疏矩阵。 2. 对于真正的超大规模面板固定效应模型，可以考虑专门为面板数据设计的 `linearmodels` 库中的 `PanelOLS`，它在计算效率和功能上更强大。 迁移到Python不是一个一蹴而就的过程，开始时可能会觉得每一步都比STATA繁琐。但一旦你熟悉了这套流程，你会发现它带来的灵活性、透明性和可扩展性，会让你的研究效率和质量都上一个台阶。最关键的是，你对自己研究的掌控力更强了，每一个数字、每一个结果你都清楚地知道它是怎么来的。这种踏实感，是单纯点鼠标、敲命令无法给予的。