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

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

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

与 R 和 Python 版本不同，Stata 实现使用 Mata（Stata 的矩阵编程语言）处理反射法迭代中的矩阵运算。Mata 支持矩阵乘法、转置和逐元素操作，能高效完成二元矩阵 M 的构建与迭代计算。

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

指标背景

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

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

计算步骤

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

1. 数字经济专利识别：通过专利唯一标识 newipzlid 与数字经济专利库进行匹配，筛选出属于数字经济领域的专利；
2. 样本筛选：剔除建筑业、金融业、房地产业企业，剔除 ST、*ST、PT 股票；
3. 专利去重与存量构建：对专利进行两轮去重（公开公告号去重 + 申请号去重），然后以公司-专利-IPC 组合为基本单位，构建累计至每年的专利存量池（含未授权专利）；
4. RCA（显性比较优势）计算：对每个企业-技术组合计算 RCA 指标，判断企业在特定技术领域是否具有比较优势；
5. 反射法迭代：基于 RCA 构建二元矩阵 M，在 Mata 中通过 20 次迭代计算企业层面的数字技术复杂度，并检测迭代坍缩；
6. 标准化：将复杂度指标标准化至 [0,1] 区间，便于跨年比较和回归分析。

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

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

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

核心公式

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

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

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

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

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

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

具体迭代过程：

第 0 轮：数一数每家企业会多少种技术（多样性），以及每种技术有多少家企业会（稀有性）。

第 1 轮：用技术的稀有性重新算企业——如果你会的技术恰好都是”稀有技术”，你的得分就高；如果你会的都是”烂大街的技术”，得分就低。

第 2 轮：用更新后的企业得分重新算技术——如果被高分企业掌握，这项技术就是”稀有技术”；如果被低分企业掌握，这项技术就不算稀有。

第 3～20 轮：重复上述过程，直到数值稳定。

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

数据准备

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

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

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

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

本次演示使用 2010～2016 年的子集数据。首先看一下数据结构：

library(haven)
library(tidyverse)

# 读取样本数据的前几行
read_dta("数字经济筛选标准/基于主分类号筛选的数字经济专利newipzlid.dta",
         n_max = 5)

# 读取样本 CSV
read_csv("上市公司与专利样本数据.csv", n_max = 10,
         col_types = cols(.default = col_character()))

数据包含以下字段：

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

详细计算代码

完整的 Stata 执行文件为 main.do，可直接在 Stata 中运行。

下面逐步讲解完整的计算过程。

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

首先从 DTA 文件中提取所有数字经济专利的 newipzlid 列表，并保存为临时文件供后续 merge 使用：

*- ── 配置 ──
local dta_path       "数字经济筛选标准/基于主分类号筛选的数字经济专利newipzlid.dta"
local csv_path       "上市公司与专利样本数据.csv"
local industry_dta   "2000～2024年国泰安上市公司行业信息.dta"
local min_year       2010
local max_year       2016
local n_iter_k       20       /*- 企业复杂度迭代次数 */
local n_iter_t       19       /*- 技术复杂度迭代次数 */
local collapse_thresh 1.0

*- ── Step 1: 读取数字经济专利 DTA ──
use "`dta_path'", clear
keep newipzlid

*- 统一 newipzlid 为字符串
*- CRITICAL: 使用 %15.0f 避免科学记数法导致不同 ID 映射到相同字符串（merge 误匹配）
format newipzlid %15.0f
capture tostring newipzlid, replace force
replace newipzlid = strtrim(newipzlid)
drop if missing(newipzlid) | newipzlid == ""

*- 去重
duplicates drop newipzlid, force

*- 保存为临时文件（后续 merge 使用）
tempfile digital_ids
save `digital_ids', replace

说明：capture tostring 先尝试将 newipzlid 转为字符串（如果 DTA 中存储为数值型），strtrim() 去除首尾空白。保存为 tempfile 后，后续可以使用 merge m:1 newipzlid usingdigital_ids’` 进行精确匹配。

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

参照论文的做法，剔除建筑业（E）、金融业（J）、房地产业（K）的企业，以及曾被标记为 ST、*ST、PT 的公司：

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

说明：ST/PT 识别使用正则表达式 ustrregexm(股票简称, "ST|PT", 1)，大小写不敏感。行业筛选取每个公司最新的行业代码首字母判断。只要公司在任何一年被标记为 ST/*ST/PT，该公司所有年份数据都会被剔除。

Step 3：读取 CSV 并筛选数据

*- ── Step 3: 读取 CSV ──
*- stringcols(2): newipzlid 作为字符串读入，避免 double 精度丢失
import delimited using "`csv_path'", encoding("utf-8") stringcols(2) clear

*- 筛选年份（2010-2016）
keep if 年份 >= `min_year' & 年份 <= `max_year'

*- 剔除 ST/PT 和 E/J/K 行业公司
merge m:1 股票代码 using `exclude_stocks', keep(master) nogen

*- 保留所有专利申请（论文使用申请专利，不筛选授权状态）

*- newipzlid 清理空格
replace newipzlid = strtrim(newipzlid)

*- merge 匹配数字经济专利
merge m:1 newipzlid using `digital_ids', keep(match) nogen

Step 4：专利去重

*- ── Step 3: 专利去重 ──
*- 此处代码需下载讲义材料查看～

说明：ustrregexra() 是 Stata 的正则替换函数，支持 Unicode。[A-Za-z]$ 匹配末尾的单个字母（如 A、B、U、S）。

Step 5：提取 IPC 四位码

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

*- ── Step 4: 提取 IPC 四位码 ──
gen ipc4 = substr(IPC主分类, 1, 4)
keep if !missing(ipc4) & ipc4 != "" & length(ipc4) >= 4

*- 清理不再需要的列
drop 公开公告号 公开公告号_clean 申请号 专利类型 IPC主分类 授权公告日

Step 6：构建专利存量池

每个公司-专利-IPC 组合只计一次，取最早年份作为进入专利池的时间：

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

说明：collapse (min) 年份, by(...) 相当于 R 中的 df_raw[, .(entry_year = min(年份)), by = c("股票代码", "newipzlid", "ipc4")]。

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

这是核心步骤。对每一年 t，在 Stata 和 Mata 中交替执行。Mata- 负责矩阵运算，Stata- 负责数据整理。

流程为：Stata 数据准备 → Mata 矩阵迭代 → Stata 结果收集。

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

Step 8：合并结果

*- ── Step 7: 合并 ──
clear
foreach f of local year_files {
    capture append using `f'
}
duplicates drop 股票代码 年份, force

Step 9：多种标准化方式

本讲义同时计算四种标准化方式以供对比：

• 方式 1（默认）：逐年 Min-Max，配合年份固定效应使用
• 方式 2：逐年百分位排名（均值恒为 0.5，丢失时间趋势）
• 方式 3：全局 Min-Max（所有年份共用同一尺度）
• 方式 4：面板 z-score → 缩放到 [0,1]

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

Step 10：描述性统计

结果合并后、保存前，查看整体和各年的统计特征：

*- ── 2010-2016 年 complexity 总体统计 ──
summarize complexity, detail

*- ── 逐年统计 ──
bysort 年份: summarize complexity

Step 11：保存结果

*- ── 保存 ──
*- 股票代码格式化为 6 位（前补 0）
tostring 股票代码, replace
replace 股票代码 = string(real(股票代码), "%06.0f")

label variable complexity         "数字技术复杂度 (Min-Max 标准化)"
label variable complexity_minmax  "数字技术复杂度 (逐年 Min-Max)"
label variable complexity_rank    "数字技术复杂度 (逐年百分位排名)"
label variable complexity_global  "数字技术复杂度 (全局 Min-Max)"
label variable complexity_panel   "数字技术复杂度 (面板 z-score → [0,1])"
label variable complexity_z       "数字技术复杂度 (面板 z-score)"
label data "数据处理：微信公众号 RStata"

save "2010~2016年上市公司数字技术复杂度_iter20.dta", replace
export delimited 股票代码 年份 complexity complexity_minmax complexity_rank ///
    complexity_global complexity_panel complexity_z complexity_raw ///
    using "2010~2016年上市公司数字技术复杂度_iter20.csv", replace

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

现象描述

在使用逐年 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.01），这就是迭代坍缩。

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

逐年 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. 应对建议

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

方案	具体做法	优点	缺点
A. 删除坍缩年份	从回归样本中剔除 2013、2015 年	干净，无测量误差	损失 2 年数据（2010-2016 子样本中占比较高）
B. 使用全样本（2008-2024）	坍缩被更多年份”稀释”，年度跳跃相对不明显	与论文做法一致，样本量大	需要重新跑全量计算
C. 加入坍缩哑变量交互	在回归中加入 complexityi,t×I(collapset)	控制坍缩年份的系统性偏差	消耗自由度，解释复杂
D. 用熵权法结果做稳健性检验	以熵权法计算的复杂度作为替代指标重新跑回归	不受迭代坍缩影响，结论更可靠	指标经济含义与反射法略有不同

推荐策略：

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

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

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