名师讲堂｜使用 Stata 处理上市公司专利互引矩阵及测算上市公司供应链协同创新程度（二）

上次课程主页： https://rstata.duanshu.com/#/brief/course/e58dd2c837c443469d65cf37c26d8b39

继续上次课程的内容，今天我们继续讲解根据上市公司专利互引信息测算供应链协同创新程度指标的方法。

上市公司供应链协同创新程度指标的构造方法是参考附件中的文献「共同股东与供应链协同创新——基于专利互引的视角_汤旭东.pdf」

上次的课程中我们已经构造好了上市公司专利互引的数据，接下来把供应商和客户数据分别和该数据匹配，得到了供应商和客户互引的数据，再统计数量就可以得到上市公司供应链协同创新程度指标了。

首先回顾上次的上市公司专利互引数据：

use 上市公司之间专利互引信息.dta, clear
*> (数据处理：微信公众号 RStata)

供应商和客户数据使用的是之前分享的这个：

2001~2023 年上市公司前5大供应商、客户与工商注册信息匹配结果（含经纬度及所处的省市区县）: https://rstata.duanshu.com/#/brief/course/c0b17ab7ab2945ffaea746dec9380a8e

为了数据不至于太大，我保留了一些关键变量和 2010～2023 年的数据：

• 2010～2023年上市公司供应商数据.dta
• 2010～2023年上市公司客户数据.dta

仅仅保留是上市公司的供应商和客户：

use "2010～2023年上市公司供应商数据", clear
*> (数据处理：微信公众号 RStata。注意经纬度解析是没法保证百分百准确的，所以结果中会存在一些错误的解析结果，介意勿用。)
keep if !mi(公司股票代码)
*> (147,221 observations deleted)
keep gysid 供应商公司ID 公司股票代码 供应商名称
gen code1 = ustrregexs(0) if ustrregexm(公司股票代码, "d{6}")
*> (72 missing values generated)
replace 公司股票代码 = subinstr(公司股票代码, code1, "", .)
*> (2,172 real changes made)
gen code2 = ustrregexs(0) if ustrregexm(公司股票代码, "d{6}")
*> (2,204 missing values generated)
replace 公司股票代码 = subinstr(公司股票代码, code2, "", .)
*> (40 real changes made)
drop 公司股票代码
gather code*
*> gysid 供应商公司ID 供应商名称
drop if mi(value)
*> (2,276 observations deleted)
drop var
ren value 股票代码
duplicates drop _all, force
*> Duplicates in terms of gysid 供应商公司ID 供应商名称 股票代码
*> (0 observations are duplicates)
save 供应商数据, replace
*> file 供应商数据.dta saved
list in 1/10
*>      +---------------------------------------------------------------------+
*>      |      gysid   供应商~D                         供应商名称   股票代码 |
*>      |---------------------------------------------------------------------|
*>   1. | 2010000001     103612           中国东方航空股份有限公司     600115 |
*>   2. | 2010000045   10515296           浙江物产燃料集团有限公司     603071 |
*>   3. | 2010000078   10107077           苏州纽威阀门股份有限公司     603699 |
*>   4. | 2010000090     105971                   同方股份有限公司     600100 |
*>   5. | 2010000095     102735   内蒙古霍林河露天煤业股份有限公司     002128 |
*>      |---------------------------------------------------------------------|
*>   6. | 2010000119   10145920       广州广电计量测试技术有限公司     002967 |
*>   7. | 2010000173     103149         东方日升新能源股份有限公司     300118 |
*>   8. | 2010000183     101498               沈阳化工股份有限公司     000698 |
*>   9. | 2010000190     105897               海南航空股份有限公司     600221 |
*>  10. | 2010000190     105897               海南航空股份有限公司     900945 |
*>      +---------------------------------------------------------------------+
use "2010～2023年上市公司客户数据", clear
*> (数据处理：微信公众号 RStata。注意经纬度解析是没法保证百分百准确的，所以结果中会存在一些错误的解析结果，介意勿用。)
keep if !mi(公司股票代码)
*> (203,437 observations deleted)
keep khid 客户公司ID 公司股票代码 客户名称
gen code1 = ustrregexs(0) if ustrregexm(公司股票代码, "d{6}")
*> (354 missing values generated)
replace 公司股票代码 = subinstr(公司股票代码, code1, "", .)
*> (4,109 real changes made)
gen code2 = ustrregexs(0) if ustrregexm(公司股票代码, "d{6}")
*> (4,287 missing values generated)
replace 公司股票代码 = subinstr(公司股票代码, code2, "", .)
*> (176 real changes made)
drop 公司股票代码
gather code*
*> khid 客户公司ID 客户名称
drop if mi(value)
*> (4,641 observations deleted)
drop var
ren value 股票代码
duplicates drop _all, force
*> Duplicates in terms of khid 客户公司ID 客户名称 股票代码
*> (0 observations are duplicates)
save 客户数据, replace
*> file 客户数据.dta saved
list in 1/10
*>      +---------------------------------------------------------------+
*>      |       khid   客户公~D                     客户名称   股票代码 |
*>      |---------------------------------------------------------------|
*>   1. | 2010000012   10134538         海南灵康制药有限公司     603669 |
*>   2. | 2010000024     105019         东风汽车股份有限公司     600006 |
*>   3. | 2010000128     101326         本钢板材股份有限公司     000761 |
*>   4. | 2010000128     101326         本钢板材股份有限公司     200761 |
*>   5. | 2010000154     101267         新兴铸管股份有限公司     000778 |
*>      |---------------------------------------------------------------|
*>   6. | 2010000169     101752     重庆长安汽车股份有限公司     000625 |
*>   7. | 2010000169     101752     重庆长安汽车股份有限公司     200625 |
*>   8. | 2010000208     102498     浙江双箭橡胶股份有限公司     002381 |
*>   9. | 2010000227      10973   九州通医药集团股份有限公司     600998 |
*>  10. | 2010000252     103367     中国神华能源股份有限公司     601088 |
*>      +---------------------------------------------------------------+

为了计算供应商和客户的互引数量，我们需要分别计算供应商引用客户的数量以及客户引用供应商的数量。

首先计算供应商引用客户的，也就是把供应商和股票代码_原专利匹配，客户与股票代码_被引专利匹配：

*- 供应商引用客户
use 供应商数据, clear
*> (数据处理：微信公众号 RStata。注意经纬度解析是没法保证百分百准确的，所以结果中会存在一些错误的解析结果，介意勿用。)
ren 股票代码 股票代码_原专利
gen year = substr(gysid, 1, 4)
keep year 供应商公司ID 股票代码_原专利
duplicates drop _all, force
*> Duplicates in terms of 供应商公司ID 股票代码_原专利 year
*> (265 observations deleted)
destring year, replace
*> year: all characters numeric; replaced as int
save 供应商数据1, replace
*> file 供应商数据1.dta saved
*- 客户专利的年份不重要，就删除了
use 客户数据, clear
*> (数据处理：微信公众号 RStata。注意经纬度解析是没法保证百分百准确的，所以结果中会存在一些错误的解析结果，介意勿用。)
drop khid
ren 股票代码 股票代码_被引专利
duplicates drop 股票代码_被引专利, force
*> Duplicates in terms of 股票代码_被引专利
*> (3,257 observations deleted)
save 客户数据1, replace
*> file 客户数据1.dta saved

关于年份的处理也要格外小心，前面也提到过，对于专利被引用量来说，随着年份的延长，肯定是增加的，所以通常不会直接统计专利的被引用量，而是统计当年或者三年内的，所以这里我们仅仅保留当年引用当年的：

use 上市公司之间专利互引信息.dta, clear
*> (数据处理：微信公众号 RStata)
gen year = substr(newipzlid_原专利, 1, 4)
gen year2 = substr(newipzlid_被引专利, 1, 4)
destring year, replace
*> year: all characters numeric; replaced as int
save 上市公司之间专利互引信息1.dta, replace
*> file 上市公司之间专利互引信息1.dta saved
use 供应商数据1, clear
*> (数据处理：微信公众号 RStata。注意经纬度解析是没法保证百分百准确的，所以结果中会存在一些错误的解析结果，介意勿用。)
joinby 股票代码_原专利 year using 上市公司之间专利互引信息1.dta
*- 匹配客户
joinby 股票代码_被引专利 using 客户数据1
duplicates drop 供应商公司ID year newipzlid_原专利 newipzlid_被引专利 股票代码_被引专利, force
*> Duplicates in terms of 供应商公司ID year newipzlid_原专利 newipzlid_被引专利 股票代码_被引专利
*> (1,376 observations deleted)
ren year 年份

由于专利数据中同时包含了专利的引用与授权公告，直接删除授权公告又担心有些授权专利的申请公告没有包含在数据里面：

统计的时候可以先去除公开公告号里面的 A、B、U、S。其中 A 代表发明专利的申请公开，B 代表发明专利的授权公告，U 代表实用新型专利的授权公告，S 代表外观设计专利的授权公告。

*- 去除重复专利
replace 公开公告号_原专利 = subinstr(公开公告号_原专利, "A", "", .)
*> (12,707 real changes made)
replace 公开公告号_原专利 = subinstr(公开公告号_原专利, "B", "", .)
*> (6,939 real changes made)
replace 公开公告号_原专利 = subinstr(公开公告号_原专利, "U", "", .)
*> (52 real changes made)
replace 公开公告号_原专利 = subinstr(公开公告号_原专利, "S", "", .)
*> (0 real changes made)
replace 公开公告号_被引专利 = subinstr(公开公告号_被引专利, "A", "", .)
*> (6,846 real changes made)
replace 公开公告号_被引专利 = subinstr(公开公告号_被引专利, "B", "", .)
*> (3,896 real changes made)
replace 公开公告号_被引专利 = subinstr(公开公告号_被引专利, "U", "", .)
*> (8,953 real changes made)
replace 公开公告号_被引专利 = subinstr(公开公告号_被引专利, "S", "", .)
*> (3 real changes made)
duplicates drop 供应商公司ID 公开公告号_原专利 公开公告号_被引专利 股票代码_被引专利, force
*> Duplicates in terms of 供应商公司ID 公开公告号_原专利 公开公告号_被引专利 股票代码_被引专利
*> (9,855 observations deleted)
*- 统计各年供应商引用客户专利的数量
contract 年份 供应商公司ID 客户公司ID
order 年份
ren _freq 供应商引用客户专利数量
save 供应商引用客户专利数量, replace
*> file 供应商引用客户专利数量.dta saved

类似的方法再统计客户引用供应商的：

*- 客户引用供应商
use 客户数据, clear
*> (数据处理：微信公众号 RStata。注意经纬度解析是没法保证百分百准确的，所以结果中会存在一些错误的解析结果，介意勿用。)
ren 股票代码 股票代码_原专利
gen year = substr(khid, 1, 4)
keep year 客户公司ID 股票代码_原专利
duplicates drop _all, force
*> Duplicates in terms of 客户公司ID 股票代码_原专利 year
*> (1,611 observations deleted)
destring year, replace
*> year: all characters numeric; replaced as int
save 客户数据2, replace
*> file 客户数据2.dta saved
*- 供应商专利的年份不重要，就删除了，否则后面的数据量会指数上涨
use 供应商数据, clear
*> (数据处理：微信公众号 RStata。注意经纬度解析是没法保证百分百准确的，所以结果中会存在一些错误的解析结果，介意勿用。)
drop gysid
ren 股票代码 股票代码_被引专利
duplicates drop 股票代码_被引专利, force
*> Duplicates in terms of 股票代码_被引专利
*> (1,391 observations deleted)
save 供应商数据2, replace
*> file 供应商数据2.dta saved
use 客户数据2, clear
*> (数据处理：微信公众号 RStata。注意经纬度解析是没法保证百分百准确的，所以结果中会存在一些错误的解析结果，介意勿用。)
joinby 股票代码_原专利 year using 上市公司之间专利互引信息1.dta
*- 匹配供应商
joinby 股票代码_被引专利 using 供应商数据2
duplicates drop 客户公司ID year newipzlid_原专利 newipzlid_被引专利 股票代码_被引专利, force
*> Duplicates in terms of 客户公司ID year newipzlid_原专利 newipzlid_被引专利 股票代码_被引专利
*> (1,455 observations deleted)
ren year 年份
*- 去除重复专利
replace 公开公告号_原专利 = subinstr(公开公告号_原专利, "A", "", .)
*> (14,937 real changes made)
replace 公开公告号_原专利 = subinstr(公开公告号_原专利, "B", "", .)
*> (8,567 real changes made)
replace 公开公告号_原专利 = subinstr(公开公告号_原专利, "U", "", .)
*> (64 real changes made)
replace 公开公告号_原专利 = subinstr(公开公告号_原专利, "S", "", .)
*> (0 real changes made)
replace 公开公告号_被引专利 = subinstr(公开公告号_被引专利, "A", "", .)
*> (8,197 real changes made)
replace 公开公告号_被引专利 = subinstr(公开公告号_被引专利, "B", "", .)
*> (4,650 real changes made)
replace 公开公告号_被引专利 = subinstr(公开公告号_被引专利, "U", "", .)
*> (10,715 real changes made)
replace 公开公告号_被引专利 = subinstr(公开公告号_被引专利, "S", "", .)
*> (6 real changes made)
duplicates drop 客户公司ID 公开公告号_原专利 公开公告号_被引专利 股票代码_被引专利, force
*> Duplicates in terms of 客户公司ID 公开公告号_原专利 公开公告号_被引专利 股票代码_被引专利
*> (11,983 observations deleted)
*- 统计各年客户引用供应商专利的数量
contract 年份 客户公司ID 供应商公司ID
order 年份
ren _freq 客户引用供应商专利数量
save 客户引用供应商专利数量, replace
*> file 客户引用供应商专利数量.dta saved

最后匹配两部分数据即可计算协同创新程度指标了：

*- 匹配两部分数据
use 供应商引用客户专利数量, clear
*> (数据处理：微信公众号 RStata。注意经纬度解析是没法保证百分百准确的，所以结果中会存在一些错误的解析结果，介意勿用。)
merge 1:1 年份 客户公司ID 供应商公司ID using 客户引用供应商专利数量
*>     Result                      Number of obs
*>     -----------------------------------------
*>     Not matched                         1,652
*>         from master                       718  (_merge==1)
*>         from using                        934  (_merge==2)
*>     Matched                               211  (_merge==3)
*>     -----------------------------------------
drop _m
replace 供应商引用客户专利数量 = 0 if mi(供应商引用客户专利数量)
*> (934 real changes made)
replace 客户引用供应商专利数量 = 0 if mi(客户引用供应商专利数量)
*> (718 real changes made)
*- 供应商公司ID 和 客户公司ID 是按照同一规则编码的，相同的时候表示是同一家公司，应该删除
drop if 供应商公司ID == 客户公司ID
*> (625 observations deleted)
gen CoCite = log((供应商引用客户专利数量 + 客户引用供应商专利数量)/2 + 1)
label data "数据计算：微信公众号 RStata"
sum CoCite
*>     Variable |        Obs        Mean    Std. dev.       Min        Max
*> -------------+---------------------------------------------------------
*>       CoCite |      1,238    .6137653    .4277335   .4054651   3.314186
save 2010～2023年供应链协同创新程度面板数据, replace
*> file 2010～2023年供应链协同创新程度面板数据.dta saved
list in 1/10
*>      +-----------------------------------------------------------+
*>      | 年份   供应商~D   客户公~D   供应~量   客户~量     CoCite |
*>      |-----------------------------------------------------------|
*>   1. | 2010     105971     101553         1         0   .4054651 |
*>   2. | 2010     104710     102088         3         2   1.252763 |
*>   3. | 2010     105971     102144         1         0   .4054651 |
*>   4. | 2010   10131196     102620         1         0   .4054651 |
*>   5. | 2010     101969     103276         1         0   .4054651 |
*>      |-----------------------------------------------------------|
*>   6. | 2011     106457     101553         3         4   1.504077 |
*>   7. | 2011   10345164     101553         3         0   .9162908 |
*>   8. | 2011     101550     101969         1         0   .4054651 |
*>   9. | 2011     103002     103367         1         0   .4054651 |
*>  10. | 2011     102088     103740         4         0   1.098612 |
*>      +-----------------------------------------------------------+

这样我们就计算得到了这一指标。

点击这里跳转到 RStata 短书平台获取附件：名师讲堂｜使用 Stata 处理上市公司专利互引矩阵及测算上市公司供应链协同创新程度（二）