今天给大家分享使用 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）”的企业。
计算步骤概述 整体计算分为以下几个步骤：
读取与清洗：读取各年工商注册数据，仅保留制造业、有效城市代码、必需变量，并界定进入/退出年份；
构建存续面板：对每个目标年份 t，筛选存活企业，按”城市 × 行业”汇总规模，得到城市×行业×年份存续企业规模面板；
计算专业化指数：对每一年，按公式计算各城市的 spec 指数（注册资本口径 + 企业数量口径）；
输出与可视化：导出结果 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（仅首次） # 设置 CRAN 镜像（knit 时 R 处于非交互模式，不会自动选择镜像） options(repos = c(CRAN = "https://mirrors.tuna.tsinghua.edu.cn/CRAN/")) # 仅在尚未安装时才安装，避免每次 knit 都重装 if (!requireNamespace("reticulate", quietly = TRUE)) { install.packages("reticulate") message("reticulate 安装完成！") } else { message("reticulate 已安装，版本：", packageVersion("reticulate")) }
虚拟环境初始化原理（已在 setup chunk 中完成） 本文档的 setup chunk（隐藏运行）包含如下逻辑：
library( reticulate) .venv_name <- ".venv" .venv_python <- virtualenv_python( .venv_name) if ( ! file.exists( .venv_python) ) { virtualenv_create( .venv_name) .venv_python <- virtualenv_python( .venv_name) } Sys.setenv( RETICULATE_PYTHON = .venv_python) use_virtualenv( .venv_name, required = TRUE )
这样做的关键在于：knitr 在处理第一个 {python} chunk 时，reticulate 已经通过 RETICULATE_PYTHON 环境变量知道要使用 .venv，不会再去碰 Anaconda。
在虚拟环境中安装 Python 包（仅首次） py_pkgs <- c ( "numpy" , "pandas" , "matplotlib" , "geopandas" ) installed <- py_list_packages( ".venv" ) $ package need_install <- setdiff( py_pkgs, installed) if ( length ( need_install) > 0 ) { virtualenv_install( ".venv" , packages = need_install) message( "已安装缺失的包：" , paste( need_install, collapse = ", " ) ) } else { message( "所有 Python 包已就绪，无需安装" ) }
验证激活状态 # 验证当前绑定的 Python 路径（应指向 .venv 目录） py_config()
查看已安装的包 pkgs <- py_list_packages( ".venv" ) key_pkgs <- c ( "numpy" , "pandas" , "matplotlib" , "geopandas" ) pkgs[ pkgs$ package %in% key_pkgs, c ( "package" , "version" ) ]
虚拟环境管理常用命令 # 查看所有已创建的虚拟环境 virtualenv_list() # 删除虚拟环境（当不再需要时） # virtualenv_remove(".venv") # 升级某个包 # virtualenv_install(".venv", packages = "matplotlib", ignore_installed = TRUE)
详细计算代码 下面按步骤完整展示 Python 代码（基于 pandas / numpy），每一段均可在 reticulate 管理的虚拟环境中直接运行。代码与同目录下的 01_测算地区产业专业化指标.py / 02_可视化.py 一致。
0. 路径与参数 首先设定工程目录、输出目录、数据目录，以及制造业行业代码与所需变量。
import os from pathlib import Path import numpy as np import pandas as pd # ---- 路径与参数 ---------------------------------------------------------- # proj_dir / out_dir / data_dir：knitr 的 {python} chunk 工作目录即 Rmd 所在目录 PROJ_DIR = Path(os.getcwd()) OUT_DIR = PROJ_DIR / "输出" # 本地化数据：项目内已裁剪为 9 个必需变量的样本（工商注册信息2025-sample ） DATA_DIR = PROJ_DIR / "工商注册信息2025-sample" FILE_YEARS = list (range (2000, 2006)) # 使用的原始数据文件（按成立年份分年存储） TARGET_YEARS = list (range (2000, 2006)) # 需要测算专业化指数的年份 OUT_DIR.mkdir (parents=True, exist_ok=True) # 制造业大类：C13-C43，剔除 C39（计算机、通信和其他电子设备制造业） mfg_codes = [f"C{n:02d}" for n in range (13, 44) if n != 39] # 仅读取所需列，显著降低内存占用 need_cols = ["注册资本" , "实缴资本" , "行业门类" , "行业大类代码" , "经营状态" , "成立年份" , "核准日期" , "市" , "市代码" ]
1. 读取并清洗单个年份文件 read_one_year() 负责把一年的 CSV 读入并清洗为”企业级”明细：
仅保留行业门类 == “制造业”且大类在 C13–C43 且非 C39 的记录；
仅保留市代码为 6 位数字行政区划码的有效城市；
把成立年份、核准日期解析为年份，按经营状态判定是否退出企业并得到退出年份；
剔除成立年份或注册资本缺失、注册资本非正的样本，以及”退出早于成立”的逻辑异常样本。
# 此处代码需要下载讲义材料查看～ print("步骤 1/4：读取与清洗原始数据 ……") firms = pd.concat([read_one_year(yr) for yr in FILE_YEARS], ignore_index=True) print(f" 清洗后制造业企业记录：{len(firms):,}") firms
将企业级明细缓存为 pickle，便于后续复算或单独调试（等价于 R 的 write_rds）：
# 缓存企业级明细，便于复算（pickle 无需额外依赖） firms.to_pickle(OUT_DIR / "firms_manufacturing_2000_2005.pkl")
2. 构建城市×行业×年份存续企业规模面板 build_year_scale(t) 对目标年份 t 筛选”存续企业”（成立年份 ≤ t，且未退出或退出年份 > t），按城市 × 行业汇总该年的注册资本总额（output_reg）与企业数（n_firm）。逐年份构建后即得到城市×行业×年份规模面板。
print ("步骤 2/4：构建城市×行业×年份存续企业规模面板 ……" )def build_year_scale (t ): sub = firms[(firms["entry_year" ] <= t) & (firms["exit_year" ].isna() | (firms["exit_year" ] > t))] res = (sub.groupby(["city_code" , "city" , "industry" ], dropna=False ) .agg(output_reg=("cap_reg" , "sum" ), n_firm=("cap_reg" , "size" )) .reset_index()) res.insert(0 , "year" , t) return res city_ind_year = pd.concat([build_year_scale(t) for t in TARGET_YEARS], ignore_index=True ) city_ind_year.to_csv(OUT_DIR / "城市_行业_年份_制造业规模.csv" , index=False ) city_ind_year
3. 计算 Krugman 式地区产业专业化指数 compute_spec_one_year(df, value_col) 是核心：对某一年的城市×行业规模数据，计算各城市的专业化指数。
关键点：
计算每个城市该口径的总规模 city_total，并求各行业份额 share；
通过 笛卡尔积 补全”城市 × 行业”全网格，该城市该行业无企业时份额记为 0；
对每个行业求全部城市份额之和 sum_share，则”其他城市平均份额”为 (sum_share - share)/(J-1)；
计算差的平方 (share - other_avg)^2，跨行业求和即得该城市的 spec。
print("步骤 3/4：计算地区产业专业化指数 ……") # 此处代码需要下载讲义材料查看～ # 合并结果 spec_all = spec_reg.merge(spec_cnt[["year", "city_code", "spec_count"]], on=["year", "city_code"], how="left") spec_all = spec_all[["year", "city_code", "city", "n_city", "spec_reg", "spec_count"]] spec_all = spec_all.sort_values(["year", "spec_reg"], ascending=[True, False]).reset_index(drop=True) spec_all.to_csv(OUT_DIR / "地区产业专业化指数_2000_2005.csv", index=False) print("步骤 4/4：输出结果 ……") print(" 各年城市数(J)：") print(spec_all.groupby("year")["n_city"].first().to_string()) print("\n2005 年产业专业化指数最高的 10 个城市：") top = spec_all[spec_all["year"] == 2005].nlargest(10, "spec_reg") print(top[["city", "spec_reg", "spec_count"]].to_string(index=False)) print(f"\n完成！结果已写入：{OUT_DIR}")
结果预览 读取结果，展示 2005 年专业化指数最高的若干城市，便于核对：
spec_all = pd.read_csv( OUT_DIR / "地区产业专业化指数_2000_2005.csv" ) print( "=== 2005 年产业专业化指数最高的 15 个城市 ===" ) print( spec_all[ spec_all[ "year" ] == 2005 ] .nlargest( 15 , "spec_reg" ) [[ "year" , "city" , "n_city" , "spec_reg" , "spec_count" ] ] .to_string( index= False) )
数据可视化 下面使用 Python 的 matplotlib 绘制三张图表，直观展示专业化指数的整体趋势、年份分布与头部城市。
图1：各城市平均专业化指数随年份变化趋势 计算全部城市在各年份的平均值与中位数，并叠加”平均值 ± 1 倍标准差”阴影带，刻画整体趋势。
此处代码需要下载讲义材料查看～
图2：各年专业化指数分布（箱线图） PALETTE = cmp.get_discrete_colors("acton", n=6) fig, ax = plt.subplots(figsize=(10, 6)) years2 = sorted(spec_all["year"].unique()) data = [spec_all[spec_all["year"] == y]["spec_reg"].values for y in years2] bp = ax.boxplot(data, patch_artist=True, flierprops=dict(alpha=0.3)) ax.set_xticks(range(1, len(years2) + 1)) ax.set_xticklabels(years2) for patch, c in zip(bp["boxes"], PALETTE): patch.set_facecolor(c) patch.set_alpha(0.85) for med in bp["medians"]: med.set_color("black"); med.set_linewidth(1.0) ax.set_xlabel("") ax.set_ylabel("产业专业化指数 spec") ax.grid(axis="y", linestyle="--", linewidth=0.5, alpha=0.4) cmp.add_title_and_subtitle(ax, "各年份城市产业专业化指数分布 (2000-2005)", "箱线图展示全部城市的专业化指数分布及其演变", title_fontsize=15, subtitle_fontsize=9, title_pad=32) cmp.add_caption(ax, "数据处理 & 绘图：微信公众号 RStata", fig=fig, fontsize=8) fig.savefig(OUT_DIR / "图2_专业化指数分布.png", dpi=300, bbox_inches="tight", facecolor="white") plt.close(fig)
图3：2005 年专业化程度最高的 20 城市 CMAP_SEQ = cmp.get_scico_colors_cmap("acton") top20 = (spec_all[spec_all["year"] == 2005] .nlargest(20, "spec_reg") .sort_values("spec_reg")) fig, ax = plt.subplots(figsize=(10, 7)) vals = top20["spec_reg"].values norm = (vals - vals.min()) / (vals.max() - vals.min() + 1e-12) bar_colors = [CMAP_SEQ(0.25 + 0.7 * t) for t in norm] ax.barh(top20["city"], top20["spec_reg"], color=bar_colors, alpha=0.95) ax.set_xlabel("产业专业化指数 spec") ax.set_ylabel("") ax.grid(axis="x", linestyle="--", linewidth=0.5, alpha=0.4) cmp.add_title_and_subtitle(ax, "2005 年制造业产业专业化程度最高的 20 个城市", "颜色越深表示专业化指数越高", title_fontsize=15, subtitle_fontsize=9, title_pad=32) cmp.add_caption(ax, "数据处理 & 绘图：微信公众号 RStata", fig=fig, fontsize=8) fig.savefig(OUT_DIR / "图3_2005专业化top20.png", dpi=300, bbox_inches="tight", facecolor="white") plt.close(fig)
点击这里跳转到 RStata 短书平台获取附件：名师讲堂｜使用 Python 测算地区产业专业化指标
评论