名师讲堂｜使用 Stata 计算企业技术相似度：交叉数据集、矩阵运算和 Mata 方法

「经济研究」2023 年第 6 期中有篇论文<共同股东网络与国有企业的创新知识溢出>中使用了企业技术相似度的指标。具体来说是使用如下公式计算的：

该指标的取值范围是 [0, 1]，越接近 1 表示两个企业的技术空间越接近。其中 Tj 表示企业 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 转换成矩阵，这个矩阵的含义会是下面这样的：

如果我们把这个矩阵和其转置矩阵相乘，可以得到：

下面我们开始写程序：

首先使用 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 短书平台获取附件：名师讲堂｜使用 Stata 计算企业技术相似度：交叉数据集、矩阵运算和 Mata 方法