为了让大家更好的理解本文内容，欢迎各位名师讲堂会员参加明晚 8 点的直播课：「使用 Stata 计算 Shapley-Shubik 权力指数：以上市公司前 10 大股东份额数据为例」。
Shapley-Shubik 指数是一种衡量投票权力的指标，它基于博弈论中的 Shapley 值的概念。Shapley-Shubik 指数可以用来分析在一个决策过程中，每个参与者对最终决策的贡献程度。
具体来说，Shapley-Shubik 指数反映了一个参与者在一个决策过程中成为“关键”参与者的概率。所谓“关键”参与者，是指如果该参与者不参与，原本可以通过的决议就无法通过。
Shapley-Shubik 指数的计算步骤如下:
列出所有可能的参与者排列组合。 对于每一种排列组合，计算该参与者是否为“关键”参与者。 统计每个参与者成为“关键”参与者的次数。 将每个参与者成为“关键”参与者的次数除以所有排列组合的总数，得到该参与者的 Shapley-Shubik 指数。
Shapley-Shubik 指数的取值范围是 [0, 1]。指数越大，表示该参与者在决策过程中的影响力越大。
在上市公司相关的研究中，Shapley-Shubik 指数可以用来衡量各股东的权力大小。本次课程中我们将以上市公司前 10 大股东份额数据为例，讲解如何在 Stata 中计算 Shapley-Shubik 指数。
一个简单的例子
考虑有四个团体，分别占据 40，39，11 和 11 个席位。总共是 101 个席位，决议通过的规则是 >= 51 票。计算每个团体的 Shapley-Shubik 权力指数。
一个简单的方法是借助这个在线应用：https://homepages.warwick.ac.uk/~ecaae/ssdirect.html
计算结果：
从 Shapley-Shubik 指数的计算过程中可以看到，如果增加一个团体，计算过程的耗时会增加 1 倍，所以我们还是使用 Mata 来提升计算速度。
mata:
y = (40 \ 39 \ 11 \ 11)
quota = 51
levels = colshape((1..rows(y)), 1)
n = factorial(rows(levels))
info = cvpermutesetup(levels, unique = 1)
i = 0
P = J(n, rows(levels), .)
while ((levels = cvpermute(info)) != J(0, 1, .)) {
++i
P[i,] = levels'
}
P
y[P[10,]]
// 循环所有的组合
R = J(n, 1, .)
for (i = 1; i <= n; i ++) {
R[i, 1] = P[i, colsum(runningsum(y[P[i,]]) :< quota) + 1]
}
R
Power = J(rows(y), 1, .)
for (i = 1; i <= rows(y); i ++) {
Power[i, 1] = colsum(R :== i)
}
Power :/ factorial(4)
end
cvpermute() 是用来获取所有的排列的，使用方法是这样的：
mata:
levels = (1 \ 2 \ 3 \ 4)
info = cvpermutesetup(levels, unique = 1)
cvpermute(info)
cvpermute(info)
cvpermute(info)
end
可以看到每次 cvpermute(info) 都会得到不同的排列，所以上面的代码就是把不同的排列组合成了一个 24 行的矩阵。
runningsum() 函数是用来计算矩阵的累加的：
mata: runningsum(1 \ 2 \ 3)
*> 1
*> +-----+
*> 1 | 1 |
*> 2 | 3 |
*> 3 | 6 |
*> +-----+
然后我们就可以编写一个函数封装这个过程了：
mata:
real colvector ShapleyShubik(real colvector y, real scalar quota) {
levels = colshape((1..rows(y)), 1)
n = factorial(rows(levels))
info = cvpermutesetup(levels, unique = 1)
i = 0
P = J(n, rows(levels), .)
while ((levels = cvpermute(info)) != J(0, 1, .)) {
++i
P[i,] = levels'
}
R = J(n, 1, .)
for (i = 1; i <= n; i ++) {
R[i, 1] = P[i, colsum(runningsum(y[P[i,]]) :< quota) + 1]
}
Power = J(rows(y), 1, .)
for (i = 1; i <= rows(y); i ++) {
Power[i, 1] = colsum(R :== i)
}
return(Power :/ n)
}
y = (40 \ 39 \ 11 \ 11)
ShapleyShubik(y, 51)
end
*> 1
*> +---------------+
*> 1 | .5 |
*> 2 | .1666666667 |
*> 3 | .1666666667 |
*> 4 | .1666666667 |
*> +---------------+
计算上市公司前 5 大股东的权力指数
为了节约时间，这里我只计算了前 5 大股东的。data.dta 中存储了一些上市公司的前 10 大股东份额数据：
use data.dta, clear
*- 计算前 5 大股东的 Shapley-Shubik 指数
drop top6-top10
*- 删除前 5 大股东总份额不足 50% 的
egen total = rowtotal(top*)
drop if total < 0.5
drop total
*- 缺失的都替换成 0
forval i = 1/5 {
replace top`i' = 0 if mi(top`i')
}
*- 创建一些空变量
forval i = 1/5 {
gen Power`i' = .
}
*- quota 都是 0.5
mata:
stock_mat = st_data(., 3..7)
Power_mat = J(rows(stock_mat), cols(stock_mat), .)
// 循环所有行
for (j = 1; j <= rows(stock_mat); j++) {
j
tempmat = stock_mat[j,]
res_mat = ShapleyShubik(colshape(tempmat, 1), 0.5)
Power_mat[j,] = rowshape(res_mat, 1)
}
Power_mat
// 存储到 dta 中
st_store(., 8..12, Power_mat)
end
save 计算结果, replace
编译 Mata 程序
由于 Mata 是编译型语言，所以我们也可以把 Mata 程序编译好再使用：
首先把下面的代码存储到 ShapleyShubik.mata 文件里面：
mata:
real colvector ShapleyShubik(real colvector y, real scalar quota) {
levels = colshape((1..rows(y)), 1)
n = factorial(rows(levels))
info = cvpermutesetup(levels, unique = 1)
i = 0
P = J(n, rows(levels), .)
while ((levels = cvpermute(info)) != J(0, 1, .)) {
++i
P[i,] = levels'
}
R = J(n, 1, .)
for (i = 1; i <= n; i ++) {
R[i, 1] = P[i, colsum(runningsum(y[P[i,]]) :< quota) + 1]
}
Power = J(rows(y), 1, .)
for (i = 1; i <= rows(y); i ++) {
Power[i, 1] = colsum(R :== i)
}
return(Power :/ n)
}
end
然后在 Stata 中运行下面的代码：
clear all
do ShapleyShubik.mata
lmbuild lShapleyShubik.mlib, replace
这样就会编译得到一个 lShapleyShubik.mlib 文件，默认会把这个文件存储到 Personal 文件夹里面。这样就不需要每次都运行下 ShapleyShubik() 函数的源代码了，可以直接使用了：
clear all
mata: ShapleyShubik((40 \ 39 \ 11 \ 11), 51)
*> 1
*> +---------------+
*> 1 | .5 |
*> 2 | .1666666667 |
*> 3 | .1666666667 |
*> 4 | .1666666667 |
*> +---------------+
也可以把这个 mlib 文件发送给别人使用，而且不会暴露源代码，这也是代码加密的一种方法。
*- 拷贝到当前文件夹里面
copy "`c(sysdir_personal)'/lShapleyShubik.mlib" lShapleyShubik.mlib, replace
直播信息
为了让大家更好的理解本文内容，欢迎各位名师讲堂会员参加明晚 8 点的直播课：「使用 Stata 计算 Shapley-Shubik 权力指数：以上市公司前 10 大股东份额数据为例」。
直播地址：腾讯会议(需要报名 RStata 培训班参加) 讲义材料：需要购买 RStata 名师讲堂会员，详情可阅读：一起来学习 R 语言和 Stata 啦！学习过程中遇到的问题也可以随时提问！
更多关于 RStata 会员的更多信息可添加微信号 r_stata 咨询：
附件下载（点击文末的阅读原文即可跳转）：
https://rstata.duanshu.com/#/brief/course/8165b097bc3d4e29975966cade18f830
