附件中的文献《生产性服务业集聚何以赋能技术扩散》中提到了城市的生产性服务业专业化集聚指标和多样化集聚指标，今天的课程中我们将会讲解如何使用 Stata 语言根据城市统计年鉴数据来测算这两个指标。
城市统计年鉴数据在这里：
1999～2025 年中国城市统计年鉴面板数据整理结果: https://rstata.duanshu.com/#/brief/course/c6b0aaf3bcba494faadbacacb012dcef
这里选择的是 2017～2019 年的样本作为演示使用。也就是附件中的 城市统计年鉴样本数据.dta。
一、指标来源与计算公式
1.1 专业化集聚指数 SP（公式4）
1.2 多样化集聚指数 DV（公式5）
二、数据与变量说明
2.1 数据来源
城市统计年鉴数据在这里：
1999～2025 年中国城市统计年鉴面板数据整理结果: https://rstata.duanshu.com/#/brief/course/c6b0aaf3bcba494faadbacacb012dcef
这里选择的是 2017～2019 年的样本作为演示使用。也就是附件中的 城市统计年鉴样本数据.dta。
由于计算所需的变量主要集中在 2003～2019 年，所以完整数据仅可计算得到 2003～2019 年的。
2.2 生产性服务业行业界定
参考 Ke et al.（2014）和韩峰和阳立高（2020），生产性服务业包括 7个细分行业：
| 代码 |
行业名称 |
说明 |
| s1 |
交通运输、仓储和邮政业 |
全市口径 |
| s2 |
金融业 |
全市口径 |
| s3 |
科学研究和技术服务业 |
2003-2016含地质勘探，2017-2019不含 |
| s4 |
租赁和商务服务业 |
全市口径 |
| s5 |
信息传输、计算机服务和软件业 |
全市口径 |
| s6 |
批发和零售业 |
全市口径 |
| s7 |
水利、环境和公共设施管理业 |
全市口径 |
注意：科学研究和技术服务业在2003-2016年统计口径含”地质勘探业”，2017-2019年不含。处理时用 coalesce() 合并两套口径。
2.3 就业数据说明
- 数据来自《中国城市统计年鉴》城镇单位从业人员（全市口径）
- 单位：人
- 总就业 = 城镇单位从业人员期末人数（全市）
三、Stata 代码实现
3.1 读取数据与变量重命名
cd "~/Desktop/使用 Stata 测算各城市生产性服务业集聚水平（Mata 版本）"
use "城市统计年鉴样本数据.dta", clear
rename 年份 year
rename 省 prov
rename 省代码 prov_code
rename 市 city
rename 市代码 cityid
rename 城镇单位从业人员期末人数_人_全市 emp_total
rename 第三产业_交通运输仓储和邮政业_人_全市 emp_s1
rename 第三产业_金融业_人_全市 emp_s2
rename 第三产业_科学研究和技术服务和地质勘探业_人_全市 emp_s3a
rename 第三产业_科学研究和技术服务业_人_全市 emp_s3b
rename 第三产业_租赁和商务服务业_人_全市 emp_s4
rename 第三产业_信息传输计算机服务和软件业_人_全市 emp_s5
rename 第三产业_批发和零售业_人_全市 emp_s6
rename 第三产业_水利环境和公共设施管理业_人_全市 emp_s7
gen emp_s3 = emp_s3a
replace emp_s3 = emp_s3b if missing(emp_s3) & !missing(emp_s3b)
drop emp_s3a emp_s3b
|
3.2 数据清洗与筛选
*- 此处代码需下载讲义材料查看～
disp "样本量（城市×年份）: " _N
summ year, meanonly
disp "年份范围: " r(min) " ~ " r(max)
|
3.3 异常值处理
summ emp_total if city == "嘉峪关市" & year == 2013, meanonly
local jyg_raw = r(mean)
summ emp_total if city == "嘉峪关市" & inlist(year, 2012, 2014), meanonly
local jyg_fix = r(mean)
disp "嘉峪关市 2013 年 emp_total 修复："
disp " 原始值: " `jyg_raw'
disp " 替换值（2012/2014均值）: " `jyg_fix'
replace emp_total = `jyg_fix' if city == "嘉峪关市" & year == 2013
disp "(鹤壁市2013年不在2017-2019样本，跳过剔除）"
local n_before = _N
drop if city == "鹤壁市" & year == 2013
local n_after = _N
disp " 剔除前: `n_before' 条 → 剔除后: `n_after' 条"
disp "异常值处理后样本量: " _N
summ year, meanonly
disp "年份范围: " r(min) " ~ " r(max)
save "cleaned_data.dta", replace
disp "已保存：cleaned_data.dta"
|
3.4 计算全国层面基准（按年份）
use "cleaned_data.dta", clear
collapse (sum) emp_total emp_s1 emp_s2 emp_s3 emp_s4 emp_s5 emp_s6 emp_s7, by(year)
egen nat_pbs_total = rowtotal(emp_s1 emp_s2 emp_s3 emp_s4 emp_s5 emp_s6 emp_s7)
gen nat_pbs_share = nat_pbs_total / emp_total
rename emp_total E_nat
foreach s in 1 2 3 4 5 6 7 {
rename emp_s`s' nat_emp_s`s'
}
drop nat_pbs_total
save "national_benchmark.dta", replace
disp "已保存：national_benchmark.dta"
summ year, meanonly
disp "全国层面基准计算完成，年份数: " _N
|
3.5 测算专业化集聚指数 SP（公式4）
*- ------------------------------------------------------------
*- 四、测算专业化集聚指数 SP（公式4）
*-
*- SP_i = ∑_s(E_is/E_i) / ∑_s(E_s/E)
*- = [城市生产性服务业7行业就业 / 城市总就业] /
*- [全国生产性服务业7行业就业 / 全国总就业]
*-
*- 经济含义：
*- SP > 1：城市生产性服务业就业比重 > 全国平均水平（相对专业化集聚）
*- SP < 1：城市生产性服务业就业比重 < 全国平均水平
*- 数据处理：微信公众号 RStata
*- ------------------------------------------------------------
*- 此处代码需下载讲义材料查看～
|
3.6 测算多样化集聚指数 DV（公式5）
*- ------------------------------------------------------------
*- 五、测算多样化集聚指数 DV（公式5）
*-
*- 论文公式(5)的精确结构：
*-
*- DV_i = ∑_s (E_{i,s}/E_i) × [ 1/∑_{s'≠s}^n [ E_{i,s'} / (E_i - E_{i,s}) ]^2 ]
*- ÷ [ 1/∑_{s'≠s}^n [ E_{s'} / (E - E_s) ]^2 ]
*-
*- 关键：公式中的份额是"其余行业在剩余就业中的份额"，而非"占总就业份额"
*- - 城市层面：E_{i,s'} / (E_i - E_{i,s}) = 其余行业就业 / (城市总就业 - 行业s就业)
*- - 全国层面：E_{s'} / (E - E_s) = 其余行业就业 / (全国总就业 - 行业s就业)
*-
*- 使用 Mata 实现 DV 计算（逐城市×年份循环）
*- 数据处理：微信公众号 RStata
*- ------------------------------------------------------------
*- 此处代码需下载讲义材料查看～
|
Mata 代码详解（mata_dv.do）
Mata 是 Stata 的矩阵编程语言，适合处理复杂的循环计算。以下是 mata_dv.do 文件的完整代码：
*-------------------------------------------------------------
*- Mata 代码：计算多样化集聚指数 DV（公式5）
*- 这是一个独立的 do 脚本，由主脚本调用
*- 数据处理：微信公众号 RStata
*-------------------------------------------------------------
*- 此处代码需下载讲义材料查看～
|
Mata 代码关键点说明：
- 数据读取：st_view() 函数将 Stata 数据读取到 Mata 矩阵中
- 双重循环：外层循环遍历每个观测（城市×年份），内层循环遍历每个行业
- others 向量：使用 select() 函数创建排除当前行业s的其他行业索引向量
- 维度问题修复：使用 city_emp[1, others] 确保提取的是行向量而非列向量
- 结果写回：st_addvar() 创建新变量，st_view() 获取变量视图，然后赋值
3.7 缩尾处理与结果保存
use "sp_results.dta", clear
merge 1:1 year prov prov_code city cityid using "dv_results.dta", nogen keep(match)
sort year cityid
_pctile DV, p(1 99)
local dv_p1 = r(r1)
local dv_p99 = r(r99)
disp "DV Winsorize（1%）："
summ DV, meanonly
disp " 处理前范围: " r(min) " ~ " r(max)
disp " 1% 分位数: " `dv_p1'
disp " 99% 分位数: " `dv_p99'
count if DV < `dv_p1' | DV > `dv_p99'
local n_winsorized = r(N)
disp " 受影响观测数: `n_winsorized'"
gen 生产性服务业多样化集聚DV = DV
replace 生产性服务业多样化集聚DV = `dv_p1' if DV < `dv_p1'
replace 生产性服务业多样化集聚DV = `dv_p99' if DV > `dv_p99'
drop DV
summ 生产性服务业多样化集聚DV, detail
rename year 年份
rename prov 省
rename prov_code 省代码
rename city 市
rename cityid 市代码
rename SP 生产性服务业专业化集聚SP
label data "数据处理：微信公众号 RStata"
save "2017-2019年各城市生产性服务业集聚水平.dta", replace
disp "已保存：2017-2019年各城市生产性服务业集聚水平.dta"
|
点击这里跳转到 RStata 短书平台获取附件：名师讲堂｜使用 Stata 测算各城市生产性服务业集聚水平（Mata 版本）
评论