名师讲堂｜使用 Python 测算各城市数字产业集聚程度

3589 字 约 12 分钟 名师讲堂 区域经济Python

今天给大家分享使用 Python 测算各城市数字产业集聚程度的方法。该方法参考屠西伟、史丹（2025）《数字产业集聚与企业能源效率改进》，通过区位熵来综合测度城市的数字产业集聚水平。

附件中提供了该参考文献的 PDF 文件，感兴趣的小伙伴可以阅读原文。

指标来源与计算原理

数字产业集聚度（Location Quotient）

区位熵的经济含义

  • DL > 1：该城市数字产业集聚度高于全国平均水平，具有相对专业化优势
  • DL = 1：与全国平均水平相当
  • DL < 1：低于全国平均水平

两种测算方法

本文介绍两种测算方式，主要区别在于分子分母的衡量单位不同：

方法 Xct Sct 优点 局限
注册资本版（论文方法） 数字产业注册资本（万元） 全部企业注册资本（万元） 反映资本密度，与论文一致 大城市分母稀释效应明显
企业数量版（备选方法） 数字产业企业数量（家） 全部企业数量（家） 不受极值影响，城市间对比更直观 无法区分大企业与小企业的贡献

计算步骤概述

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

  1. 读取行业分类：加载《数字经济及其核心产业统计分类（2021）》代码表
  2. 构建注销查找表：从注销企业 CSV 中提取 newgcid → exit_year 映射
  3. 单年聚合：逐年读取注册企业 CSV，先过滤已注销企业，再按城市聚合
  4. 面板累计：跨年累加，得到各城市各年的存量企业指标
  5. 计算区位熵：按公式计算 DL，并进行 Winsorize 极端值处理
  6. 输出结果：保存为 .dta 文件

数据说明

数据来源

  • 工商注册信息：

1949～2023 年工商企业注册信息数据（含经纬度及其所属的省市区县）(版本2)：https://rstata.duanshu.com/#/brief/course/6d38a3f10cdb467492f3204d1ebdd313

  • 注销企业信息：

1970～2023 年各年各省市区县、各行业注销公司工商信息及数量统计面板数据：https://rstata.duanshu.com/#/brief/course/bcdf21ad0e614645b8449e69342e0851

  • 数字经济核心产业分类：数字经济及其核心产业统计分类.dta，提取自《数字经济及其核心产业统计分类（2021）》。

注销企业处理逻辑

本文采用个体层面过滤的方法处理注销企业：

第 t 年存量企业 = t 年及之前注册的 且 t 年及之前未注销的企业

具体实现：预先从注销企业 CSV 提取 newgcid → exit_year 查找表，缓存为 pickle 文件，在每年聚合前关联该表，直接过滤掉 exit_year <= t 的企业，再对存活企业进行聚合。

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

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

重要说明（避免”已初始化”报错）：reticulate 在 R 会话中只能绑定一次 Python——一旦某个 {python} 代码块运行，Python 解释器就被锁定，之后再调用 use_virtualenv() 会报错。因此，虚拟环境的激活必须在所有 {python} 代码块之前完成。本文档的解决方案是在 setup chunk 中通过 Sys.setenv(RETICULATE_PYTHON = ...) 提前锁定 Python 路径。

安装 reticulate（仅首次）

# 设置 CRAN 镜像
options(repos = c(CRAN = "https://mirrors.tuna.tsinghua.edu.cn/CRAN/"))

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)
}

# 通过环境变量抢先锁定 Python（优先级最高）
Sys.setenv(RETICULATE_PYTHON = .venv_python)
use_virtualenv(.venv_name, required = TRUE)

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

py_pkgs <- c("numpy", "pandas", "matplotlib")
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 包已就绪，无需安装")
}

验证激活状态

py_config()

步骤一：读取数字经济产业行业代码

加载分类标准

import pandas as pd
import numpy as np

digi_dta = pd.read_stata("数字经济及其核心产业统计分类.dta")
digi_dta = digi_dta[digi_dta['国民经济行业代码'].notna() & (digi_dta['国民经济行业代码'] != "")]

# 提取 3 位前缀（匹配 CSV 行业中类代码后 3 位数字）
digi_prefixes = set(digi_dta['国民经济行业代码'].str[:3].unique())
# 构建完整 4 位代码集合（备用于小类匹配）
digi_codes_4d = set(digi_dta['国民经济行业代码'].unique())

《数字经济及其核心产业统计分类（2021）》使用 4 位行业代码，而工商注册信息 CSV 中的行业代码字段格式为 I641（中类）或 I6411（小类），均带有字母前缀。

代码的匹配策略如下：

行业字段 提取规则 与分类标准对比
行业小类代码（4位数字） code[1:5] 与 digi_codes_4d 完全匹配
行业中类代码（3位数字） code[1:4] 与 digi_prefixes 前缀匹配

优先使用小类代码，兜底使用中类代码，确保最大覆盖率。

步骤二：构建注销查找表

注销企业处理是整个计算中最关键的一步。

预构建查找表

此处代码需下载讲义材料查看～

代码要点：

  • usecols=[‘newgcid’, ‘退出日期’]：只读两列，大幅节省内存
  • drop_duplicates(subset=[‘newgcid’], keep=’first’)：每个企业只保留最早的注销记录
  • to_pickle(exit_pkl)：缓存为 pickle 文件，多进程工人从磁盘按需加载，避免跨进程序列化开销

步骤三：单年聚合函数

函数设计

这是计算流程的核心函数，逻辑为：读取 → 过滤注销 → 识别数字产业 → 企业缩尾（注册资本版）→ 按城市聚合。

注册资本版与企业数量版的唯一区别在于：注册资本版在聚合前需做企业层面 Winsorize，企业数量版（计数）无需缩尾。以下以注册资本版为例说明。

此处代码需下载讲义材料查看～

五个步骤解析：

步骤 核心代码 说明
① 读取 pd.read_csv(..., usecols=...) 只读必要列，单年文件可达 200MB+，列选择可节省 60% 内存
② 过滤注销 dt[dt['exit_year'].isna() | (dt['exit_year'] > year_val)] NaN 表示未注销；> year_val 表示当年末尚未退出
③ 行业识别 dt['行业小类_数字'].isin(digi_codes_4d) isin() 等价于 R 的 %in%，布尔值自动转 0/1
④ 企业缩尾 .clip(lower=p05, upper=p95) 先于聚合，在企业个体层面截断极端注册资本，防止天价壳公司拉偏城市
⑤ 聚合 dt['注册资本'] * dt['is_digi'].astype(int) 利用布尔→0/1转换，sum(注册资本 * is_digi) 即数字产业注册资本

并行执行

此处代码需下载讲义材料查看～

ProcessPoolExecutor.map() 等价于 R 的 future_map()，自动将任务分发到多个子进程并行执行。partial() 将固定参数绑定为偏函数，只将年份作为可变参数传入 map()。

关键区别（Python vs R）：

R Python 说明
future::plan(multisession) ProcessPoolExecutor(max_workers=N) 并行引擎
furrr::future_map() executor.map() 并行分发
saveRDS() / readRDS() to_pickle() / read_pickle() 大对象缓存
options(future.globals.maxSize) pickle 磁盘缓存 + 进程内按需加载 避免跨进程序列化

步骤四：构建面板与计算累计存量

补全面板并计算累计值

注册数据是流量（某年新注册的企业），而区位熵需要存量（截至该年仍存活的企业总量）。由于已在聚合前过滤了注销企业，这里只需按城市做累计加总即可。

此处代码需下载讲义材料查看～

关键函数对照：

R 代码 Python 代码 含义
CJ(a, b) pd.MultiIndex.from_product([a, b]) 笛卡尔积，生成所有组合的完整面板
set(dt, i, j, value) df[col].fillna(0) 填充缺失值
cumsum(), by=city_code groupby('city_code')[col].cumsum() 按城市分组计算逐年累积

步骤五：计算区位熵与极端值处理

计算 DL

# 每年全国总存量
# 此处代码需下载讲义材料查看～

np.where() 是 Python 的向量化条件函数，等价于 R 的 fifelse()。这里对零分母做安全保护，避免产生 Inf。

Winsorize 处理极端值

本项目采用两层 Winsorize 策略：

  1. 企业层面（步骤三已完成）：对每个年份内的企业注册资本做 5%/95% 截断
  2. 城市层面：对最终 DL 指标做 5%/95% 截断
dl_valid = panel_full['DL'].dropna()
dl_valid = dl_valid[np.isfinite(dl_valid)]

p05 = dl_valid.quantile(0.05)
p95 = dl_valid.quantile(0.95)

panel_full['DL_raw'] = panel_full['DL']   # 保留原始值供对比
panel_full['DL'] = panel_full['DL'].clip(lower=p05, upper=p95)

.clip(lower=p05, upper=p95) 等价于 R 的 pmax(pmin(DL, p95), p05)，一步到位完成 Winsorize 截断。

步骤六：合并城市名称与输出

从注册文件反查城市名称

city_map = pd.DataFrame(columns=['city_code', 'city_name'])
for y in sorted(reg_years, reverse=True):
    fpath = os.path.join(dir_reg, f"{y}.csv")
    if not os.path.exists(fpath):
        continue
    dt = pd.read_csv(fpath, usecols=['市', '市代码'], dtype={'市代码': str},
                     na_values=[''], keep_default_na=True)
    dt.columns = ['city_name', 'city_code']
    dt = dt[dt['city_code'].notna() & (dt['city_code'] != "") &
            dt['city_name'].notna() & (dt['city_name'] != "") & (dt['city_name'] != " ")]
    dt = dt.drop_duplicates(subset=['city_code'], keep='first')
    dt = dt[~dt['city_code'].isin(city_map['city_code'])]
    if len(dt) > 0:
        city_map = pd.concat([city_map, dt], ignore_index=True)

output = output.merge(city_map, left_on='citycode', right_on='city_code', how='left')

保存为 DTA 文件

# Stata 15+ (version=118) 支持 Unicode 标签
output['city_name'] = output['city_name'].fillna('')
output.to_stata(
    "数字产业集聚度_注册资本_各城市.dta",
    version=118,
    data_label='数据处理：微信公众号 RStata',
    variable_labels={
        'citycode': '市代码',
        'city_name': '市',
        'year': '年份',
        'stock_cap_total': '存量注册资本_万元',
        'stock_cap_digi': '存量数字企业注册资本_万元',
        'DL_raw': 'DL_原始',
        'DL': '数字产业集聚度',
    },
    write_index=False
)

Python 的 to_stata() 使用 version=118（Stata 15+ 格式），支持 Unicode 变量标签。variable_labels 参数设置中文标签，在 Stata 中通过 describe 可以看到。列名必须为 ASCII，因此用英文列名 + 中文标签的组合。

以上所有计算步骤对应的完整代码参见附件中的 计算数字产业集聚度_注册资本版.py 和 计算数字产业集聚度_企业数量版.py，可直接在终端运行。

结果展示（2010–2012 示例）

以下直接读取已生成的结果文件进行展示。运行前请先执行上述 Python 脚本生成 DTA 文件。

读取计算结果

import pandas as pd
import numpy as np

# 读取已生成的结果文件
dt_cap = pd.read_stata("数字产业集聚度_注册资本_各城市.dta")
dt_cnt = pd.read_stata("数字产业集聚度_企业数量_各城市.dta")
print(f"注册资本版：{len(dt_cap)} 行，{dt_cap['citycode'].nunique()} 个城市，年份 {int(dt_cap['year'].min())}-{int(dt_cap['year'].max())}")
print(f"企业数量版：{len(dt_cnt)} 行，{dt_cnt['citycode'].nunique()} 个城市，年份 {int(dt_cnt['year'].min())}-{int(dt_cnt['year'].max())}")

描述性统计

print("========== 注册资本版 描述性统计 ==========")
dl_cap = dt_cap['DL'].dropna()
dl_cap = dl_cap[np.isfinite(dl_cap) & (dl_cap >= 0)]
print(f"Mean: {dl_cap.mean():.4f}  Median: {dl_cap.median():.4f}  SD: {dl_cap.std():.4f}")
print(f"Min: {dl_cap.min():.4f}  Max: {dl_cap.max():.4f}  N: {len(dl_cap)}")

print("\n========== 企业数量版 描述性统计 ==========")
dl_cnt = dt_cnt['DL'].dropna()
dl_cnt = dl_cnt[np.isfinite(dl_cnt) & (dl_cnt >= 0)]
print(f"Mean: {dl_cnt.mean():.4f}  Median: {dl_cnt.median():.4f}  SD: {dl_cnt.std():.4f}")
print(f"Min: {dl_cnt.min():.4f}  Max: {dl_cnt.max():.4f}  N: {len(dl_cnt)}")

Top-10 城市（2012 年）

print("========== 注册资本版 Top-10（2012 年）==========")
top10_cap = dt_cap[dt_cap['year'] == 2012].nlargest(10, 'DL')
for _, row in top10_cap.iterrows():
    print(f"  {row['city_name']}  DL={row['DL']:.3f}  注册资本_亿元={row['stock_cap_total']/1e4:.2f}")

print("\n========== 企业数量版 Top-10（2012 年）==========")
top10_cnt = dt_cnt[dt_cnt['year'] == 2012].nlargest(10, 'DL')
for _, row in top10_cnt.iterrows():
    print(f"  {row['city_name']}  DL={row['DL']:.3f}  企业数={int(row['stock_n_total'])}")

两种方法的 Top-10 城市排名存在差异，主要原因如下：

  • 注册资本版更受大体量企业影响——若某城市少数几家高注册资本的数字企业占全国注册资本比例高，区位熵会被抬高
  • 企业数量版更反映数字产业在城市中的普及度，对城市规模中性

区域对比

# 四大区域分类（国家统计局标准）
# 此处代码需下载讲义材料查看～

与论文结果的对比

参考论文（屠西伟、史丹 2025）的结论，数字产业集聚度应呈现 东部 > 中部 > 西部 > 东北 的梯度格局。以示例数据计算的结果中：

  • 一致之处：东部排第一、中部排第二，与论文完全一致
  • 差异之处：东北与西部的排序可能与论文相反，需要全量数据 + GDP 加权才能最终确认

可视化：区域对比

此处代码需下载讲义材料查看～

可视化：Top-20 城市趋势

# 取最新年份 Top-20 城市
latest = int(dt_cnt['year'].max())
top_cities = dt_cnt[dt_cnt['year'] == latest].nlargest(20, 'DL')['citycode'].tolist()
trend_dt = dt_cnt[dt_cnt['citycode'].isin(top_cities)]

fig, ax = plt.subplots(figsize=(11, 7))

for city in top_cities:
    sub = trend_dt[trend_dt['citycode'] == city].sort_values('year')
    name = sub['city_name'].iloc[0] if 'city_name' in sub.columns else city
    ax.plot(sub['year'], sub['DL'], label=name, linewidth=1.1, alpha=0.85)
    ax.scatter(sub['year'], sub['DL'], s=18, alpha=0.85)

ax.set_xlabel('年份', fontproperties=font_prop)
ax.set_ylabel('数字产业集聚度 DL', fontproperties=font_prop)
ax.grid(alpha=0.3)
# 主标题用 fig.suptitle，字号可控
fig.suptitle('数字产业集聚度 Top-20 城市趋势（企业数量版，2010–2012）', fontsize=24, fontproperties=font_prop, y=0.99)
ax.legend(loc='upper right', fontsize=12, ncol=2, prop=font_prop)
fig.text(0.5, 0.01, '数据爬取&绘制：微信公众号 RStata', fontsize=8, ha='center', va='bottom', color='gray', fontproperties=font_prop)
plt.subplots_adjust(top=0.92, bottom=0.10)
plt.show()

点击这里跳转到 RStata 短书平台获取附件：名师讲堂｜使用 Python 测算各城市数字产业集聚程度

相关文章

评论