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

在《Nature》论文：「Papers and patents are becoming less disruptive over time」 中使用了 CD 指数数据，该指标的计算公式如下：

其中当专利 j 引用了专利 i 时 fit 取1，否则取 0；若专利 j 引用了专利 i 的后向引用专利时 bit 取 1，否则取 0；

附件中也提供了该论文的 pdf 文件。

论文中也提供了一个示意图帮助我们理解这个指数：

看起来很复杂，其实只需要关注两个：

1. 若后面的专利既引用专利 i 又引用专利 i 的后向引用专利，那么专利 i 的该项引用的 CD1 指数 -1；
2. 若后面的专利只引用专利 i 但未引用专利 i 的后向引用专利，那么专利i的该项引用的 CDI 指数为 1；
3. 其他的情况都是 0 ；

这样直接计算得到的是每个专利引用关系的 CD 指数，平均之后就得到了某个专利的 CD 指数。

今天的课程中我们将以上市公司专利为例进行讲解。

在计算之前我们需要准备两组数据：

1. 一个是专利的引用关系，也就是每个专利引用了其他的哪些专利；
2. 第二个是专利的属性信息，也就是每个公司每年申请了哪些专利。

这两个数据我们都分享过。全部专利的引用与被引用关系数据在这里：

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

由于全部的数据非常大，所以附件中仅仅提供了 2010 年前的数据作为样本，通过类似如下代码即可处理得到：

这部分代码就不用运行了，“全部专利引用与被引用信息分年.rds” 是一个超过 40 GB 的大数据，代码需要电脑性能足够好才能运行。

library(tidyverse)
read_rds("全部专利引用与被引用信息分年.rds") -> df1

# 根据公开公告号去重
df1 %>%
  mutate(公开公告号_original = str_remove(公开公告号_original, "[A-Z]$"),
         公开公告号 = str_remove(公开公告号, "[A-Z]$")) %>%
  distinct(公开公告号_original, 公开公告号, .keep_all = T) -> df2

# 引用信息
bind_rows(
  df2 %>%
    filter(引用或被引用 %in% c("他引信息", "自引信息")) %>%
    set_names(c("newipzlid1", "申请日1", "公开公告号1",
                "省1", "省代码1", "市1", "市代码1", "县1", "县代码1",
                "引用或被引用",
                "newipzlid2", "申请日2", "公开公告号2",
                "省2", "省代码2", "市2", "市代码2", "县2", "县代码2")) %>%
    select(-引用或被引用),
  df2 %>%
    filter(引用或被引用 %in% c("被他引信息", "被自引信息")) %>%
    set_names(c("newipzlid2", "申请日2", "公开公告号2",
                "省2", "省代码2", "市2", "市代码2", "县2", "县代码2",
                "引用或被引用",
                "newipzlid1", "申请日1", "公开公告号1",
                "省1", "省代码1", "市1", "市代码1", "县1", "县代码1")) %>%
    select(-引用或被引用)
) %>%
  distinct() -> df0a

df0a

# 保存
df0a %>%
  write_rds("全部专利引用信息（去重后）.rds")

# 选择 2010 年前的样本
read_rds("全部专利引用信息（去重后）.rds") %>%
  filter(year(申请日1) <= 2010) %>%
  write_csv("2010年及之前专利引用信息.csv")

专利的属性信息就是上市公司专利匹配结果：

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

同样这个数据也是非常巨大，附件中提供了 “2010年上市公司与专利数据匹配结果.csv” 样本。

需要注意，如果想要计算 2010 年的专利 CD 指数，需要使用全部的 2010 年及之前的专利引用信息。因为这里需要专利引用的引用信息，2010 年的专利会引用 2010 年之前的专利，然后这些专利还会引用更早期的专利。

下面我们开始使用 Stata 处理，首先准备数据：

cd "~/Desktop/使用 Stata 测算专利的创新突破度：CD指数/"
*- 导入 2010 年及之前专利引用信息
import delimited "2010年及之前专利引用信息.csv", clear
rename (newipzlid1 newipzlid2) (newipzlid newipzlid2)
keep newipzlid newipzlid2
save patent_citations.dta, replace

*- 获取所有不重复的专利编号
use patent_citations.dta, clear
keep newipzlid
rename newipzlid value
save temp1.dta, replace

use patent_citations.dta, clear
keep newipzlid2
rename newipzlid2 value
append using temp1.dta
duplicates drop
rename value newipzlid
save unique_newipzlid.dta, replace

*- 导入上市公司专利数据
import delimited "2010年上市公司与专利数据匹配结果.csv", clear
rename (股票代码 年份 公开公告号 申请号) (firm_id year 公开公告号 申请号)
merge m:1 newipzlid using unique_newipzlid.dta
drop if _m == 2

*- 生成标识变量
gen value = 1 if _m == 3
replace value = 0 if missing(value)

*- 清理公开公告号（去除末尾字母）
gen 公开公告号_clean = regexr(公开公告号, "[A-Z]$", "")

*- 按公开公告号去重，但是尽可能保存 _m == 3 的
bysort firm_id year 公开公告号_clean: egen max_value1 = max(value)
keep if value == max_value1
gsort firm_id year 公开公告号_clean
duplicates drop firm_id year 公开公告号_clean, force
drop 公开公告号_clean

*- 按申请号去重，但是尽可能保存 _m == 3 的
bysort firm_id year 申请号: egen max_value2 = max(value)
keep if value == max_value2
gsort firm_id year 申请号
duplicates drop firm_id year 申请号, force
drop 申请号

*- 清理并保存专利信息
keep firm_id year newipzlid
rename newipzlid patent_id
drop if mi(firm_id)
save patent_info.dta, replace

准备后向引用关系，也就是每个专利引用了哪些专利：

use patent_citations.dta, clear
rename (newipzlid newipzlid2) (cited_patent backward_cited)
save backward_refs.dta, replace

计算每次引用的 CD1 值：

*- 当专利 j 引用了专利 i 时 fit 取1，否则取 0；
*- 若专利 j 引用了专利 i 的后向引用专利时 bit 取 1，否则取 0；
*- 若后面的专利既引用专利 i 又引用专利 i 的后向引用专利，那么专利 i 的该项引用的 CD1 指数 -1
*- 若后面的专利只引用专利 i 但未引用专利 i 的后向引用专利，那么专利i的该项引用的 CDI 指数为 1
use patent_citations.dta, clear
rename (newipzlid newipzlid2) (citing_patent check_cited)
save temp2, replace

rename (citing_patent check_cited) (citing_patent cited_patent)
*- 合并被引用专利的后向引用专利
joinby cited_patent using backward_refs.dta, unmatched(master)
drop _m

*- 检查是否也引用了后向引用专利
joinby citing_patent using temp2, unmatched(master)
drop _m
save temp3, replace

*- 这里运行下面这句就明白上面的代码在做什么准备了：
*- keep if backward_cited == check_cited
*- 标记是否同时引用了后向引用专利
gen cites_backward = 0
replace cites_backward = 1 if backward_cited == check_cited & ///
    backward_cited != cited_patent & ///
    !missing(backward_cited)

*- 计算 CD1
bysort citing_patent cited_patent: egen any_backward = max(cites_backward)
gen cd1 = 1
replace cd1 = -1 if any_backward == 1

*- 保存 CD1 计算结果
keep citing_patent cited_patent cd1
duplicates drop
save cd1_values.dta, replace

cd1_values 是专利引用对的结果，分组汇总就得到了专利层面的结果：

use cd1_values.dta, clear
bysort cited_patent: egen cd2 = mean(cd1)
bysort cited_patent: gen n_citations = _N
keep cited_patent cd2 n_citations
duplicates drop
save cd2_values.dta, replace

再分企业平均就可以计算企业年度的 CD 指数：

use cd2_values.dta, clear
keep if cd2 == -1

*- 计算企业年度的 CD 指数
use cd2_values, clear
ren cited_patent patent_id
joinby patent_id using patent_info.dta

collapse (mean) CD = cd2 (count) n_patents = patent_id, by(firm_id year)

*- 导出结果
export excel using "2010年上市公司CD指数计算结果.xlsx", firstrow(variables) replace

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