名师讲堂｜使用 Python 测算上市公司数字技术复杂度

由于存在坍缩问题，所以这个就不打算分享处理好的数据了，大家可以参考代码自行决定如何计算。

今天给大家讲解如何使用 Python 测算上市公司的数字技术复杂度指标。该指标来源于宋晴等（2026）发表在《数量经济技术经济研究》上的论文「数字技术复杂度与企业市场价值——基于策略性创新的视角」。附件中提供了该论文的 PDF 文件。

本文在 R 语言中通过 reticulate 包来调用 Python 代码。setup chunk 会创建并锁定 Python 虚拟环境 .venv，确保所有 {python} 代码块在正确的环境中执行。

数字技术复杂度指标的来源与计算方法

指标背景

数字技术复杂度（Digital Technology Complexity）是衡量企业在数字经济相关技术领域创新能力的综合性指标。与传统的专利数量指标不同，技术复杂度不仅考虑企业拥有的专利数量，还考虑了企业所涉及技术的稀缺性和复杂性——即企业是否集中在少数高门槛技术领域还是分散在普遍性技术领域。

该指标的计算基于 Hidalgo 和 Hausmann（2009）提出的反射法（Method of Reflections），最初用于测算国家出口产品的经济复杂度。宋晴等（2026）将其创新性地应用于企业层面的数字技术专利分析。

计算步骤

整个计算过程可以分为以下 6 个主要步骤：

1. 数字经济专利识别：通过专利唯一标识 newipzlid 与数字经济专利库进行匹配，筛选出属于数字经济领域的专利；
2. 样本筛选：剔除建筑业、金融业、房地产业企业，剔除 ST、*ST、PT 股票（只要有一年出现即剔除该公司所有年份数据）；
3. 专利去重与 IPC 提取：对专利进行两轮去重（公开公告号去重 + 申请号去重），提取 IPC 主分类号四位码作为技术分类基本单位；
4. RCA（显性比较优势）计算：对每个企业-技术组合计算 RCA 指标，判断企业在特定技术领域是否具有比较优势；
5. 反射法迭代：基于 RCA 构建二元矩阵 M，通过迭代计算企业层面的数字技术复杂度；
6. 标准化：将复杂度指标标准化至 [0,1] 区间，便于跨年比较和回归分析。

关于”忽增忽降”问题的重要说明

在使用反射法计算复杂度指标时，可能会遇到一个看似异常的现象：逐年标准化后的复杂度均值在某些年份突然出现大幅上升或下降（例如从 0.6 跳到 1.0，次年又跌回 0.4）。

这并不是计算错误，而是反射法在特定数据条件下的固有特征。本讲义将在后文（”计算结果解读”一节）中详细解释其数学原理和数据层面的根本原因。

核心公式

)

直观理解：反射法是什么？

核心公式涉及较多数学符号，这里用一个通俗类比帮助理解反射法的本质。

问题意识：假设你要评价哪家企业的数字技术”厉害”，最简单的办法是数专利数量——谁专利多谁厉害。但这样有问题：一个企业拥有 1000 件低端专利，未必比拥有 10 件尖端专利的企业更”复杂”。反射法的出发点正是剥离专利规模，只看技术结构。

一个类比：想象一个”谁有稀缺能力”的互评游戏：

• 企业说：我厉害，因为我掌握的技术别人都不会（技术越稀有，我就越厉害）
• 技术说：我厉害，因为会我的都是厉害企业（企业越厉害，我就越稀有）

这两句话互相”反射”——企业的得分取决于它掌握的技术有多稀有，而技术的稀有程度又取决于掌握它的企业有多厉害。信息在企业 ↔ 技术之间来回反射，经过多轮迭代后逐渐收敛到一个稳定值，这就是”反射法”（Method of Reflections）名字的由来。

具体迭代过程：

一句话总结：反射法 = 用”你拥有的东西有多稀缺”来定义”你有多厉害”，并且让这个定义来回迭代直到收敛。它本质上是一种基于网络位置的排序算法，和 Google PageRank 的思想有异曲同工之处——一个企业”厉害”，不仅取决于它自己，还取决于它关联的技术的”质量”。

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

在 R 中通过 reticulate 包来调用 Python，最好的实践是为项目创建一个专属的 Python 虚拟环境，将所需依赖隔离到独立空间。

重要说明（避免”已初始化”报错）：reticulate 在 R 会话中只能绑定一次 Python——一旦某个 {python} 代码块运行，Python 解释器就被锁定。因此，虚拟环境的激活必须在所有 {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（隐藏运行）包含如下逻辑：

# 以下展示 setup chunk 中已执行的逻辑（eval=FALSE，仅供参考，不重复执行）

library(reticulate)

.venv_name <- ".venv"
.venv_python <- virtualenv_python(.venv_name)

# 虚拟环境不存在时自动创建
if (!file.exists(.venv_python)) {
  virtualenv_create(.venv_name, packages = c("numpy", "pandas", "scipy", "pyreadstat"))
  .venv_python <- virtualenv_python(.venv_name)
}

# 通过环境变量抢先锁定 Python（优先级最高，早于任何 {python} chunk）
Sys.setenv(RETICULATE_PYTHON = .venv_python)
use_virtualenv(.venv_name, required = TRUE)

检查 Python 包是否已就绪

# 检查关键包是否已安装在虚拟环境中
py_pkgs <- c("numpy", "pandas", "scipy", "pyreadstat")
tryCatch({
  installed <- py_list_packages(".venv")$package
  need_install <- setdiff(py_pkgs, installed)
  if (length(need_install) > 0) {
    message("以下包缺失，正在自动安装：", paste(need_install, collapse = ", "))
    virtualenv_install(".venv", packages = need_install)
    message("安装完成！")
  } else {
    message("所有 Python 包已就绪 (numpy, pandas, scipy, pyreadstat)")
  }
}, error = function(e) {
  message("包检查跳过：", conditionMessage(e))
})

验证激活状态

# 验证当前绑定的 Python 路径
py_config()

查看已安装的包

# 列出当前 Python 环境中已安装的关键包
tryCatch({
  pkgs <- py_list_packages(".venv")
  key_pkgs <- c("numpy", "pandas", "scipy", "pyreadstat")
  pkgs[pkgs$package %in% key_pkgs, c("package", "version")]
}, error = function(e) {
  message("包列表查询跳过（Python 尚未初始化）")
  invisible(NULL)
})

虚拟环境管理常用命令

# 以下命令仅供参考（eval=FALSE，不执行）
# 查看所有已创建的虚拟环境
virtualenv_list()

# 删除虚拟环境（当不再需要时）
# virtualenv_remove(".venv")

# 升级某个包
# virtualenv_install(".venv", packages = "scipy", ignore_installed = TRUE)

---

数据准备

我们需要准备以下三组数据：

1. 数字经济专利库：基于主分类号筛选的数字经济专利newipzlid.dta，包含所有被判定为数字经济领域的专利 newipzlid；
2. 上市公司与专利数据：1985～2024年上市公司与专利数据匹配结果.csv，包含每家上市公司每年申请的所有专利信息；
3. 上市公司行业信息：2000～2024年国泰安上市公司行业信息.dta，包含每个公司的行业代码和股票简称，用于样本筛选。

1985～2024 年上市公司与专利数据匹配结果（版本3，含申请、授权信息）：https://rstata.duanshu.com/#/brief/course/04100321f88b411f90429be934bff

1985～2024 年数字经济产业相关专利筛选结果：https://rstata.duanshu.com/#/brief/course/d5dfb9ca0858457ebc4f176fce9fee80

本讲义使用 2010～2016 年的子样本进行演示。附件中提供了提取好的样本数据 上市公司与专利样本数据.csv，其变量与完整数据相同：

变量名	含义
股票代码	上市公司的股票代码
newipzlid	专利唯一标识 ID
年份	专利申请年份
专利类型	专利类型（发明申请、发明授权、实用新型、外观设计）
公开公告号	专利公开公告号
申请号	专利申请号
IPC主分类	专利的 IPC 主分类号
授权公告日	专利授权公告日

首先来看一下样本数据的基本结构：

import pandas as pd

# 读取样本数据
df_sample = pd.read_csv("上市公司与专利样本数据.csv", dtype={'股票代码': str})
df_sample

样本筛选说明

参照论文的做法，对样本进行以下筛选：

1. 行业筛选：根据证监会 2012 版行业代码，剔除建筑业（代码首字母 E）、金融业（J）、房地产业（K）的企业。这些行业与数字技术创新关联较弱，且专利行为与制造业有本质差异；
2. ST 筛选：剔除曾经被标记为 ST、*ST 或 PT 的公司。这些公司处于财务异常状态，其专利申请行为可能不代表正常的创新决策。只要该公司在任何一年被标记为 ST/*ST/PT，就剔除该公司所有年份的数据。

行业信息 DTA 的关键变量：

变量名	含义
股票代码	上市公司股票代码
股票简称	股票简称（含 ST/*ST/PT 标记）
行业代码	证监会 2012 版行业分类代码
统计截止日期	行业信息统计截止日期

详细计算代码

下面我们逐步讲解完整的计算过程。Python 代码需要加载以下包：

import pandas as pd
import numpy as np
from scipy.sparse import coo_matrix
from datetime import datetime
import sys
import os

Step 1：读取数字经济专利 DTA

首先从数字经济专利 DTA 文件中提取所有数字经济专利的 newipzlid 列表：

# 配置路径与参数
dta_path  = "数字经济筛选标准/基于主分类号筛选的数字经济专利newipzlid.dta"
csv_path  = "上市公司与专利样本数据.csv"
industry_dta = "2000～2024年国泰安上市公司行业信息.dta"
min_year  = 2010
max_year  = 2016
n_iter_k  = 20    # 企业复杂度迭代次数
n_iter_t  = 19    # 技术复杂度迭代次数
collapse_thresh = 0.06  # Kc range 小于此值视为迭代坍缩

# 行业剔除：证监会2012版行业代码首字母
#   E = 建筑业, J = 金融业, K = 房地产业
exclude_industry_letters = ["E", "J", "K"]

print("=" * 60)
print("  计算上市公司历年数字技术复杂度")
print(f"  方法：反射法（Method of Reflections）")
print(f"  参考：宋晴等（2026）")
print(f"  年份范围：{min_year}-{max_year}")
print(f"  迭代次数：k={n_iter_k}, t={n_iter_t}")
print(f"  坍缩检测阈值：Kc range < {collapse_thresh}")
print(f"  行业剔除：{','.join(exclude_industry_letters)}")
print(f"  ST/PT 剔除：是")
print("=" * 60 + "\n")

# 读取 DTA，只提取 newipzlid 列
# 使用 columns 参数等价于 R 的 read_dta(col_select = "newipzlid")
dta_df = pd.read_stata(dta_path, columns=['newipzlid'])

# 提取唯一 newipzlid
digital_patent_ids = dta_df['newipzlid'].dropna().unique()
del dta_df

print(f"  DTA 唯一 newipzlid 数量: {len(digital_patent_ids):,}")

Step 2：读取行业信息，构建剔除名单

从行业信息 DTA 中识别需要剔除的公司：

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

上市公司行业信息来源于之前分享的这个数据：

2000～2024 年上市公司注册地址与办公地址（含经纬度及其所处的省市区县）
：https://rstata.duanshu.com/#/course/746d4f595eba41a2a116b5a6edc3ac0a

Step 3：读取 CSV 并筛选数据

读取上市公司与专利数据，依次进行年份筛选、公司筛选和数字经济专利匹配。

说明：论文使用上市公司申请的数字技术专利（不筛选授权状态），因此保留所有专利申请记录。

print("\n【Step 3】读取 CSV 并筛选数据...")

# 只读取需要的列
cols_use = ["股票代码", "newipzlid", "年份", "专利类型",
            "公开公告号", "申请号", "IPC主分类", "授权公告日"]

# newipzlid 以字符串读取，避免 pandas 推断为 int64 导致与 DTA str 类型不匹配
# 授权公告日以字符串读取，避免 pandas 自动推断日期类型
df_raw = pd.read_csv(csv_path, usecols=cols_use,
                     dtype={'newipzlid': str, '授权公告日': str, '股票代码': str})

print(f"  CSV 总行数: {len(df_raw):,}")

# --- 3a: 筛选年份（2010-2016）---
df_raw = df_raw[(df_raw['年份'] >= min_year) & (df_raw['年份'] <= max_year)]
print(f"  年份筛选后: {len(df_raw):,}")

# --- 3b: 剔除 ST/*ST/PT 和指定行业公司 ---
n_before = df_raw['股票代码'].nunique()
df_raw = df_raw[~df_raw['股票代码'].isin(exclude_stocks)]
n_after = df_raw['股票代码'].nunique()
print(f"  公司筛选（剔除 ST/PT + 建筑/金融/房地产）:")
print(f"    剔除前公司数: {n_before}")
print(f"    剔除后公司数: {n_after}（剔除 {n_before - n_after} 家）")

# --- 3c: 不筛选授权状态（论文使用申请专利，保留所有专利申请）---
print(f"  保留所有专利申请（含未授权）: {len(df_raw):,}")
print(f"  专利申请类型分布:")
print(df_raw.groupby('专利类型').size().sort_values(ascending=False).to_string())

# --- 3d: 通过 newipzlid 匹配数字经济专利 ---
digital_set = set(digital_patent_ids)
df_raw = df_raw[df_raw['newipzlid'].isin(digital_set)].copy()
del digital_patent_ids, digital_set
print(f"  newipzlid 匹配（数字经济专利）后: {len(df_raw):,}")

Step 4：专利去重

由于专利数据中存在重复的情况，所以分两步进行去重：

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

Step 5：提取 IPC 四位码

以专利主分类号的 IPC 四位码作为技术分类的基本单位：

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

Step 6：逐年计算数字技术复杂度（反射法）

这是核心步骤。对每一年 t，执行以下操作：

1. 从存量数据中筛选 年份 <= t 的专利（简化存量池：同一 newipzlid 不会跨年出现，直接用 df_raw[年份 <= t]）
2. 计算每个企业-技术组合的 RCA
3. 构建 RCA≥1 的二元矩阵 M
4. 使用反射法迭代 20 次（企业）/ 19 次（技术）
5. 得到每个企业的数字技术复杂度 Kc
6. 坍缩检测：检查 Kc 的 range 是否小于阈值（默认 0.06）

★ 关于”迭代坍缩”的说明：
当企业度分布 CV 接近 1.0 时，20 次迭代后所有企业的 Kc 值趋同（range < 0.06），
导致逐年 Min-Max 标准化后均值出现异常跳跃。详见后文”忽增忽降问题的深入解读”章节。

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

Step 7：合并结果

将各年份的计算结果合并为完整的数据框：

print("\n【Step 7】合并结果...")

results = pd.concat(results_list.values(), ignore_index=True)
del results_list

print(f"  总观测数: {len(results):,}")
print(f"  年份范围: {results['年份'].min()}-{results['年份'].max()}")
print(f"  唯一企业数: {results['股票代码'].nunique()}")

results.head(20)

Step 8：多种标准化方式对比

反射法得到的原始 Kc 值在不同年份之间的量纲差异较大，需要进行标准化。本节对比多种标准化方式，并重点解释逐年 Min-Max 标准化中”忽增忽降”现象的成因。

方式 1：逐年 Min-Max（论文使用方式）

此处内容需下载讲义材料查看～

方式 2：逐年百分位排名

此处内容需下载讲义材料查看～

方式 3：全局 Min-Max

此处内容需下载讲义材料查看～

方式 4：面板 z-score 缩放到 [0,1]

此处内容需下载讲义材料查看～

代码实现

此处内容需下载讲义材料查看～

Step 9：描述性统计

print("\n【Step 9】描述性统计...")

print(f"\n── {min_year}～{max_year} 年 complexity 总体统计 ──")
import numpy as np

desc = results['complexity'].agg(['count', 'mean', 'min', 'median', 'max', 'std']).round(4)
desc.index = ['样本量', '均值', '最小值', '中位数', '最大值', '标准差']
print(desc.to_string())

print(f"\n── 逐年统计 ──")
yearly = results.groupby('年份').agg(
    企业数 = ('complexity', 'count'),
    均值   = ('complexity', lambda x: round(x.mean(), 4)),
    中位数 = ('complexity', lambda x: round(x.median(), 4)),
    标准差 = ('complexity', lambda x: round(x.std(), 4))
).reset_index()
print(yearly.to_string(index=False))

“忽增忽降”问题的深入解读

现象描述

在使用逐年 Min-Max 标准化后，你可能会观察到复杂度指标的逐年均值出现大幅波动。例如：

• 2012 年均值 = 0.91（偏高）
• 2013 年均值 = 0.52（大幅下跌）
• 2014 年均值 = 0.94（大幅回升）
• 2015 年均值 = 0.33（再次大幅下跌）

这种”忽增忽降”并非计算错误，而是以下两个因素共同作用的必然结果。

根本原因一：迭代坍缩（Iterative Collapse）

反射法的数学本质是让复杂度值在二分图（企业 ↔ IPC 技术）上不断扩散。当 M 矩阵中企业的度分布（即每个企业拥有的 RCA≥1 的 IPC 种类数）趋于均匀时，迭代会使所有企业的 Kc 值收敛到几乎相同的值。

我们用 2013 年（坍缩年）和 2012 年（正常年）的对比来说明：

指标	2012 年（正常）	2013 年（坍缩）
企业数	2,466	2,806
技术数（IPC四位码）	124	125
M 矩阵密度	0.031	0.034
企业度 CV	1.053	1.011
Kc range	166.4	0.057

企业度 CV（变异系数）= sd / mean，衡量企业间专利多样性的差异程度。当 CV 接近 1.0 时，所有企业的专利组合高度相似，反射法无法区分谁更复杂，迭代后 Kc 值全部挤在一起（range < 0.06），这就是迭代坍缩。

根本原因二：逐年标准化的独立尺度

逐年 Min-Max 每年独立地拉伸数据到 [0,1]。当某年发生迭代坍缩时，Kc 的 range 极小，几乎所有企业的 Kc 值相同，Min-Max 标准化后所有企业都被映射到相同或非常接近的值（通常是 0 或 1 附近），导致该年均值出现异常。

反之，当某年未发生坍缩、Kc 分布较宽时，Min-Max 后能保留较好的区分度，均值处于正常区间。

数据层面的深层原因

为什么企业度 CV 会在某些年份跌破 1.0（趋于均匀）？这反映了产业层面的真实结构变化：

1. 企业数激增：2010 年 1,442 家 → 2015 年 3,284 家（翻倍），大量新进入企业初期的专利组合往往集中在少数热门 IPC 上（如 G06F、H04L）；
2. 技术数基本停滞：2010 年 118 种 → 2015 年 130 种（仅 +10%），但企业数翻倍 → 每项技术被更多企业同时持有，M 矩阵密度上升但结构变”平”；
3. 专利同质化涌入：2013 年前后数字经济政策密集出台，大量企业同时申请类似的数字技术专利，企业间差异缩小。

这不是数据错误，而是反射法算法对产业结构突变的敏感性表现。

解决方案

方案	做法	优点	缺点
A. 减少迭代次数	改用 5 次迭代（而非 20 次）	所有年份都有真实区分度，不会坍缩	与论文的 20 次迭代设定不完全一致
B. 检测坍缩年并自动降级	range < 1.0 时对该年改用 5 次迭代	非坍缩年与论文一致	跨年不可比（不同年用不同迭代次数）
C. 使用 log(Kc) + Min-Max	对原始 Kc 取对数后再标准化	压缩右尾，减轻坍缩影响	改变指标的经济含义

推荐做法：在回归分析中使用逐年 Min-Max 标准化 + 年份固定效应。年份固定效应会吸收所有年份层面的宏观因素（包括产业政策冲击），复杂度系数的识别完全来自”同年企业间的横截面差异”，这是论文的基准做法。

对后续计量回归的影响

上述”忽增忽降”问题是否会污染回归结果，是实证研究中必须正视的问题。以下从计量经济学角度进行系统分析。

1. 测量误差与衰减偏误（Attenuation Bias）

2. 与年份固定效应的交互作用

逐年 Min-Max 标准化配合年份固定效应是论文的基准设定，但坍缩年份存在特殊问题：

年份固定效应吸收的是所有企业共同面对的年份冲击，复杂度系数的识别理论上完全来自”同年企业间的横截面差异”。但当某年横截面差异本身被算法扭曲时，该年份观测对整体估计的贡献是失真的。

3. 动态面板设定下的额外风险

此处内容需下载讲义材料查看～

4. 应对建议

针对上述风险，提出以下四种处理思路：

推荐策略：

1. 基准回归使用全样本（2008-2024）跑反射法复杂度 + 年份固定效应，此时坍缩的影响被 17 个年份大幅稀释，2013 和 2015 只是其中 2 年；
2. 稳健性检验一：删除 2013 和 2015 年重新跑回归，将结果报告于 Robustness Table；
3. 稳健性检验二：用熵权法计算的复杂度作为替代指标，验证核心结论是否依然成立；
4. 若坚持使用 2010-2016 子样本：务必在回归中报告 Cluster-Robust SE，并在附录中披露坍缩年份的描述性统计。

小结：论文本身使用 2008-2024 全样本并成功发表，说明审稿人接受了这一测量方式。若你的研究聚焦于 2010-2016 子样本（如受某政策冲击的影响），则需要在论文中主动讨论坍缩问题，并用上述稳健性检验证明结论不受其干扰。

Step 10：保存结果

将股票代码统一转为 6 位字符串（不足前补 0），然后保存为 CSV 文件：

print("\n【Step 10】保存结果...")

# 股票代码转 6 位字符串，不足前补 0
results['股票代码'] = results['股票代码'].astype(str).str.zfill(6)

# 保存 CSV（用于检查）
output_csv = f"{min_year}～{max_year}年上市公司数字技术复杂度_iter{n_iter_k}_py.csv"
output_cols = ['股票代码', '年份', 'complexity',
              'complexity_minmax', 'complexity_rank',
              'complexity_global', 'complexity_panel', 'complexity_z']

results[output_cols].to_csv(output_csv, index=False, encoding='utf-8-sig')
print(f"  CSV 已保存至: {output_csv}")
print("  数据处理：微信公众号 RStata")

# 显示前几行
results[output_cols].head(20)

关键注意事项

点击这里跳转到 RStata 短书平台获取附件：名师讲堂｜使用 Python 测算上市公司数字技术复杂度