该课程已经讲解完，感兴趣的小伙伴可以文末的阅读原文跳转到 RStata 平台观看录播～
「经济研究」2023 年第 6 期中有篇论文<共同股东网络与国有企业的创新知识溢出>中使用了企业技术相似度的指标。具体来说是使用如下公式计算的：
该指标的取值范围是 [0, 1]，越接近 1 表示两个企业的技术空间越接近。其中 表示企业 j 在每个专利类别中的专利比例构成向量。
今天的课程中我们将以工企专利数据为例来演示如何在 Stata 中编程计算该指标，以及如何使用 Mata 程序来大大提升计算效率。
根据工企专利匹配结果计算每个企业的专利申请类别比例
附件中提供了一份 工企专利匹配2013(完整版).dta 文件，是 2013 年工企专利匹配结果，首先我们从这个数据出发计算每个企业的专利申请类别比例。这里我们使用专利分类号中的“部”来表示专利的申请类别。也就是主分类号的首字母：
A 部－人类生活必需（农、轻、医） B 部－作业、运输 C 部－化学、冶金 D 部－纺织、造纸 E 部－固定建筑物（建筑、采矿） F 部－机械工程 G 部－物理 H 部－电学
use "工企专利匹配2013(完整版).dta", clear
*- 保留需要的一些变量
keep 企业匹配唯一标识码 企业名称 行业大类代码 zlid 主分类号 专利类型
*- 去除主分类号后面的年份标志
replace 主分类号 = ustrregexs(1) if ustrregexm(主分类号, "(.*)\(")
*- 删除设计专利
drop if 专利类型 == "设计型专利"
*- 生成部变量
gen 部 = substr(主分类号, 1, 1)
codebook 部
*- 统计每部的专利占比
drop 主分类号 专利类型
collapse (count) n = zlid, by(企业匹配唯一标识码 企业名称 部 行业大类代码)
bysort 企业匹配唯一标识码 企业名称: egen sum = sum(n)
gen ratio = n / sum
drop n sum
*- 安装 tidy: ssc install tidy
spread 部 ratio
foreach i of varlist A-H {
replace `i' = 0 if missing(`i')
}
save "专利比例构成", replace
要计算的公司数量不多的情况
如果要计算的公司数量不多，我们可以直接在 dta 数据里面进行操作，这里我们选择前 100 家公司为例。如果你对矩阵运算比较难以理解可以采用这种方法。
为了构建“公司对”，我们需要准备两个内容相同但是变量名不相同的数据，data1 和 data2：
use 专利比例构成, clear
drop 行业大类代码
*- 使用前 100 家公司演示
keep in 1/100
gen 公司序号 = _n
order 公司序号
gather A-H
ren var 部
ren value ratio
drop 企业匹配唯一标识码 企业名称
save data1, replace
list in 1/10
*> +-----------------------+
*> | 公司序号 部 ratio |
*> |-----------------------|
*> 1. | 1 A 0 |
*> 2. | 1 B 0 |
*> 3. | 1 C 0 |
*> 4. | 1 D 0 |
*> 5. | 1 E 0 |
*> |-----------------------|
*> 6. | 1 F 0 |
*> 7. | 1 G 0 |
*> 8. | 1 H 1 |
*> 9. | 2 A 0 |
*> 10. | 2 B 1 |
*> +-----------------------+
data2 只是把 data1 的变量名改一下：
foreach i of varlist _all {
ren `i' `i'2
}
save data2, replace
然后使用 cross using 把两个数据交叉起来，这样就可以生成一一对应的“公司对”数据了：
use data1, clear
cross using data2
里面是有重复的，例如 公司1-公司2 和 公司2-公司1 是重复的，公司1-公司1 这样的我们也是不需要的：
*- 删除掉重复的企业对
keep if 公司序号 <= 公司序号2 & 部 == 部2
分别计算 Tii、Tjj、Tij：
gen tii = ratio * ratio
gen tjj = ratio2 * ratio2
gen tij = ratio * ratio2
然后再分组汇总即可：
collapse (sum) tii (sum) tjj (sum) tij, by(公司序号 公司序号2)
按照最开头的公式：
gen tech = tij / (sqrt(tii) * sqrt(tjj))
drop tii tjj tij
save "技术相似度", replace
这里我们使用的是 100 家公司，最终创建的交叉数据集里面有 64w 个观测值，可以想象，如果我们使用全部的 4 万个企业，最终创建的交叉数据集里面会有超过 102400000000，也就是 1024 亿条观测值。因此这种方法只能用于企业数量很少的时候。
要计算的公司数量较多的情况
下面我们再看如何使用矩阵运算的方法来计算这一指标。
保留 1000 家公司：
use 专利比例构成, clear
drop 行业大类代码
keep in 1/1000
如果要计算的企业数量很多，需要设置更大的变量数上限，例如计算 4 万家公司的：
set maxvar 42000
*- 通常 maxvar 的默认值是 5000，所以这里并不需要设置
在写代码之前我们先研究下思路，首先我们可以把数据中的 A-H 转换成矩阵，这个矩阵的含义会是下面这样的：
如果我们把这个矩阵和其转置矩阵相乘，可以得到：
这样的话，企业 i 和 企业 j 的技术相似度就是 了。下面我们开始写程序：
首先使用 mkmat 可以把 dta 里面的变量转换成矩阵：
mkmat A - H, mat(totalmat)
然后计算乘积：
mat z = totalmat * totalmat'
创建一个空矩阵用于存储技术相似度：
*- 创建一个 1000x1000 的空矩阵
mat tech = J(1000, 1000, .)
使用一个双层的循环来填充里面的结果：
forval i = 1/1000 {
di "`i'"
forval j = 1/1000 {
mat tech[`i', `j'] = z[`i', `j'] / (sqrt(z[`i', `i']) * sqrt(z[`j', `j']))
}
}
svmat 命令可以把矩阵再转换成 dta 数据：
svmat tech
drop A - H
gen 公司序号 = _n
order 企业匹配唯一标识码 企业名称 公司序号
save 技术相似度2, replace
如果想计算每个公司与其他公司的平均技术相似度可以使用行运算：
use 技术相似度2, clear
*- 计算 行均值
egen techmean = rowtotal(tech*)
*- 减去本身的
replace techmean = techmean - 1
*- 除以（公司数量 - 1）
replace techmean = techmean / (`=_N' - 1)
drop tech1 - tech`=_N'
save 平均技术相似度, replace
上面在计算 tech 的时候也可以更快点，只填补需要的一半，不过这种就不便于求平均技术相似度了：
mat tech = J(1000, 1000, .)
forval i = 1/1000 {
di "`i'"
forval j = `=`i' + 1'/1000 {
mat tech[`i', `j'] = z[`i', `j'] / (sqrt(z[`i', `i']) * sqrt(z[`j', `j']))
}
}
尽管这种方法提速很多，但是计算全部 4 万家还是非常困难的，幸运的是，更多的时候我们并不会直接计算所有公司的，而是计算每个行业内的。
对于这种情况，如果企业较多可以把数据拆分成多个文件，然后再逐个文件的计算：
use 专利比例构成, clear
cap mkdir "分行业大类"
levelsof 行业大类代码, local(codelist)
use 专利比例构成, clear
foreach i in `codelist' {
preserve
keep if 行业大类代码 == `i'
save "分行业大类/`i'", replace
restore
}
然后逐个文件计算即可。首先查看企业最多的行业：
use 专利比例构成, clear
gen id = _n
collapse (count) id, by(行业大类代码)
gsort -id
最多的行业有 5178 个公司，所以这里：
clear all
set maxvar 5200
索引所有的文件可以使用上面的循环，也可以使用下面的方式：
这个还是先别运行，会挺耗费时间
local files: dir "分行业大类" files "*.dta"
di `"`files'"'
cap mkdir "分行业技术相似度技术指标"
foreach f in `files' {
use "分行业大类/`f'", clear
di "`f' 文件公司数量：`=_N'"
mkmat A - H, mat(totalmat)
mat z = totalmat * totalmat'
mat tech = J(`=_N', `=_N', .)
forval i = 1/`=_N' {
forval j = 1/`=_N' {
mat tech[`i', `j'] = z[`i', `j'] / (sqrt(z[`i', `i']) * sqrt(z[`j', `j']))
}
}
svmat tech
drop A - H
gen 公司序号 = _n
order 企业匹配唯一标识码 企业名称 行业大类代码 公司序号
egen techmean = rowtotal(tech*)
replace techmean = techmean - 1
replace techmean = techmean / (`=_N' - 1)
drop tech1 - tech`=_N' 公司序号
save "分行业技术相似度技术指标/`f'", replace
}
合并计算结果：
local files: dir "分行业技术相似度技术指标" files "*.dta"
di `"`files'"'
use 分行业技术相似度技术指标/38.dta, clear
drop in 1/`=_N'
foreach f in `files'{
append using "分行业技术相似度技术指标/`f'"
}
save "计算结果", replace
使用该方法尽管高效了很多，但是还是挺慢的。下面我们再看如何使用 Mata 提速。
使用 Mata 程序计算企业的技术相似度
大家通常认为 Mata 是一种矩阵语言，Stata 公司也是这么定位的。更准确地说，Mata 是一中跨平台的便携式的、编译的编程语言。值得注意的是，Mata 是编译型语言，这意味着程序可以运行的很快。而前面用的 ado 则是一种解释型语言，也就是说 ado 很慢。总体上来说，Mata 比 ado 快 10~40 倍，下面的对比结果也印证了这一结论。
这里我们同样也适用 1000 家企业为例：
clear all
use 专利比例构成, clear
keep in 1/1000
mata:
totalmat = st_data(., ("A", "B", "C", "D", "E", "F", "G", "H"))
z = totalmat * totalmat'
n = st_nobs()
tech = J(n, n, .)
for (i = 1; i <= n; i++) {
for (j = 1; j <= n; j++) {
tech[i, j] = z[i, j] / (sqrt(z[i, i]) * sqrt(z[j, j]))
}
}
st_matrix("tech", tech)
end
mat list tech
drop A-H
svmat tech
可以对比下 Stata 中的矩阵运算和 Mata 中的矩阵运算耗时：
clear all
use 专利比例构成, clear
keep in 1/2000
*- Mata 算法
timer on 1
qui {
mata:
totalmat = st_data(., ("A", "B", "C", "D", "E", "F", "G", "H"))
z = totalmat * totalmat'
n = st_nobs()
tech = J(n, n, .)
for (i = 1; i <= n; i++) {
for (j = 1; j <= n; j++) {
tech[i, j] = z[i, j] / (sqrt(z[i, i]) * sqrt(z[j, j]))
}
}
st_matrix("tech", tech)
end
}
timer off 1
*- ado 算法
timer on 2
qui {
mkmat A - H, mat(totalmat)
mat z = totalmat * totalmat'
local n = `=_N'
mat tech = J(`n', `n', .)
forval i = 1/`n' {
forval j = 1/`n' {
mat tech[`i', `j'] = z[`i', `j'] / (sqrt(z[`i', `i']) * sqrt(z[`j', `j']))
}
}
}
timer off 2
timer list
ret list
*> scalars:
*> r(N) = 1000
*> r(t1) = .324
*> r(nt1) = 1
*> r(t2) = 15.614
*> r(nt2) = 1
di r(t2) / r(t1)
*> 48.191358
在我的电脑上 Mata 程序仅用时 0.32s，ado 程序却用了超过 15.61s，相差大概 50 倍，如果把企业数设置成 2000 这一比例会变成 58 倍。因此如果企业数量较多还是使用 Mata 程序更恰当。
最后我们再循环处理分行业的结果：
clear all
local files: dir "分行业大类" files "*.dta"
di `"`files'"'
cap mkdir "分行业技术相似度技术指标"
*- 由于 mata 程序里面有 end，所以不能直接用在循环里面，首先定义一个 mata 函数，void 表示这是个无返回值的函数。
mata:
void matfun() {
totalmat = st_data(., ("A", "B", "C", "D", "E", "F", "G", "H"))
z = totalmat * totalmat'
n = st_nobs()
tech = J(n, n, .)
for (i = 1; i <= n; i++) {
for (j = 1; j <= n; j++) {
tech[i, j] = z[i, j] / (sqrt(z[i, i]) * sqrt(z[j, j]))
}
}
st_matrix("tech", tech)
}
end
foreach f in `files' {
qui use "分行业大类/`f'", clear
di "`f' 文件公司数量：`=_N'"
mata: matfun()
qui {
svmat tech
drop A - H
gen 公司序号 = _n
order 企业匹配唯一标识码 企业名称 行业大类代码 公司序号
egen techmean = rowtotal(tech*)
replace techmean = techmean - 1
replace techmean = techmean / (`=_N' - 1)
drop tech1 - tech`=_N' 公司序号
save "分行业技术相似度技术指标/`f'", replace
}
}
*- 合并结果
local files: dir "分行业技术相似度技术指标" files "*.dta"
di `"`files'"'
use 分行业技术相似度技术指标/38.dta, clear
drop in 1/`=_N'
foreach f in `files'{
append using "分行业技术相似度技术指标/`f'"
}
save "分行业技术相似度技术指标计算结果", replace
gsort -techmean
在掌握 mata 方法的基础上还是建议大家用 mata 方法，确实快了很多。
直播信息
该课程已经讲解完，感兴趣的小伙伴可以文末的阅读原文跳转到 RStata 平台观看录播～
直播地址：腾讯会议(需要报名 RStata 培训班参加) 讲义材料：需要购买 RStata 名师讲堂会员，详情可阅读：一起来学习 R 语言和 Stata 啦！学习过程中遇到的问题也可以随时提问！
更多关于 RStata 会员的更多信息可添加微信号 r_stata 咨询：
附件下载（点击文末的阅读原文即可跳转）：
https://rstata.duanshu.com/#/brief/course/11ec07a861384c218f5f77582875b0a2
