名师讲堂｜使用 Stata 测算各城市虚拟集聚程度

指标来源与背景

虚拟集聚度（Virtual Agglomeration Index，简称 Vag）是一种用于衡量城市间经济关联程度的指标。该指标由宋林等学者提出（参见《虚拟集聚与城市经济韧性》），其核心思想是：一个城市的经济发展不仅取决于自身条件，还受到周边城市的辐射影响。

数据来源

本教程使用的数据包括：

新增企业数据：1949～2023年各省市区县、行业新增企业数量统计.dta
注销企业数据：1970～2023年各年各省市区县、各行业注销公司数量面板数据.dta
行政区划数据：2021行政区划/市.shp（用于计算城市质心坐标）

数据处理基于工商注册信息，由 RStata 数据中心整理。

Stata 计算代码

环境设置与数据读取

*- 设置工作目录
cd "/Users/ac/Desktop/使用 Stata 测算各城市虚拟集聚程度"

clear
set more off
set maxvar 5000

第一步：计算城市质心坐标

使用 shp2dta 命令的 gencentroids() 选项直接提取城市的几何质心坐标：

*- 读取市级行政区划 shp 文件并转换为 Stata 格式
shp2dta using "2021行政区划/市.shp", ///
    data("市_data.dta") ///
    coor("市_coord.dta") ///
    genid(cityid) ///
    gencentroids(city) ///
    replace

*- 提取城市质心坐标
use "市_data.dta", clear
rename 市 cityname
keep cityid cityname x_city y_city
destring cityid, replace
save "city_centroids.dta", replace

说明：

gencentroids(city) 会生成 x_city 和 y_city 变量，表示城市质心的经纬度坐标
genid(cityid) 为每个城市生成唯一 ID

第二步：处理新增企业数据

*- 读取新增企业数据
use "1949～2023年各省市区县、行业新增企业数量统计.dta", clear
drop if missing(成立年份)
rename 新增注册企业数量 newcount

*- 按年份-省份-城市-行业汇总
collapse (sum) newcount, by(成立年份 省份 城市 行业门类)

rename 成立年份 year
rename 城市 city
rename 省份 province
rename 行业门类 industry
rename newcount cum_new

*- 计算累计新增企业数
sort province city industry year
by province city industry: gen cumsum_new = sum(cum_new)
save "new_summary.dta", replace

第三步：处理注销企业数据

此处代码需下载讲义材料查看～

第四步：合并计算存续企业数

*- 合并新增和注销数据
use "new_summary.dta", clear
merge 1:1 province city industry year using "quit_summary.dta", nogenerate

replace cumsum_new = 0 if missing(cumsum_new)
replace cumsum_quit = 0 if missing(cumsum_quit)

*- 计算存续企业数（累计新增 - 累计注销，下限为0）
gen surviving = cumsum_new - cumsum_quit
replace surviving = 0 if surviving < 0
save "firms_merged.dta", replace

第五步：汇总 IT 行业数据

local it_sector "信息传输、软件和信息技术服务业"

use "firms_merged.dta", clear
*- IT企业数 = IT行业存续企业数（非二元指示变量）
gen IT_firms = (industry == "`it_sector'") * surviving

*- 按年份-城市汇总
collapse (sum) total_surv = surviving (sum) IT_count = IT_firms, ///
    by(year city province)
save "city_year_data.dta", replace

第六步：计算区位熵

此处代码需下载讲义材料查看～

第七步：匹配城市坐标

use "city_year_lq.dta", clear
rename city cityname

*- 与城市质心坐标匹配
merge m:1 cityname using "city_centroids.dta", ///
    keepusing(x_city y_city) update nogenerate
drop if missing(x_city)
drop if year >= 9999
save "city_year_lq_coord.dta", replace

第八步：构建距离矩阵

*- 加载城市质心坐标
use "city_centroids.dta", clear
rename cityid cityid_i
rename cityname cityname_i
rename x_city lon_i
rename y_city lat_i

*- 使用 cross 命令生成所有城市对的笛卡尔积
cross using "city_centroids.dta"

*- 重命名 using 变量
rename cityid cityid_j
rename cityname cityname_j
rename x_city lon_j
rename y_city lat_j

*- 使用 geodist 命令计算球面距离
geodist lat_i lon_i lat_j lon_j, sphere gen(dist_km)

*- 距离倒数作为权重（自距离设为1）
gen weight = 1 / dist_km
replace weight = 1 if dist_km == 0 | missing(dist_km)
save "distance_matrix.dta", replace

第九步：计算虚拟集聚度

use "city_year_lq_coord.dta", clear
levelsof year, local(years)

*- 此处代码需下载讲义材料查看～

*- 整理输出格式
use "vag_final.dta", clear
drop if missing(city) | missing(year)
destring year national_IT national_total, replace
format year national_IT national_total total_surv IT_count %12.0g

rename year 年份
rename city 城市
rename total_surv 总企业数
rename IT_count IT企业数
rename national_IT 全国IT
rename national_total 全国总
rename nat_IT_ratio 全国IT集聚度
rename city_IT_ratio 城市IT集聚度
rename LQ 区位熵
rename Vag 虚拟集聚度

gsort 年份 城市
save "vag_results.dta", replace
export delimited using "vag_results.csv", replace

Stata 绘图代码

环境设置与配色

clear all
set more off

cd "/Users/ac/Desktop/使用 Stata 测算各城市虚拟集聚程度"
set scheme lightrstata

cap mkdir "picdir"

local caption_src = `"数据来源：RStata 数据中心，根据工商注册信息计算得到"'
local subti_txt "数据处理 & 绘制：微信公众号 RStata"
local mapdata_nhb "`c(pwd)'/map-data"
local mapdata "`c(pwd)'/map-data"

*- 配色方案
local c1 = "254 212 57"
local c2 = "112 154 225"

图1：各城市虚拟集聚度历年趋势

use "vag_results.dta", clear
keep if 年份 >= 1990

*- 标记直辖市
gen is_muni = 0
replace is_muni = 1 if inlist(城市, "北京市", "上海市", "天津市", "重庆市")

*- 获取所有城市列表
levelsof 城市, local(all_cities) sep(" ")

*- 构建绘图命令（其他城市用细灰线）
local cmd "tw"
foreach city in `all_cities' {
    local cmd `"`cmd' (line 虚拟集聚度 年份 if 城市 == `"`city'"', sort lc(gs12) lw(*0.2))"'
}

*- 高亮四个直辖市
此处代码需下载讲义材料查看～

gr export "picdir/01_各城市虚拟集聚度历年趋势.png", replace width(4800)

#| fig-cap: "各城市虚拟集聚度历年趋势（四个直辖市高亮显示）"
knitr::include_graphics("picdir/01_各城市虚拟集聚度历年趋势.png")

图2：2023年城市分布地图

use "vag_results.dta", clear
keep if 年份 == 2023

*- 创建城市匹配标签（去掉"市"后缀）
gen city_match = 城市
replace city_match = "北京" if 城市 == "北京市"
replace city_match = "上海" if 城市 == "上海市"
replace city_match = "天津" if 城市 == "天津市"
replace city_match = "重庆" if 城市 == "重庆市"
replace city_match = subinstr(city_match, "市", "", .)
save /tmp/city_data_2023.dta, replace

*- 加载地图数据并匹配
use "`mapdata'/chinacity2021mini_db.dta", clear
gen city_match = subinstr(市, "市", "", .)
merge 1:m city_match using /tmp/city_data_2023.dta, ///
    keepusing(虚拟集聚度) keep(master match) nogenerate

replace 虚拟集聚度 = -1000 if missing(虚拟集聚度)

*- 计算分位数断点
_pctile 虚拟集聚度 if 虚拟集聚度 != -1000, n(8)
local r1 = string(`r(r1)', "%9.3f")
local r2 = string(`r(r2)', "%9.3f")
local r3 = string(`r(r3)', "%9.3f")
local r4 = string(`r(r4)', "%9.3f")
local r5 = string(`r(r5)', "%9.3f")
local r6 = string(`r(r6)', "%9.3f")
local r7 = string(`r(r7)', "%9.3f")
sum 虚拟集聚度 if 虚拟集聚度 != -1000
local maxv = string(`=r(max)+0.1', "%9.3f")

*- 绘制带邻国的城市地图
此处代码需下载讲义材料查看～

gr export "picdir/02_2023年各城市虚拟集聚度分布地图_带邻国.png", replace width(4800)

#| fig-cap: "2023年各城市虚拟集聚度分布地图"
knitr::include_graphics("picdir/02_2023年各城市虚拟集聚度分布地图_带邻国.png")

点击这里跳转到 RStata 短书平台获取附件：名师讲堂｜使用 Stata 测算各城市虚拟集聚程度