名师讲堂｜使用 Stata 计算地区间技术互补指数

在附件文献「地区间技术互补与创新驱动型经济增长_郑江淮」中提到了地区间技术互补指数的指标：

今天的课程中我们将一起使用 Stata 复现这个指标。

使用的专利数据就是之前给大家分享的这个：

1985～2024 年专利申请与授权数据（版本 3，含申请人所处的省市区县）：https://rstata.duanshu.com/#/brief/course/2397451274c546d3a36e156ffc865988

在附件中我提供了 2001 年的专利数据作为示例。

首先处理专利数据，主要是去除重复专利、筛选发明专利、拆分 IPC 等：

clear all
mata: mata clear
*- 读取数据
use 2001, clear
*> (数据处理：微信公众号 RStata)

由于专利数据中同时包含了专利的申请和授权公告，所以有重复的，在统计数量前需要进行去重。通常公开公告号的结尾字母用以区分公告类别，例如 A 代表发明专利的申请公开，B 代表发明专利的授权公告，U 代表实用新型专利的授权公告，S 代表外观设计专利的授权公告。因此可以通过去除公开公告号结尾的字母进行去重，另外申请号也存在重复的，也需要进行去重。

*- 去除重复专利
gen 公开公告号_clean = ustrregexra(公开公告号, "[A-Z]$", "")
duplicates drop 公开公告号_clean, force
*> Duplicates in terms of 公开公告号_clean
*> (833 observations deleted)
duplicates drop 申请号, force
*> Duplicates in terms of 申请号
*> (13,567 observations deleted)
*- 保留发明专利和需要的变量
keep if index(专利类型, "发明")
*> (105,139 observations deleted)
ren 年份 year
ren IPC ipc
ren 市 city
drop if missing(city) | missing(ipc)
*> (6 observations deleted)
keep newipzlid year city ipc
save df1, replace
*> file df1.dta saved

IPC 分类号使用分号分隔，可以使用 split + gather 进行拆分和转置：

*- 拆分IPC并提取小类
use df1, clear
*> (数据处理：微信公众号 RStata)
replace ipc = ustrregexra(ipc, "[s.]", "")
*> (21,196 real changes made)
split ipc, p(";") gen(ipc_)
*> variables created as string:
*> ipc_1   ipc_3   ipc_5   ipc_7   ipc_9   ipc_11  ipc_13  ipc_15  ipc_17  ipc_19  ipc_21  ipc_23  ipc_25  ipc_27  ipc_29  ipc_31  ipc_33  ipc_35  ipc_37
*> ipc_2   ipc_4   ipc_6   ipc_8   ipc_10  ipc_12  ipc_14  ipc_16  ipc_18  ipc_20  ipc_22  ipc_24  ipc_26  ipc_28  ipc_30  ipc_32  ipc_34  ipc_36
drop ipc
gather ipc*
*> newipzlid year city
drop if mi(value)
*> (926,318 observations deleted)
drop var
gen ipc = substr(value, 1, 4)
drop value
save df2, replace
*> file df2.dta saved

gather 是外部命令，可以使用下面的代码安装：

ssc install tidy

下面就可以开始一步步的计算了。

步骤1：计算每个城市-年份-IPC的专利数量

use df2, clear
*> (数据处理：微信公众号 RStata)
contract year city ipc, freq(patent_count)
save patent_count, replace
*> file patent_count.dta saved

步骤2：计算每个城市每年所有 IPC 的专利总数

use patent_count, clear
*> (数据处理：微信公众号 RStata)
collapse (sum) total_city = patent_count, by(year city)
save city_year_total, replace
*> file city_year_total.dta saved

步骤3：计算每年每个IPC的全国专利总数

use patent_count, clear
*> (数据处理：微信公众号 RStata)
collapse (sum) total_ipc = patent_count, by(year ipc)
save ipc_year_total, replace
*> file ipc_year_total.dta saved

步骤4：计算 RPCA（显性专利比较优势指数）

use patent_count, clear
*> (数据处理：微信公众号 RStata)
merge m:1 year city using city_year_total, nogen
*>     Result                      Number of obs
*>     -----------------------------------------
*>     Not matched                             0
*>     Matched                            11,158
*>     -----------------------------------------
merge m:1 year ipc using ipc_year_total, nogen
*>     Result                      Number of obs
*>     -----------------------------------------
*>     Not matched                             0
*>     Matched                            11,158
*>     -----------------------------------------
gen share_city = patent_count / total_city
egen total_national = total(patent_count), by(year)
gen share_national = total_ipc / total_national
gen RPCA = share_city / share_national
gen x_id = cond(RPCA > 1, 1, 0)
save rpc_data, replace
*> file rpc_data.dta saved

步骤5：构建技术关联矩阵 Φ（基于所有年份的共现）

注意这里是需要把所有年份的发明专利数据放在一起处理，不过这里只提供了 2001 年的样本，所以还是以 2001 年的数据为例进行演示。

clear all
use df2, clear
*> (数据处理：微信公众号 RStata)
keep newipzlid ipc
duplicates drop _all, force
*> Duplicates in terms of newipzlid ipc
*> (41,568 observations deleted)
*- 删除只有一个 IPC 小类的，这种不存在共现问题
bysort newipzlid: gen ipc_count = _N
keep if ipc_count > 1
*> (15,144 observations deleted)
drop ipc_count

在 Mata 中生成每个专利的 IPC 组合：

mata:
    patent = st_sdata(., "newipzlid")
    ipc = st_sdata(., "ipc")

    uniq_patent = uniqrows(patent)
    printf("总专利数量：%f\n", rows(uniq_patent))

    // 以某个专利为例
    patent1 = "2001000007"
    idx1 = selectindex(patent :== patent1)
    ipclist = ipc[idx1,]
    ipclist

    // 生成组合
    fullcombinations = J(0, 3, "")
    for (r = 1; r <= rows(ipclist); r++) {
        for (c = 1; c <= rows(ipclist); c++) {
            fullcombinations = fullcombinations \ (patent1, ipclist[r,], ipclist[c,])
        }
    }
    fullcombinations

    ....
*- 这里的代码需要下载讲义材料查看
end

这段代码又复杂，效率又低，所以我还是觉得应该放弃 Mata，这部分操作可以调用一些 R 语言的代码快速实现。

关于如何在 Stata 中调用 R 语言代码，可以预先学习下平台上的课程「名师讲堂｜Stata 中文文本分析」的第一次课。

名师讲堂｜Stata 中文文本分析: https://rstata.duanshu.com/#/brief/course/b6a9efd94e5a48c2bba52dc9fdfd4291

安装所需的 R 包：

rcall vanilla: install.packages("tidyverse", repos = "https://mirrors.ustc.edu.cn/CRAN/", dependencies = T)
rcall vanilla: install.packages("haven", repos = "https://mirrors.ustc.edu.cn/CRAN/", dependencies = T)

使用 Stata 去除只有一个 IPC 小类的：

clear all
use df2, clear
*> (数据处理：微信公众号 RStata)
keep newipzlid ipc
duplicates drop _all, force
*> Duplicates in terms of newipzlid ipc
*> (41,568 observations deleted)
bysort newipzlid: gen ipc_count = _N
keep if ipc_count > 1
*> (15,144 observations deleted)
drop ipc_count
save patent_ipc, replace
*> file patent_ipc.dta saved

然后调用少量的 R 语言代码实现这种操作：

rcall vanilla: ///
    ....
*- 这里的代码需要下载讲义材料查看

use patent_ipc_pairs.dta, clear

这部分代码里面，“city_ipc.dta” 是要输入的数据，“city_ipc_pairs.dta”是保存的结果，city 是分组变量，ipc 是要组合的变量。使用这段代码就可以按照 city 变量进行分组，获取每个组内 ipc 变量的所有组合。

其实上面的 ipc 拆分和转置使用 R 语言效率也更高些：

rcall vanilla: install.packages("tidytext", repos = "https://mirrors.ustc.edu.cn/CRAN/", dependencies = T)

use df1, clear
rcall vanilla: ///
    library(tidyverse); ///
    haven::read_dta("df1.dta") %>% ///
    tidytext::unnest_tokens(input = "ipc", output = "ipc", ///
                            token = stringr::str_split, ///
                           pattern = "; ", to_lower = F) %>% ///
    haven::write_dta("df1a.dta")

use df1a, clear
replace ipc = substr(ipc, 1, 4)
save df2, replace

再继续计算 IPC 共现矩阵：

use patent_ipc_pairs.dta, clear
contract ipc1 ipc2, freq(C_ij)
save ipc_pairs, replace
*> file ipc_pairs.dta saved
*- 计算每个 IPC 的专利总数（用于归一化）
use df2, clear
*> (数据处理：微信公众号 RStata)
duplicates drop newipzlid ipc, force
*> Duplicates in terms of newipzlid ipc
*> (41,568 observations deleted)
keep newipzlid ipc
contract ipc, freq(total_patents)
drop if ipc == "  C1"
*> (1 observation deleted)
save ipc_total_patents, replace
*> file ipc_total_patents.dta saved

构建技术关联矩阵 Φ：

use ipc_pairs, clear
ren ipc1 ipc
merge m:1 ipc using ipc_total_patents, nogen
*> (variable ipc was str4, now str20 to accommodate using data's values)
*>     Result                      Number of obs
*>     -----------------------------------------
*>     Not matched                            16
*>         from master                         0
*>         from using                         16
*>     Matched                             8,890
*>     -----------------------------------------
rename total_patents total_i
ren ipc ipc1
ren ipc2 ipc
merge m:1 ipc using ipc_total_patents, nogen
*> (variable ipc was str4, now str20 to accommodate using data's values)
*>     Result                      Number of obs
*>     -----------------------------------------
*>     Not matched                            32
*>         from master                        16
*>         from using                         16
*>     Matched                             8,890
*>     -----------------------------------------
rename total_patents total_j
ren ipc ipc2
drop if mi(C_ij)
*> (32 observations deleted)
egen max_patent = rowmax(total_i total_j)
gen phi_ij = C_ij / max_patent
keep ipc1 ipc2 phi_ij
*- 标准化
egen min = min(phi_ij)
egen max = max(phi_ij)
replace phi_ij = (phi_ij - min) / (max - min)
*> (8,890 real changes made)
drop min max
save phi_mat, replace
*> file phi_mat.dta saved
*- 转换成矩阵
spread ipc2 phi_ij

foreach i of varlist _all {
    cap replace `i' = 0 if mi(`i')
}
save phi_mat_wide, replace

步骤6：构建每个城市每年的技术优势向量 M_mat

use phi_mat, clear
keep ipc1
ren ipc1 ipc
duplicates drop _all, force
*> Duplicates in terms of ipc
*> (8,303 observations deleted)
save ipclist, replace
*> file ipclist.dta saved
use rpc_data, clear
*> (数据处理：微信公众号 RStata)
keep year city ipc x_id
joinby ipc using ipclist
spread ipc x_id

foreach i of varlist _all {
    cap replace `i' = 0 if mi(`i')
}
save M_mat, replace

步骤7：计算两个城市 c 和 d 之间的技术互补指数 comp_cd

*- 假设我们选择年份 2001，北京市和上海市
use phi_mat_wide, clear
drop ipc1
mkmat _all, mat(phi_mat)

*- 使用 Mata
use M_mat, clear
mata:
    citylist = st_sdata(., "city")
    M_mat = st_data(., 3..st_nvar())
    phi_mat = st_matrix("phi_mat")
    city_c = "北京市"
    city_d = "上海市"

    ....
    *- 这里的代码需要下载讲义材料查看

    // 计算分子：M_c_tilde * Phi * M_d_tilde
    numerator = M_c_tilde * phi_mat * M_d_tilde'

    // 计算分母：sqrt(|M_c_tilde| * |M_d_tilde|)
    denominator = sqrt(sum(M_c_tilde) * sum(M_d_tilde))

    comp_cd = numerator / denominator
    comp_cd

    // .0537670076
end

由此就可以编写一个计算两个城市间技术互补指数的函数了：

use M_mat, clear
mata:
    citylist = st_sdata(., "city")
    M_mat = st_data(., 3..st_nvar())
    phi_mat = st_matrix("phi_mat")

    comp_cd_mat = J(0, 3, "")
    for (c = 1; c <= rows(citylist); c++){
        printf("%f\n", c)
        for (d = 1; d <= rows(citylist); d++){
            ....
            *- 这里的代码需要下载讲义材料查看

            comp_cd = numerator / denominator
            comp_cd_mat = comp_cd_mat \ (city_c, city_d, strofreal(comp_cd))
        }
    }
    st_matrix("obs", rows(comp_cd_mat))
    stata("clear")
    st_addvar("str100", "city1")
    st_addvar("str100", "city2")
    st_addvar("str100", "comp_cd")
    stata("set obs `=obs[1,1]'")
    st_sstore(., ("city1", "city2", "comp_cd"), comp_cd_mat)
end

compress
destring, replace
drop if city1 == city2
sum comp_cd
*>     Variable |        Obs        Mean    Std. dev.       Min        Max
*> -------------+---------------------------------------------------------
*>      comp_cd |    106,712    .0257762      .02747          0   .5128884

save 2001年各城市城市间技术互补指数, replace

如果想要获取所有年份的数据，循环年份即可，不过需要注意技术关联矩阵是需要基于所有年份的专利数据重新计算。

点击这里跳转到 RStata 短书平台获取附件：名师讲堂｜使用 Stata 计算地区间技术互补指数