名师讲堂｜使用 Stata 测算专利的创新突破度 CD 指数：以上市公司专利数据为例（更新）

3374 字 约 12 分钟 名师讲堂 Stata专利上市公司创新与知识

上次课程的代码里面计算的时候忽略了”如果某个专利 j 没有引用专利 i 但是引用了其后项引用，则 cd = 0” 的情况，导致计算结果里面 1 和 -1 出现的太多。所以赶紧更正了下。同时本次使用 Mata 重写了核心计算函数，大幅提升了运算速度。

在《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 的该项引用的 CD 指数为 -1；
  2. 若后面的专利只引用专利 i 但未引用专利 i 的后向引用专利，那么专利 i 的该项引用的 CD 指数为 1；
  3. 若后面的专利只引用专利 i 的后向引用专利，但是没有引用专利 i，那么专利 i 的该项引用的 CD 指数为 0；

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

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

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

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

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

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

由于全部的数据非常大，所以附件中提供了已经由 R 语言版本课程处理好的三个文件，可以直接使用：

  • patent_citations.dta：全部专利间的引用关系，newipzlid 引用 newipzlid2；
  • patent_citations2.dta：在 patent_citations 基础上，只保留引用间隔在 5 年以内且 date - date2 >= 0 的记录；
  • 具有引用与被引用信息的专利unique_newipzlid.dta：patent_citations 数据中所有出现过的专利编号。

这三个文件并不是样本，而是可以计算所有专利 CD 指数的全部数据。

如果需要自己从 CSV 文件转换为 dta 格式，可以参考如下代码（附件中已提供 dta 文件，这部分可以跳过）：

如果电脑内存较小，可以参考这个课程读取超大文件：Stata 如何读取超大的 dta 和 csv 文件：https://rstata.duanshu.com/#/brief/course/fbf32233112146fdad3f9788a6f8881f

*- csv -> dta
import delimited using "patent_citations.csv", clear

gen tempdate = date(date, "YMD")
format tempdate %tdCY-N-D
drop date
ren tempdate date

gen tempdate2 = date(date2, "YMD")
format tempdate2 %tdCY-N-D
drop date2
ren tempdate date2
save patent_citations, replace

import delimited using "patent_citations2.csv", clear

gen tempdate = date(date, "YMD")
format tempdate %tdCY-N-D
drop date
ren tempdate date

gen tempdate2 = date(date2, "YMD")
format tempdate2 %tdCY-N-D
drop date2
ren tempdate date2
save patent_citations2, replace

import delimited using "具有引用与被引用信息的专利unique_newipzlid.csv", clear
save "具有引用与被引用信息的专利unique_newipzlid.dta", replace

附件中提供的 patent_citations2.dta 文件就是根据上述代码处理得到的，保留了引用间隔 5 年内的专利引用记录。

单专利计算示例（Stata 基础版）

由于需要考虑这种情况：”若后面的专利只引用专利 i 的后向引用专利，但是没有引用专利 i，那么专利 i 的该项引用的 CD 指数为 0”，所以必须循环所有的专利来计算。

以某个专利为例，演示完整的计算逻辑：

*- 以某个专利为例：
local id = 20190809613
use patent_citations2, clear
keep if newipzlid2 == `id'
local date2 = date2[1]
save testpat, replace

*- 所有可能的 j
*- 准备 i 的后向引用
use patent_citations2, clear
keep if newipzlid == `id'
keep newipzlid2
if _N == 0 {
  cap drop newipzlid2
  set obs 1
  gen newipzlid2 = 0
}
cap duplicates drop
save backrefs, replace

*- 申请日期晚于 i
use patent_citations2, clear
keep if date >= `date2'
save part0, replace

*- 引用 i 的
use part0, clear
keep if newipzlid2 == `id'
keep newipzlid
if _N == 0 {
  cap drop newipzlid
  set obs 1
  gen newipzlid = 1
}
cap duplicates drop
save part1, replace

*- 引用 i 后向引用的
use part0, clear
merge m:1 newipzlid2 using backrefs
keep if _m == 3
keep newipzlid
if _N == 0 {
  cap drop newipzlid
  set obs 1
  gen newipzlid = 2
}
cap duplicates drop
save part2, replace

*- 由此就可以分出三种情况：
*- 此处代码需要下载讲义材料查看～

这样就得到了一个专利的 CD 指数。循环所有的即可。

但是这种逐步读写临时文件的方式效率极低，不适合大规模计算。

使用 Mata 大幅提速

显然单专利的 dta 文件流式计算非常耗时。Stata 内置的 Mata 语言可以直接操作内存矩阵，配合 asarray（哈希表），速度比循环合并快几十倍。

asarray 简介

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

*- 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 版）

use patent_citations2, clear

mata:

data = st_data(., ("newipzlid", "newipzlid2", "date", "date2"))
N = rows(data)

// 选择某个专利尝试
id = 20190809613

// 该专利的申请日期
idx = selectindex(data[., 2] :== id)
date_i = data[idx[1], 4]

date_i

// 专利 i 的后向引用
idx_back = selectindex(data[., 1] :== id)
if (rows(idx_back) > 0) {
    backrefs = uniqrows(data[idx_back, 2])
} else {
    backrefs = J(0, 1, .)
}

backrefs

// 所有可能的 j
idx_part0 = selectindex(data[., 3] :>= date_i)
part0 = data[idx_part0, .]

part0[1,.]
npart0 = rows(part0)
npart0

// part1：引用 i 的
idx_p1 = selectindex(part0[., 2] :== id)
if (rows(idx_p1) > 0) {
    part1 = uniqrows(part0[idx_p1, 1])
} else {
    part1 = J(0, 1, .)
}
part1

// part2：引用 i 后向引用的
// 构造 backrefs 的哈希表
Aback = asarray_create("real")
for (k = 1; k <= rows(backrefs); k++) asarray(Aback, backrefs[k,1], 1)

// 循环查找 part0 里面的被引专利是否在 backrefs 里面：
// part0[r,2] 表示专利 j 的引用专利
idx_p2 = J(npart0, 1, .)
cnt_p2 = 0
for (r = 1; r <= npart0; r++) {
    if (asarray_contains(Aback, part0[r,2])) {
        cnt_p2++
        idx_p2[cnt_p2] = r
    }
}
if (cnt_p2 > 0) {
    part2 = uniqrows(part0[idx_p2[1..cnt_p2], 1])
} else {
    part2 = J(0, 1, .)
}

part2

// 比较 part1 和 part2 的关系
// 构造 part2 的哈希表
A2 = asarray_create("real")
for (k = 1; k <= rows(part2); k++) asarray(A2, part2[k,1], 1)

sum_cd = 0
cnt    = 0

// 如果 part1 的元素在 part2 里面 -> 引用 i，引用后向 -> -1
// 如果 part1 的元素不在 part2 里面 -> 引用 i，没有引用后向 -> 1
for (k = 1; k <= rows(part1); k++) {
    if (asarray_contains(A2, part1[k,1])) {
        sum_cd = sum_cd + (-1)
    } else {
        sum_cd = sum_cd + 1
    }
    cnt++
}

// 构造 part1 的哈希表
A1 = asarray_create("real")
for (k = 1; k <= rows(part1); k++) asarray(A1, part1[k,1], 1)

// 如果 part2 的元素不在 part1 的里面 -> 没有引用 i，引用后向 -> 0，数量 +1，总值不变
for (k = 1; k <= rows(part2); k++) {
    if (!asarray_contains(A1, part2[k,1])) cnt++
}

// 最后，CD 指数就是均值
mean_cd = (cnt > 0 ? sum_cd / cnt : .)

printf("cited_pat = %12.0g | cd = %9.6f (n = %g)\n", id, mean_cd, cnt)
end

这样就得到了单个专利的 CD 指数。接下来把它封装成一个可复用的函数。

封装 compute_cd 函数

*- 这样我们就可以编写一个函数计算某个专利的 CD 指数了：
mata:

real scalar compute_cd(real scalar id, real matrix data)
{
    real scalar date_i, npart0, cnt_p2, cnt_tmp, sum_cd, cnt, mean_cd, r, k
    real matrix idx, idx_back, backrefs, idx_part0, part0, idx_p1, part1
    real matrix idx_p2, part2, tmp
    transmorphic Aback, A1, A2

    /* 获取 id 的申请日期 */
    idx = selectindex(data[., 2] :== id)
    if (rows(idx) == 0) return(.)
    date_i = data[idx[1], 4]

    /* backrefs：id 引用的专利集合 */
    idx_back = selectindex(data[., 1] :== id)
    if (rows(idx_back) > 0) {
        backrefs = uniqrows(data[idx_back, 2])
    } else {
        backrefs = J(0, 1, .)
    }

    /* part0：申请日期 >= date_i 的行 */
    idx_part0 = selectindex(data[., 3] :>= date_i)
    part0  = data[idx_part0, .]
    npart0 = rows(part0)

    /* part1：part0 中引用了 id 的专利 */
    idx_p1 = selectindex(part0[., 2] :== id)
    if (rows(idx_p1) > 0) {
        part1 = uniqrows(part0[idx_p1, 1])
    } else {
        part1 = J(0, 1, .)
    }

    /* part2：part0 中引用了 backrefs 的专利 */
    Aback = asarray_create("real")
    for (k = 1; k <= rows(backrefs); k++) asarray(Aback, backrefs[k,1], 1)

    idx_p2 = J(npart0, 1, .)
    cnt_p2 = 0
    for (r = 1; r <= npart0; r++) {
        if (asarray_contains(Aback, part0[r,2])) {
            cnt_p2++
            idx_p2[cnt_p2] = r
        }
    }
    if (cnt_p2 > 0) {
        part2 = uniqrows(part0[idx_p2[1..cnt_p2], 1])
    } else {
        part2 = J(0, 1, .)
    }

    /* 计算 cd */
    A2 = asarray_create("real")
    for (k = 1; k <= rows(part2); k++) asarray(A2, part2[k,1], 1)

    sum_cd = 0
    cnt    = 0
    for (k = 1; k <= rows(part1); k++) {
        if (asarray_contains(A2, part1[k,1])) {
            sum_cd = sum_cd - 1
        } else {
            sum_cd = sum_cd + 1
        }
        cnt++
    }

    A1 = asarray_create("real")
    for (k = 1; k <= rows(part1); k++) asarray(A1, part1[k,1], 1)
    for (k = 1; k <= rows(part2); k++) {
        if (!asarray_contains(A1, part2[k,1])) cnt++
    }

    mean_cd = (cnt > 0 ? sum_cd / cnt : .)
    return(mean_cd)
}

end

上述函数代码已单独存放在附件的 compute_cd.do 中，后续需要使用时 do compute_cd.do 即可加载。

调用示例：

*- 调用示例
use patent_citations2, clear
mata:
data = st_data(., ("newipzlid", "newipzlid2", "date", "date2"))
id = 20190559229
cd = compute_cd(id, data)
printf("cited_pat = %21.0g | cd = %9.6f\n", id, cd)
end

计算多个专利：

*- 计算多个
mata:
// 选择 10 个
id = uniqrows(data[34560000..34560010, 2])
r = rows(id)
cd_mat = J(r, 2, .)
for (i = 1; i <= r; i++) {
  i
  cd_mat[i, 1] = id[i,1]
  cd_mat[i, 2] = compute_cd(id[i,1], data)
}
cd_mat
end

输出结果示例：

*>                  1              2
*>    +-------------------------------+
*>  1 |    2010799085    .5714285714  |
*>  2 |    2010799103              1  |
*>  3 |    2010799112             .2  |
*>  4 |    2010799129   -.1111111111  |
*>  5 |    2010799163    .0384615385  |
*>  6 |    2010799164    .2857142857  |
*>  7 |    2010799223    .0909090909  |
*>  8 |    2010799239    .0833333333  |
*>    +-------------------------------+

筛选无需计算的专利（CD 指数恒为 1）

与 R 语言版本思路完全相同：如果某个专利 i 在 patent_citations2 里面没有引用任何其他专利，则它没有后向引用，所有后续专利的引用方向都是 +1，CD 指数必然为 1，无需计算。

*- 存在后向引用的专利
use patent_citations2, clear
keep newipzlid
ren newipzlid newipzlid2
duplicates drop
save temp, replace

*- i 不存在后向引用的：cd_index 直接赋值为 1
use patent_citations2, clear
merge m:1 newipzlid2 using temp
keep if _m == 1
keep newipzlid2
duplicates drop
gen cd_index = 1
save nobackrefs, replace

*- 存在后向引用的，才需要计算
use patent_citations2, clear
merge m:1 newipzlid2 using temp
keep if _m == 3
keep newipzlid2
duplicates drop
ren newipzlid2 newipzlid
save withbackrefs, replace

附件中也提供了 nobackrefs.dta 文件。

然后读取上市公司专利信息，筛选出真正需要计算的专利：

*- 读取上市公司的
use 2010年上市公司与专利数据匹配结果, clear

*- 去除不需要计算的
merge m:1 newipzlid using withbackrefs
keep if _m == 3
keep newipzlid
duplicates drop
save 待计算, replace

单线程示例计算（5 个专利）

不过这里我们就不完全计算了，而是选择 5 个：

do compute_cd.do
use 待计算, clear
mata:
// 选择 5 个
idall = st_data(., "newipzlid")
id = idall[1..5, .]
r = rows(id)
cd_mat = J(r, 2, .)
for (i = 1; i <= r; i++) {
  i
  cd_mat[i, 1] = id[i,1]
  cd_mat[i, 2] = compute_cd(id[i,1], data)
}
cd_mat
end

输出结果示例：

*>                 1             2
*>    +-----------------------------+
*>  1 |   2010000002            .2  |
*>  2 |   2010000004   .3333333333  |
*>  3 |   2010000005   .8620689655  |
*>  4 |   2010000006   .6666666667  |
*>  5 |   2010000091           .25  |
*>    +-----------------------------+

通过这个结果我们就可以估算计算全部需要多长时间了。

如果要计算较多专利，可以每次把结果写入 CSV 文件：

cap erase cd_results.csv
mata:
fh = fopen("cd_results.csv", "w")
fput(fh, "cited_pat,cd")
// 选择 5 个
idall = st_data(., "newipzlid")
id = idall[1..5, .]
id
r = rows(id)
for (i = 1; i <= r; i++) {
    i
    idx = id[i, 1]
    cd = compute_cd(idx, data)

    line = sprintf("%12.0f,%9.6f", idx, cd)
    fput(fh, line)
}

fclose(fh)
printf("完成，结果已保存至 cd_results.csv\n")
end

import delimited using "cd_results.csv", clear

多线程并行计算（parallel 命令）

显然单线程计算是个极其耗时的工作。Stata 可以使用 parallel 包进行多线程并行计算：

*- 最后就是多线程运算了：
*- 此处代码需要下载讲义材料查看～

其中 para_matacode.do 是每个 worker 执行的子脚本，附件中已提供。其核心逻辑为每个 worker 根据自身编号取一列数据，循环调用 compute_cd 函数并将结果写入各自的 CSV/dta 文件：

do compute_cd.do
use 待计算2, clear
mata:
worker_id = st_global("pll_instance")
worker_id
filename = "csv/para_res_" + worker_id + ".csv"
fh = fopen(filename, "rw")
z = strtoreal(worker_id)
id = st_data(1::5, z)
r = rows(id)
for (i = 1; i <= r; i++) {
    idx = id[i, 1]
    cd = compute_cd(idx, data)
    fput(fh, sprintf("%12.0f,%9.6f", idx, cd))
}
fclose(fh)
stata("import delimited using " + filename + ", clear varnames(nonames)")
stata("save dta/para_res_" + worker_id + ".dta, replace")
end

电脑性能不好的小伙伴还是建议老老实实单线程。

合并计算结果（使用附件中提供的 appendall.ado）：

*- 合并结果
appendall dta, clear
drop filename
ren v1 newipzlid
ren v2 cd_index
save cd_index, replace

与上市公司数据关联

然后我们读取和处理上市公司专利，去重的时候要尽可能保留出现在 patent_citations 里面的专利：

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

和上述 CD 指数计算结果关联：

merge m:1 newipzlid using cd_index
replace cd_index = 0 if mi(cd_index)
collapse (mean) cd_index, by(股票代码 年份)

由于上面并没有计算所有的专利 CD 指数，这里直接用 0 代替缺失的其实是不对的，但是如果你计算了全部可计算的专利，这样做是可以的。没有被引用的专利 CD 指数是 0。

分公司汇总就可以得到公司层面的结果了。

仔细学习的小伙伴也会发现，其实这个代码是可以计算全部专利 CD 指数的，不过由于运算过于耗时，还是只计算自己需要的就好。

点击这里跳转到 RStata 短书平台获取附件：名师讲堂｜使用 Stata 测算专利的创新突破度 CD 指数：以上市公司专利数据为例（更新）

相关文章

评论