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

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

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

指标背景

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

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

计算步骤

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

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

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

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

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

核心公式

以下公式来源于宋晴等（2026）论文附录”一、复杂度测算与迭代过程”。以 t 年为例，本年度共 I 家上市公司，其数字技术专利存量共涉及 P 类数字技术。

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

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

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

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

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

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

具体迭代过程：

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

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

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

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

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

数据准备

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

数字经济专利库：基于主分类号筛选的数字经济专利newipzlid.dta，包含所有被判定为数字经济领域的专利 newipzlid；
上市公司与专利数据：上市公司与专利样本数据.csv，包含每家上市公司每年申请的所有专利信息（2010～2016 年子样本）；
上市公司行业信息：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 主分类号
`授权公告日`	专利授权公告日

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

library(tidyverse)
library(data.table)
library(haven)

# 读取样本数据
fread("上市公司与专利样本数据.csv") %>%
  as_tibble()

样本筛选说明

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

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

行业信息 DTA 的关键变量：

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

详细计算代码

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

library(data.table)
library(Matrix)
library(tidyverse)
library(haven)

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 <- 1.0  # Kc range 小于此值视为迭代坍缩

# 行业剔除：证监会2012版行业代码首字母
exclude_industry_letters <- c("E", "J", "K")

# 读取 DTA，提取 newipzlid 列
dta_dt <- read_dta(dta_path)
digital_patent_ids <- dta_dt$newipzlid
rm(dta_dt)
gc(verbose = FALSE)

cat("  DTA 唯一 newipzlid 数量:",
    format(length(digital_patent_ids), big.mark = ","), "\n")

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

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

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

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

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

Step 3：读取 CSV 并筛选数据

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

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

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

df_raw <- fread(csv_path, select = cols_use, encoding = "UTF-8",
                 colClasses = c("授权公告日" = "character"))

cat("  CSV 总行数:", format(nrow(df_raw), big.mark = ","), "\n")

# --- 3a: 筛选年份（2010-2016）---
df_raw <- df_raw[年份 %between% c(min_year, max_year)]
cat("  年份筛选后:", format(nrow(df_raw), big.mark = ","), "\n")

# --- 3b: 剔除 ST/*ST/PT 和指定行业公司 ---
n_before <- df_raw[, uniqueN(股票代码)]
df_raw <- df_raw[!股票代码 %in% exclude_stocks]
n_after <- df_raw[, uniqueN(股票代码)]
cat("  公司筛选（剔除 ST/PT + 建筑/金融/房地产）:\n")
cat("    剔除前公司数:", n_before, "\n")
cat("    剔除后公司数:", n_after, "（剔除", n_before - n_after, "家）\n")
rm(exclude_stocks)
gc(verbose = FALSE)

# --- 3c: 不筛选授权状态（论文使用申请专利，保留所有专利申请）---
cat("  保留所有专利申请（含未授权）:", format(nrow(df_raw), big.mark = ","), "\n")
cat("  专利申请类型分布:\n")
print(df_raw[, .(n = .N), by = 专利类型][order(-n)])

# --- 3d: 通过 newipzlid 匹配数字经济专利 ---
df_raw <- df_raw[newipzlid %in% digital_patent_ids]
rm(digital_patent_ids)
gc(verbose = FALSE)
cat("  newipzlid 匹配（数字经济专利）后:",
    format(nrow(df_raw), big.mark = ","), "\n")

Step 4：专利去重

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

# 清洗公开公告号：去掉末尾字母
df_raw[, 公开公告号_clean := str_replace(公开公告号, "[A-Z]$", "")]

# 去重 1: 按股票代码 + 年份 + 公开公告号_clean
n_before <- nrow(df_raw)
df_raw <- df_raw[!duplicated(df_raw,
                              by = c("股票代码", "年份", "公开公告号_clean"))]
cat(sprintf("  去重 1（公开公告号）: %s → %s (-%d)\n",
            format(n_before, big.mark = ","),
            format(nrow(df_raw), big.mark = ","),
            n_before - nrow(df_raw)))

# 去重 2: 按股票代码 + 年份 + 申请号
n_before <- nrow(df_raw)
df_raw <- df_raw[!duplicated(df_raw, by = c("股票代码", "年份", "申请号"))]
cat(sprintf("  去重 2（申请号）:     %s → %s (-%d)\n",
            format(n_before, big.mark = ","),
            format(nrow(df_raw), big.mark = ","),
            n_before - nrow(df_raw)))

Step 5：提取 IPC 四位码

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

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

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

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

从存量数据中筛选年份 <= t 的专利
计算每个企业-技术组合的 RCA
构建 RCA≥1 的二元矩阵 M
使用反射法迭代 20 次（企业）/ 19 次（技术）
得到每个企业的数字技术复杂度 Kc
坍缩检测：检查 Kc 的 range 是否小于阈值（默认 1.0）

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

# --- 6h: 坍缩检测 ---
  kc_range <- max(K_c) - min(K_c)
  is_collapse <- kc_range < collapse_thresh

  if (is_collapse) {
    collapse_years <- c(collapse_years, t)
    cat(sprintf("  ★★★ %d年: 迭代坍缩！Kc range=%.4f < %.1f\n",
                t, kc_range, collapse_thresh))
  }

  # --- 6i: 保存该年结果 ---
  results_list[[as.character(t)]] <- data.table(
    股票代码      = firms,
    年份          = t,
    complexity_raw = K_c
  )

  cat(sprintf("  %d: 企业=%d, 技术=%d, Kc_mean=%.2f, Kc_range=%.2f%s  %.1fs\n",
              t, length(firms), length(techs),
              mean(K_c), kc_range,
              ifelse(is_collapse, " [坍缩]", ""),
              as.numeric(difftime(Sys.time(), t_iter_start, units = "secs"))))
}

if (length(collapse_years) > 0) {
  cat("\n  ★★★ 以下年份发生迭代坍缩（Kc range <", collapse_thresh, ") ★★★\n")
  cat("  坍缩年份:", paste(collapse_years, collapse = ", "), "\n")
  cat("  原因：M 矩阵的企业度分布趋于均匀，反射法迭代后企业间差异丧失。\n")
} else {
  cat("\n  所有年份均未发生迭代坍缩。\n")
}

Step 7：合并结果

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

results <- rbindlist(results_list)
rm(results_list); gc(verbose = FALSE)

cat("  总观测数:", format(nrow(results), big.mark = ","), "\n")
cat("  年份范围:", min(results$年份), "-", max(results$年份), "\n")
cat("  唯一企业数:", results[, uniqueN(股票代码)], "\n")

results %>%
  as_tibble()

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

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

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

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

方式 2：逐年百分位排名

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

方式 3：全局 Min-Max

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

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

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

代码实现

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

Step 9：描述性统计

cat("\n── 2010-2016 年 complexity 总体统计 ──\n")
desc <- results[, .(
  样本量 = .N,
  均值   = round(mean(complexity), 4),
  最小值 = round(min(complexity), 4),
  中位数 = round(median(complexity), 4),
  最大值 = round(max(complexity), 4),
  标准差 = round(sd(complexity), 4)
)]
print(desc)

cat("\n── 逐年统计 ──\n")
yearly <- results[, .(
  企业数   = .N,
  均值     = round(mean(complexity), 4),
  中位数   = round(median(complexity), 4),
  标准差   = round(sd(complexity), 4)
), by = 年份]
print(yearly, row.names = 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.01），这就是迭代坍缩。

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

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

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

数据层面的深层原因

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

企业数激增：2010 年 1,442 家 → 2015 年 3,284 家（翻倍），大量新进入企业初期的专利组合往往集中在少数热门 IPC 上（如 G06F、H04L）；
技术数基本停滞：2010 年 118 种 → 2015 年 130 种（仅 +10%），但企业数翻倍 → 每项技术被更多企业同时持有，M 矩阵密度上升但结构变”平”；
专利同质化涌入：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. 应对建议

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

推荐策略：

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

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

Step 10：保存结果

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

# 股票代码转 6 位字符串，不足前补 0
results[, 股票代码 := str_pad(as.character(股票代码),
                               width = 6, side = "left", pad = "0")]

# 保存 DTA
output_dta <- paste0(min_year, "~", max_year, "年上市公司数字技术复杂度_iter", n_iter_k, ".dta")
results_dt <- as_tibble(results[, .(股票代码, 年份, complexity, complexity_minmax,
                                     complexity_rank, complexity_global, complexity_panel)])
attr(results_dt, "label") <- "数据处理：微信公众号 RStata"
haven::write_dta(results_dt, output_dta)
cat("  DTA 已保存至:", output_dta, "\n")

# 保存 CSV（用于检查）
output_csv <- paste0(min_year, "~", max_year, "年上市公司数字技术复杂度_iter", n_iter_k, ".csv")
fwrite(results[, .(股票代码, 年份, complexity, complexity_minmax,
                   complexity_rank, complexity_global, complexity_panel)],
        output_csv, quote = TRUE)
cat("  CSV 已保存至:", output_csv, "\n")

关键注意事项

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

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

指标背景

计算步骤

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

核心公式

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

数据准备

样本筛选说明

详细计算代码

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

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

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

Step 3：读取 CSV 并筛选数据

Step 4：专利去重

Step 5：提取 IPC 四位码

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

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

Step 7：合并结果

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

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

方式 2：逐年百分位排名

方式 3：全局 Min-Max

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

代码实现

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

Step 9：描述性统计

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

现象描述

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

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

数据层面的深层原因

解决方案

对后续计量回归的影响

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

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

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

4. 应对建议

Step 10：保存结果

关键注意事项

相关文章

评论