今天给大家分享使用 Stata 测算地区产业专业化指标 的方法。该指标参考自杨本建、唐金汶（2022）《数字经济与区域产业布局》中的公式 (2)，其思想源于 Kalemli-Ozcan et al. (2003) 与 Du et al. (2022) 的 Krugman 式专业化指数 。
附件中提供了该参考文献的 PDF 文件（数字经济与区域产业布局.pdf），感兴趣的小伙伴可以阅读原文。
本文使用工商企业注册信息（已裁剪为测算所需的 9 个变量，置于 sample_data/ 文件夹），以 2000–2005 年为例，演示完整的 Stata 实现 。对应思路亦可参考同一项目的 R 语言版本讲义。
指标来源与计算过程 地区产业专业化指数（公式）
行业范围：制造业 C13–C43，剔除 C39 按国民经济行业分类（GB/T 4754）的制造业门类 ，取 2 位行业大类代码 C13–C43 ，并剔除 C39（计算机、通信和其他电子设备制造业） 。剔除 C39 的依据是论文第 248 页附录 1 的说明。
产业规模的代理变量 以企业的注册资本 作为行业规模的代理变量（主指标）；同时以企业数量 作为稳健性口径。两者计算逻辑完全一致，仅在汇总时替换聚合字段。
存续企业的界定（进入与退出） 参考李磊等（2023）的做法，按“进入—退出”口径统计每年各城市各行业的存续企业 ：
进入：以企业的成立年份作为进入年份；
退出：当经营状态含“注销/吊销/异常”时视为退出企业，退出年份取核准日期的年份；
存续：在年份 t 满足“成立年份 ≤ t，且（未退出 或 退出年份 > t）”的企业。
计算步骤概述
读取与清洗：读取各年工商注册数据，仅保留制造业、有效城市代码、必需变量，并界定进入/退出年份；
构建存续面板：对每个目标年份 t，筛选存活企业，按“城市 × 行业”汇总规模，得到城市×行业×年份存续企业规模面板；
计算专业化指数：对每一年，按公式计算各城市的 spec 指数（注册资本口径 + 企业数量口径）；
输出与可视化：导出结果 dta，并绘制趋势、分布与 Top 城市等图表。
详细讲解计算代码（Stata） 0. 路径与参数 首先设定工程目录、输出目录、数据目录，以及制造业行业代码（正则形式）与所需变量。
*- ---- 0. 路径与参数 ---- cd "/Users/ac/Desktop/使用 Stata 测算地区产业专业化指标/" global proj_dir = c(pwd) global out_dir "${proj_dir}/输出" global data_dir "${proj_dir}/sample_data" cap mkdir "${out_dir}" *- 目标年份 global file_years 2000 2001 2002 2003 2004 2005 global target_years 2000 2001 2002 2003 2004 2005 *- 制造业大类：C13-C43，剔除 C39（论文 248 页附录 1） *- 用正则判定行业大类落在 C13-C43 且非 C39（C39 不列入 alternation 即被排除） local mfg_pat "13|14|15|16|17|18|19|20|21|22|23|24|25|26|27|28|29|30|31|32|33|34|35|36|37|38|40|41|42|43" di "制造业大类正则（剔除 C39）：^C(`mfg_pat')$"
上述代码中：
脚本开头通过 cd “…” 切换到工程目录，再用 global proj_dir = c(pwd) 固定工程根目录，直接 do 运行即可，无需手动指定；
制造业大类通过正则表达式 ^C(13|14|…|43)$ 判定（C39 不列入 alternation 即被排除）；
因 Stata 变量名须为 ASCII，CSV 按列位置导入（varnames(nonames)），并 drop in 1 删除中文表头行，再用 rename (v1 v2 … v9) 改为英文变量名。
1. 读取并清洗单个年份文件 import delimited 读取每年 CSV，仅保留制造业且大类在 C13–C43（非 C39）、城市代码为 6 位数字，解析进入/退出年份，并剔除关键指标缺失或逻辑异常样本；最后把各年合并为企业级长表并缓存为 dta。
cap use "${out_dir}/firms_manufacturing_2000_2005.dta" , clear if _rc != 0 { tempfile all local first = 1 foreach yr of global file_years { di "读取 `yr'.csv ..." import delimited using "${data_dir}/`yr'.csv" , varnames(nonames) encoding("utf-8" ) clear drop in 1 rename (v1 v2 v3 v4 v5 v6 v7 v8 v9) (cap_reg cap_paid ind_sector ind_code op_status entry_year_raw approve_date city city_code) keep if ind_sector == "制造业" keep if regexm (ind_code, "^C(`mfg_pat')$" ) keep if regexm (city_code, "^[0-9][0-9][0-9][0-9][0-9][0-9]$" ) destring entry_year_raw, replace rename entry_year_raw entry_year gen approve_year = real (substr (approve_date, 1, 4)) destring cap_reg cap_paid, replace ignore("," ) gen exited = (regexm (op_status, "销" ) | regexm (op_status, "退" ) | regexm (op_status, "异常" )) gen exit_year = . replace exit_year = int(approve_year) if exited & !missing (approve_year) drop if missing (entry_year) | missing (cap_reg) | cap_reg <= 0 drop if !missing (exit_year) & exit_year < entry_year keep city_code city ind_code entry_year exit_year cap_reg cap_paid if `first' == 1 { save `all' , replace local first = 0 } else { append using `all' save `all' , replace } } di "清洗后制造业企业记录：`=_N' 条" save "${out_dir}/firms_manufacturing_2000_2005.dta" , replace } else { di "复用已存在的 firms_manufacturing_2000_2005.dta，共 `=_N' 条" }
要点：
keep if ind_sector == “制造业” 与 keep if regexm(ind_code, “^C(mfg_pat’)$”)` 完成行业筛选；
regexm(op_status, “销”) | regexm(op_status, “退”) | regexm(op_status, “异常”) 判定退出企业；
drop if missing(entry_year) | missing(cap_reg) | cap_reg <= 0 与 drop if !missing(exit_year) & exit_year < entry_year 剔除无效/异常样本。
2. 构建 城市×行业×年份 存续企业规模面板 对目标年份 t 筛选“存续企业”（成立年份 ≤ t，且未退出或退出年份 > t），按城市 × 行业汇总该年注册资本总额（output_reg）与企业数（n_firm），逐年份 append 后得到城市×行业×年份面板。
*- ---- 2. 构建 城市×行业×年份 存续企业规模面板 ---- *- 此处代码需要下载讲义材料查看～
3. 计算 Krugman 式地区产业专业化指数 对某一年的城市×行业规模数据，计算各城市专业化指数。核心逻辑：
求各城市该口径总规模 city_total，并计算各行业份额 share；
通过 cross 补全“城市 × 行业”全网格，该城市该行业无企业时份额记为 0；
对每个行业求全部城市份额之和 sum_share，则“其他城市平均份额”为 (sum_share - share)/(J-1)；
计算差的平方 (share - other_avg)^2，跨行业求和即得该城市 spec。
注册资本（output_reg）与企业数（n_firm）两个口径分别循环计算后合并。
*- ---- 3. 计算 Krugman 式地区产业专业化指数 ---- *- 此处代码需要下载讲义材料查看～
合并结果并预览 use "${out_dir}/spec_spec_reg.dta" , clear merge 1:1 year city_code using "${out_dir}/spec_spec_count.dta" , keep (3) nogeneratekeep year city_code city n_city spec_reg spec_countsort year city_codesave "${out_dir}/地区产业专业化指数_2000_2005.dta" , replace di "步骤 4/4：输出结果 ……" di "各年城市数(J)：" tabstat n_city, by (year) stat(mean ) format (%9.0f) nototaldi "2005 年产业专业化指数最高的 10 个城市：" use "${out_dir}/地区产业专业化指数_2000_2005.dta" , clear keep if year == 2005gsort -spec_reglist city spec_reg spec_count in 1/10, noobsdi "完成！结果已写入：${out_dir}"
运行结束后，输出/地区产业专业化指数_2000_2005.dta 包含变量：year, city_code, city, n_city, spec_reg, spec_count。
结果预览 下面用 R 的 haven 包读取 Stata 输出的 dta，展示 2005 年专业化指数最高的若干城市，便于核对（R 仅用于读表与制表，核心测算由上面的 Stata 代码完成）。
spec <- read_dta( file.path( out_dir, "地区产业专业化指数_2000_2005.dta" ) ) spec <- spec |> mutate( across( c ( spec_reg, spec_count, n_city) , as.numeric ) ) spec |> filter( year == 2005 ) |> slice_max( spec_reg, n = 15 ) |> select( year, city, n_city, spec_reg, spec_count) |> knitr:: kable( digits = 4 , caption = "2005 年产业专业化指数最高的 15 个城市（Stata 测算结果）" )
插图 数据可视化（Stata） 结果可视化由 02_可视化.do 完成，使用 Stata 的 twoway / graph box / graph hbar 绘制三张图表。
cd "/Users/ac/Desktop/使用 Stata 测算地区产业专业化指标/" global proj_dir = c(pwd )global out_dir "${proj_dir}/输出" cap mkdir "${out_dir}" use "${out_dir}/地区产业专业化指数_2000_2005.dta" , clear graph box spec_reg, over(year) ytitle("产业专业化指数 spec" ) title("各年份城市产业专业化指数分布 (2000-2005)" ) subtitle("箱线图展示全部城市的专业化指数分布及其演变" ) note ("数据处理 & 绘图：微信公众号 RStata" ) ylabel(,format (%6.1f)) graph export "${out_dir}/图2_专业化指数分布.png" , replace width(2000) height(1200)di "三张图表已导出至：${out_dir}"
图1：各城市平均专业化指数随年份变化趋势
图2：各年专业化指数分布（箱线图）
图3：2005 年专业化程度最高的 20 城市
点击这里跳转到 RStata 短书平台获取附件：名师讲堂｜使用 Stata 测算地区产业专业化指标
评论