名师讲堂|使用 Python 测算地区产业专业化指标

由于借助 AI 工具学习编程已经变得非常容易了，因此之后的课程就不再默认进行视频讲解了，如果特别需要视频讲解也可以联系李老师预约讲解～讲义材料学习过程中遇到的问题也可以及时与李老师联系。

购买 RStata 名师讲堂会员即可参加该课程啦（之前的和未来的都可以参加）！

价格：2800/年 或者 4800/长期

购买会员可以从这里下单：https://rstata.duanshu.com/#/card/list/

名师讲堂会员权益：

1. 参加每个月 3～4 次的名师讲堂课程；
2. 参加平台上的其他 R 语言和 Stata 的课程；
3. 以会员折扣价购买我们分享的数据资料（10 元/份）；
4. 课程内外的提问解答服务（课程外的尽量帮忙解决）。

* 如果发票可添加小编微信 r_stata2 （RStata 李老师）开具。如需数据资料，购买后可添加小编微信免费领取数据折扣卡。

更多关于 RStata 会员的更多信息可添加微信号 r_stata2 咨询：

课程主页（点击文末的阅读原文即可跳转）：https://rstata.duanshu.com/#/brief/course/9dfc640436054f4cbba41308f7dc065d

今天给大家分享使用 Python 测算地区产业专业化指标的方法。该方法参考自杨本建、唐金汶（2022）《数字经济与区域产业布局》中的公式 (2)，其思想源于 Kalemli-Ozcan et al. (2003) 与 Du et al. (2022) 的 Krugman 式专业化指数。

附件中提供了该参考文献的 PDF 文件（数字经济与区域产业布局.pdf），感兴趣的小伙伴可以阅读原文。

本文使用的原始数据为工商企业注册信息（已在本项目内裁剪为测算所需的 9 个变量），以 2000–2005 年为例演示完整计算过程。与 R 版本不同，本文的 Python 代码全部运行在由 reticulate 创建和管理的独立虚拟环境（.venv）中，依赖（numpy / pandas / matplotlib）与系统 Python 完全隔离。

指标来源与计算过程

地区产业专业化指数（公式）

行业范围：制造业 C13–C43，剔除 C39

按国民经济行业分类（GB/T 4754）的制造业门类（行业门类 = "制造业"），取 2 位行业大类代码 C13–C43，并剔除 C39（计算机、通信和其他电子设备制造业）。剔除 C39 的依据是论文第 248 页附录 1 的说明（该行业受数字经济影响特殊，在相关研究中通常单独处理）。

产业规模的代理变量

论文以企业的注册资本作为行业规模的代理变量（主指标）；同时以企业数量作为稳健性口径。两者计算逻辑完全一致，仅在汇总时替换聚合字段。

存续企业的界定（进入与退出）

参考李磊等（2023）的做法，按"进入—退出"口径统计每年各城市各行业的存续企业：

• 进入：以企业的成立年份作为进入年份；
• 退出：当经营状态含"注销/吊销"时视为退出企业，退出年份取核准日期的年份；
• 存续：在年份 t 满足"成立年份 ≤ t，且（未退出 或 退出年份 > t）"的企业。

计算步骤概述

整体计算分为以下几个步骤：

1. 读取与清洗：读取各年工商注册数据，仅保留制造业、有效城市代码、必需变量，并界定进入/退出年份；
2. 构建存续面板：对每个目标年份 t，筛选存活企业，按"城市 × 行业"汇总规模，得到城市×行业×年份存续企业规模面板；
3. 计算专业化指数：对每一年，按公式计算各城市的 specspecspec 指数（注册资本口径 + 企业数量口径）；
4. 输出与可视化：导出结果 CSV，并绘制趋势、分布与 Top 城市等图表。

使用 reticulate 创建与管理 Python 虚拟环境

在 R 中通过 reticulate 包来调用 Python，最好的实践是为项目创建一个专属的 Python 虚拟环境，将所需依赖隔离到独立空间，避免与系统 Python（如 Anaconda）发生版本冲突。

重要说明（避免"已初始化"报错）：reticulate 在 R 会话中只能绑定一次 Python——一旦某个 {python} 代码块运行，Python 解释器就被锁定，之后再调用 use_virtualenv() 会报错：

ERROR: The requested version of Python cannot be used, as another version has already been initialized.

因此，虚拟环境的激活必须在所有 {python} 代码块之前完成。本文档的解决方案是在 setup chunk 中通过 Sys.setenv(RETICULATE_PYTHON = ...) 提前锁定 Python 路径，这是 reticulate 选取 Python 的最高优先级入口。

安装 reticulate（仅首次）

虚拟环境初始化原理（已在 setup chunk 中完成）

本文档的 setup chunk（隐藏运行）包含如下逻辑：

这样做的关键在于：knitr 在处理第一个 {python} chunk 时，reticulate 已经通过 RETICULATE_PYTHON 环境变量知道要使用 .venv，不会再去碰 Anaconda。

在虚拟环境中安装 Python 包（仅首次）

验证激活状态

查看已安装的包

虚拟环境管理常用命令

详细计算代码

下面按步骤完整展示 Python 代码（基于 pandas / numpy），每一段均可在 reticulate 管理的虚拟环境中直接运行。代码与同目录下的 01_测算地区产业专业化指标.py / 02_可视化.py 一致。

0. 路径与参数

首先设定工程目录、输出目录、数据目录，以及制造业行业代码与所需变量。

1. 读取并清洗单个年份文件

read_one_year() 负责把一年的 CSV 读入并清洗为"企业级"明细：

• 仅保留行业门类 == "制造业"且大类在 C13–C43 且非 C39 的记录；
• 仅保留市代码为 6 位数字行政区划码的有效城市；
• 把成立年份、核准日期解析为年份，按经营状态判定是否退出企业并得到退出年份；
• 剔除成立年份或注册资本缺失、注册资本非正的样本，以及"退出早于成立"的逻辑异常样本。

将企业级明细缓存为 pickle，便于后续复算或单独调试（等价于 R 的 write_rds）：

2. 构建城市×行业×年份存续企业规模面板

build_year_scale(t) 对目标年份 ttt 筛选"存续企业"（成立年份 ≤ t，且未退出或退出年份 > t），按城市 × 行业汇总该年的注册资本总额（output_reg）与企业数（n_firm）。逐年份构建后即得到城市×行业×年份规模面板。

3. 计算 Krugman 式地区产业专业化指数

compute_spec_one_year(df, value_col) 是核心：对某一年的城市×行业规模数据，计算各城市的专业化指数。

关键点：

1. 计算每个城市该口径的总规模 city_total，并求各行业份额 share；
2. 通过 笛卡尔积补全"城市 × 行业"全网格，该城市该行业无企业时份额记为 0；
3. 对每个行业求全部城市份额之和 sum_share，则"其他城市平均份额"为 (sum_share - share)/(J-1)；
4. 计算差的平方 (share - other_avg)^2，跨行业求和即得该城市的 spec。

结果预览

读取结果，展示 2005 年专业化指数最高的若干城市，便于核对：

数据可视化

下面使用 Python 的 matplotlib 绘制三张图表，直观展示专业化指数的整体趋势、年份分布与头部城市。

图1：各城市平均专业化指数随年份变化趋势

计算全部城市在各年份的平均值与中位数，并叠加"平均值 ± 1 倍标准差"阴影带，刻画整体趋势。

图2：各年专业化指数分布（箱线图）

图3：2005 年专业化程度最高的 20 城市

如何参加课程？

购买 RStata 名师讲堂会员即可参加该课程啦（之前的和未来的都可以参加）！

价格：2800/年 或者 4800/长期

购买会员可以从这里下单：https://rstata.duanshu.com/#/card/list/

名师讲堂会员权益：

1. 参加每个月 3～4 次的名师讲堂课程；
2. 参加平台上的其他 R 语言和 Stata 的课程；
3. 以会员折扣价购买我们分享的数据资料（10 元/份）；
4. 课程内外的提问解答服务（课程外的尽量帮忙解决）。

* 如果发票可添加小编微信 r_stata2 （RStata 李老师）开具。如需数据资料，购买后可添加小编微信免费领取数据折扣卡。

更多关于 RStata 会员的更多信息可添加微信号 r_stata2 咨询：

课程主页（点击文末的阅读原文即可跳转）：https://rstata.duanshu.com/#/brief/course/9dfc640436054f4cbba41308f7dc065d