名师讲堂｜使用 Stata 测算专利的创新突破度：ACD指数和mACD指数

讲义材料已更新！

本课程介绍如何使用 Stata 配合 Mata 语言测算专利的 ACD 指数和 mACD 指数，并提供完整的代码模板。

本次课程在原始 CD 指数（Funk & Owen-Smith, 2017）的基础上，介绍徐照宜等（2023）在《金融研究》上发表的改进算法：ACD 指数和 mACD 指数。改进的核心是解决”从未被引用专利”的歧义问题，以及调整分母以更好地反映专利的突破性。

本文提供完整的 Stata + Mata 实现：用 Mata 在内存矩阵上直接运算（配合 asarray 哈希表加速），并用 Stata 官方生态的 parallel 命令做多核并行。完整可直接运行的代码见附件 main.do。

一、背景：原始 CD 指数及其问题

1.1 CD 指数公式

1.2 三种引用情况的 CD 值

类型	条件	CD 值	含义
圆形	引用目标专利，不引用后向引用	+1	突破性
正方形	既引用目标专利又引用后向引用	-1	渐进性
三角形	只引用后向引用专利	0	中性

论文中的示意图如下：

1.3 原始 CD 指数的问题

1. 从未被引用的专利 CD = 0：这与”渐进性专利”的 CD = -1 混淆（沉默专利 vs 组合型专利无法区分）。
2. 正负值相抵：多个专利的 CD 值正负相抵，汇总时相互抵消，难以加总衡量整体突破性。

二、改进的 ACD 和 mACD 指数

三、数据来源与处理过程

3.1 原始数据

本课程使用的数据来源于中国专利数据库，包含 1985-2024年 上市公司申请的专利及其引用关系数据。前面的处理（去重、5 年窗口筛选、识别无后向引用专利）已在 R 语言版本课程中完成，并导出为 .dta 格式供 Stata 直接使用。

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

1985~2024 年全部专利引用与被引用详细信息: https://rstata.duanshu.com/#/brief/course/225ad0b59a9945e1831d2e8b96ca1001

3.2 数据文件说明

文件名	行数	列名	说明
具有引用与被引用信息的专利unique_newipzlid.dta	23,330,294	newipzlid	所有具有引用与被引用信息的专利ID列表
patent_citations.dta	63,020,845	newipzlid, newipzlid2, date, date2	原始专利引用关系数据
nobackrefs.dta	8,277,416	newipzlid2, cd_index	没有后向引用的专利（cd_index = 1）
patent_citations2.dta	46,102,368	newipzlid, newipzlid2, date, date2	5年窗口筛选后的引用数据（推荐使用）
2010年上市公司与专利数据匹配结果.dta	81,332（含重复）	股票代码, 年份, newipzlid, 公开公告号, 申请号	上市公司专利数据（含重复，需按 R 脚本口径去重，见 §3.4）

变量说明：

• newipzlid：引用专利的ID（申请日的专利，引用方）
• newipzlid2：被引用专利的ID（被引用的专利，被引用方）
• date：引用专利的申请日期
• date2：被引用专利的申请日期

3.3 数据处理流程

数据从原始专利库到最终可用格式，需要经过以下处理步骤（结果已随课程提供，无需重跑）：

┌─────────────────────────────────────────────────────────────────┐
│  原始数据：全部专利引用与被引用信息分年.dta                       │
│  （包含专利引用与被引用信息，含省市区信息）                        │
└─────────────────────────────────────────────────────────────────┘
                              ↓
                    ① 根据公开公告号去重
                              ↓
┌─────────────────────────────────────────────────────────────────┐
│  patent_citations.dta                                            │
│  - 结构：newipzlid（引用方）、newipzlid2（被引用方）、           │
│          date（引用方申请日）、date2（被引用方申请日）            │
└─────────────────────────────────────────────────────────────────┘
                              ↓
              ② 保留5年内的相互引用（时间窗口筛选）
              条件：0 ≤ date - date2 ≤ 5年
                              ↓
┌─────────────────────────────────────────────────────────────────┐
│  patent_citations2.dta  （推荐使用的主数据）                     │
│  筛选理由：专利通常在申请后5年内被引用                             │
└─────────────────────────────────────────────────────────────────┘
                              ↓
         ③ 识别无后向引用的专利
             即：被其他专利引用，但自己不引用任何专利的专利
             特点：这些专利的 ACD 直接设为 2（纯突破性）
                              ↓
┌─────────────────────────────────────────────────────────────────┐
│  nobackrefs.dta                                                  │
│  newipzlid2（无后向引用的专利ID），cd_index = 1                  │
└─────────────────────────────────────────────────────────────────┘

3.4 上市公司专利去重

2010年上市公司与专利数据匹配结果.dta 中的上市公司专利存在重复记录（同一专利以相同的 公开公告号 或 申请号 多次出现）。在计算 ACD/mACD 之前必须去重，否则重复专利会被重复计入、扭曲突破度结果。

去重共三步：

1. 清洗公开公告号：去掉末尾字母（如 CN12345A → CN12345），对应 R 的 str_replace(公开公告号, “[A-Z]$”, “”)。
2. 按 股票代码 + 年份 + 公开公告号_clean 去重。
3. 再按 股票代码 + 年份 + 申请号 去重。

优先保留引用网络专利：为避免误删处于引用网络中的专利，main.do 在去重前先用 具有引用与被引用信息的专利unique_newipzlid.dta（2,330 万余条）给每个公司专利打标记 keep_cite；每一步 duplicates drop 前先按 keep_cite 降序排序，使”有引用信息的那一条”排在分组首位而被保留。这样可做到”尽可能保留”引用网络专利——只要某个重复组里存在引用网络专利，被保留的就是它，而非无引用信息的重复项。

对应 main.do 第二节代码：

*- 2.1 上市公司专利去重（参考 去除重复专利.R，并优先保留有引用信息的专利）
use "2010年上市公司与专利数据匹配结果.dta", clear
gen double nid = newipzlid
drop newipzlid
rename nid newipzlid
merge m:1 newipzlid using "具有引用与被引用信息的专利unique_newipzlid.dta", keep(match master)
gen byte keep_cite = (_merge == 3)
drop _merge
gen str13 公开公告号_clean = ustrregexra(公开公告号, "[A-Z]$", "")
gsort 股票代码 年份 公开公告号_clean -keep_cite
duplicates drop 股票代码 年份 公开公告号_clean, force
gsort 股票代码 年份 申请号 -keep_cite
duplicates drop 股票代码 年份 申请号, force
gsort 股票代码 年份 newipzlid -keep_cite
duplicates drop 股票代码 年份 newipzlid, force
keep 股票代码 年份 newipzlid
save company_patents_unique.dta, replace

3.5 变量说明（引用数据）

• newipzlid：引用专利的ID（引用方）
• newipzlid2：被引用专利的ID（被引用方）
• date：引用专利的申请日期
• date2：被引用专利的申请日期

3.6 数据读取（Stata）

读取本课程提供的处理后数据（推荐使用 patent_citations2.dta）：

*- 读取专利引用数据（5年窗口，推荐使用）
use patent_citations2.dta, clear
describe
list in 1/10

* 基本规模统计
di "专利引用数据行数: " _N

更直接的规模统计（避免嵌套宏）也可以写成：

use patent_citations2.dta, clear
di "专利引用数据行数: " _N
use nobackrefs.dta, clear
di "无后向引用专利数: " _N

专利引用数据的变量说明：

• newipzlid：引用专利的ID（引用方）
• newipzlid2：被引用专利的ID（被引用方）
• date：引用专利的申请日期
• date2：被引用专利的申请日期

四、计算示例（Stata + Mata）

4.1 单个专利的计算逻辑

以某个专利 id = 20190809613 为例，详细展示 ACD 和 mACD 的计算过程。下面先用基础 Stata 命令把每一步的中间结果写出来，便于理解；生产环境请直接用 §4.2 的 Mata 函数。本小节完整代码见附件 单个专利测试.do。

*- 以某个专利为例：
*- 此处代码需下载讲义材料查看～

该专利的运行结果如下（parta=0、partb=1、partc=0，即只有一条”只引用目标专利”的记录，属纯突破性）：

*>     +-----------------------------------------------------------+
*>     | cited_~t   parta   partb   partc   n_t   m_t   ACD   mACD |
*>     |-----------------------------------------------------------|
*>  1. | 2.02e+10       0       1       0     1     1     2      2 |
*>     +-----------------------------------------------------------+

上述步骤对应论文中的三种情况：

• parta：既引用目标专利又引用后向引用（cd = -1）
• partb：只引用目标专利（cd = +1）
• partc：只引用后向引用（cd = 0，不计入求和）

4.2 使用 Mata 函数批量封装

把上述逻辑封装成 Mata 函数 compute_acd()，直接操作内存矩阵（配合 asarray 哈希表做 O(1) 集合查找），速度提升数十倍。该函数已内置在 main.do 与附件 compute_acd_mACD.do 中：

关于 transmorphic：Mata 是静态强类型语言，每个变量都要声明元素类型（如 real matrix）。transmorphic 是其中的通配类型（万能类型，类似 C 的 void*、TypeScript 的 any），可承载任意类型的对象。这里 Aback、A1、A2 存放的是 asarray_create() 返回的哈希表句柄——按 Stata 规定，关联数组只能用 transmorphic 标量承载，写成 real scalar 会报错。

关于 asarray 简介：

asarray 是 Mata 里的哈希表（字典），核心作用：用「键 key」快速查「值 value」，速度极快（O(1)），比循环、匹配快几十倍。

*- 1. 创建哈希表
h = asarray_create("real", 1)     // 键是数字（最常用）
h = asarray_create("string", 1)   // 键是字符串

*- 2. 放入键值对
asarray(h, key, value)

*- 3. 判断某个 key 是否存在，返回 1=存在，0=不存在
asarray_contains(h, key)

*- 4. 读取键对应的值
val = asarray(h, key)

*- mata 代码中 transmorphic 表示声明通配型对象
*- 此处代码需下载讲义材料查看～

调用示例：

use patent_citations2.dta, clear
mata:
data = st_data(., ("newipzlid", "newipzlid2", "date", "date2"))
end

*- 加载函数
do compute_acd_mACD.do

*- 测试计算
mata:
id = 20190809613
result = compute_acd(id, data)
printf("ACD = %.6f, n_t = %.0f, m_t = %.0f, mACD = %.6f\n",
       result[1], result[2], result[3], result[4])
end

五、批量计算与多线程优化

5.1 筛选需要计算的专利

并非所有专利都需要逐一计算：

• 无后向引用的专利：所有引用它的专利都落在 partb，ACD = 2。Mata 函数在计算时会自动得到这一结果，无需预先拆分。
• 从未被引用的专利：compute_acd() 返回缺失值，主程序将其 ACD 设为 0。

因此 main.do 直接把”被引用过的公司专利”全部交给 Mata 计算，逻辑非常简洁：

*- 2.3 公司专利中、且至少被引用过一次的专利 = 需要计算的清单
use company_patents_unique.dta, clear
rename newipzlid cited_id
merge m:1 cited_id using cited_patents.dta
keep if _merge == 3
drop _merge
rename cited_id newipzlid
duplicates drop newipzlid, force
save patents_to_calc.dta, replace

若想单独列出”无后向引用的公司专利”（ACD 恒为 2）用于核查，可用稳健的二分法：
merge m:1 newipzlid using with_backrefs.dta, keep(match master) nogen 后 keep if _merge == 1。

5.2 单线程批量计算（默认路径）

main.do 的默认路径是单线程批量计算：把引用数据整体载入 Mata 一次，然后逐个专利调用 compute_acd()，把结果写入矩阵 R，最后用 st_store 一次性写回 Stata 数据集（避免 Mata sprintf 对大整数 ID 的格式化误差）。

*- 3.1 将引用数据整体载入 Mata（仅一次）
use patent_citations2.dta, clear
*- 此处代码需下载讲义材料查看～

format newipzlid %18.0f
save cd_results.dta, replace

为什么用 st_store 而不是 sprintf 拼字符串？本机 Mata 版本对带宽度/精度的格式（如 %15.0f、%.6f）支持有限，直接把结果矩阵写回 Stata、由 Stata 负责数值格式化最稳妥。

5.3 多线程计算（Stata 的 parallel 方案）

由于计算量较大，可使用 Stata 的 parallel 命令做多核并行。核心思路：把待计算专利拆成 K 份存入 dta1/，启动 K 个 Stata 实例，每个实例运行子程序 para_matacode.do 计算自己负责的那一份、结果写入 dta2/，最后用 appendall 合并。

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

要点说明

• para_matacode.do 内部内联了 compute_acd() 函数（拆成两个 mata: 块：函数定义块 + worker 执行块），并在每个 worker 中重新载入 patent_citations2.dta。请确保内存充足（建议 K × 1.5GB 以上空闲内存）。
• worker 读取 dta1/待计算_part<实例号>.dta，用 st_store 把结果矩阵直接写回 Stata 数据集，保存为 dta2/para_res_<实例号>.dta（已带变量名 newipzlid/acd/nt/mt/macd），因此合并时用 appendall dta2, clear 即可，无需再 rename。
• sample 0.1 仅用于快速试跑（抽取 0.1% 专利），正式全量计算时删除该行。子程序 para_matacode.do 中避免使用 /* */ 块注释——parallel 会把 dofile 包在 noisily { } 中执行，块注释会被误解析报错，故一律用 * 行注释。

para_matacode.do 的 worker 主逻辑如下：

mata:
inst = strtoreal(st_global("pll_instance"))

stata("use patent_citations2.dta, clear")
data = st_data(., ("newipzlid", "newipzlid2", "date", "date2"))

stata("use dta1/待计算_part" + sprintf("%f", inst) + ".dta, clear")
idall = st_data(., "newipzlid")
r = rows(idall)

R = J(r, 5, .)
for (i = 1; i <= r; i++) {
    idx = idall[i, 1]
    res = compute_acd(idx, data)
    if (missing(res[1])) res = (0, 0, 0, 0)
    R[i, 1] = idx
    R[i, 2] = res[1]
    R[i, 3] = res[2]
    R[i, 4] = res[3]
    R[i, 5] = res[4]
}

stata("clear")
st_addobs(rows(R))
st_addvar("double", ("newipzlid", "acd", "nt", "mt", "macd"))
st_store(., ("newipzlid", "acd", "nt", "mt", "macd"), R)
stata("save dta2/para_res_" + sprintf("%f", inst) + ".dta, replace")
printf("worker %g 完成，计算 %g 个专利\n", inst, r)
end

六、结果汇总

6.1 合并计算结果

cd_results.dta 已包含全部”被引用过的公司专利”（无后向引用者自动为 ACD = 2）。将其与全部公司专利左连，从未被引用的公司专利 ACD = 0：

use company_patents_unique.dta, clear
merge m:1 newipzlid using cd_results.dta
replace acd  = 0 if _merge == 1
replace macd = 0 if _merge == 1
replace nt   = 0 if _merge == 1
replace mt   = 0 if _merge == 1
drop _merge
save all_acd_patent.dta, replace

6.2 公司层面汇总

将专利层面的 ACD / mACD 汇总到公司-年份层面（取均值）：

collapse (mean) acd macd nt mt, by(股票代码 年份)
save company_year_acd.dta, replace
list in 1/10

6.3 统计汇总

use company_year_acd.dta, clear
di "=== ACD 统计 ==="
summarize acd
di "=== mACD 统计 ==="
summarize macd

点击这里跳转到 RStata 短书平台获取附件：名师讲堂｜使用 Stata 测算专利的创新突破度：ACD指数和mACD指数