论文解读与复现｜政治等级与区域经济发展：来自中国空间断点回归的证据（一）
引言
文章利用重庆晋升为中国省级直辖市这一政策冲击来考察城市的政治等级对经济发展的影响。重庆市升级为直辖市之后，获得了更高等级的行政权力，其结果可以从以下几个方面来看待：首先，升级后的重庆市获得了省级行政权力，可以采取完全独立于四川的行政政策。第二，重庆官员不再由四川省政府任命或评估。因此，它不需要像以前那样遵守四川省的政策，并可以享受的权力大幅增加。第三，关于财政事务，重庆不再需要与四川省政府分享财政收入，这使得重庆的财政收入留存率更高。
重庆升级为直辖市的历史事件使得本文能够利用空间断点回归设计来确定因果效应。重庆与四川接壤，生活在重庆-四川边境附近的人们有着相似的民族构成，生活在相似的地理条件下，有着相似的文化。更重要的是，他们长期属于同一个省。因此，在重庆升级时，边境两侧的大多数特征几乎是相同的。因此，四川的城镇是重庆的城镇的理想反事实，本文可以采用空间断点来比较重庆-四川边界的城镇发展情况。
制度背景
在中国，国家以下各级管辖区由四层政府组成：省一级(23个省、5个自治区、4个直辖市以及2个特别行政区)、地一级(约333个地级市)、县一级(约2853个县)和乡一级(约40497个镇)。
中国政治等级制度的一个显著特征是，下级政府完全从属于上级政府；因此，上级政府比下级政府拥有更大的权力。从宪法上来说，中国是一个单一的国家，而不是一个联邦国家，地区权力由中央政府授予。为了统治一个拥有大量人口和大陆面积的国家，中央政府使用嵌套的多层行政系统，其中上层政府管理和监督下层政府，从而将中央政府的控制扩展到所有级别的地区。因此，地方政府的决策权在很大程度上依赖于城市的等级制度。此外，下级领导人由其直接上级政府评估和任命，而不是由当地居民选举产生。在这种制度下，上级政府可以轻易地干预下级政府的政策决定，因此前者决定了后者的发展道路。最后，与行政系统相似，中国财政系统的特点是不同层级之间严格的垂直等级关系。更高一级的政府有权分配其下一级的财政收入。因此，中央政府决定与省政府分享收入，省政府反过来为各自的地级市制定规则。
重庆历史上一直是四川省的一部分。自元朝(公元1271-1368年)以来，重庆一直是四川省管辖的一个地区。1996年，四川省由24个地级市组成，包括成都、重庆、涪陵、万县和黔江。重庆当时是一个单列的地级市，其他23个都是地级市，四川省政府所在地为成都。1997年3月14日，中央政府决定将重庆与邻近的涪陵、万县、黔江等地合并，由此产生的单一区划成为重庆直辖市，其地位相当于一个省。重庆市人口3000万，分布在43个县和933个镇(不含县)。四川省的其余部分有8500万人口，分布在20个地级市、180个县和4155个镇。
根据中国官方的宣传，中央政府决定提升重庆有两个原因：首先，中央政府希望振兴重庆经济，使重庆成为推动中国西部地区增长的引擎。人们认为四川省政府不能在这样一个大省有效地提供公共产品。新成立的重庆市被认为能够制定更适合当地情况的政策，并由于其更高的行政地位而具有更大的促进经济发展的能力。第二，重庆市政府有望更好地协调三峡工程库区居民的安置。三峡工程始于1994年，位于重庆和湖北省之间的长江干流上。为了修建水库，100多万人进行了移民。
理论框架
重庆事的升级，可以从收益和成本两个角度进行考虑：
第一，升级之后，重庆可以采取完全独立于四川的政策方案，而且重庆官员不再由四川省政府任命或考核，而是直接由中央政府任命或考核；因此，他们不需要像以前那样遵守四川省的政策。由于法律上和事实上的权力大幅增加，重庆政府可以制定更贴近当地人民偏好的政策。特别是，重庆在许多经济事务中获得了更大的决策权，从批准投资项目(政府资助和私人资助)、土地分配和基础设施建设，到立法和执法。因此，重庆可以更容易地进行投资项目，对企业更有吸引力。
第二，从财政方面考虑，重庆获得了更多的公共服务责任，获得了更高的财政收入留存率。在升级之前，重庆的一部分(万县、黔江和涪陵)需要与四川省政府分享财政收入。升级之后，重庆能够直接与中央政府分享收入，不再与四川省政府分享收入。更高的税收留存率为地方政府提供了更高的激励，促使它们制定有利于发展的政策来增加税收基础。
第三，尽管在许多发展中国家，权力下放往往会导致更多的地方腐败，中国对地方领导人的奖惩制度可能会降低这种风险。此外，在升级之前，尽管重庆官员有晋升至国家政府的愿望，但他们在促进投资和吸引企业方面受到有限权力的限制。升级后权力下放，有助于鼓励官员促进经济发展。
最后，重庆的升级也可能会增加行政成本。然而，正如Reingewertz (2012)所建议的那样，规模经济减少所产生的成本可能只有在小城市激增时才有可能。由于重庆是一个人口超过3000万的大城市，它不太可能遭受规模经济的缩减。而且重庆案例中，省级单位只增加一个；因此，增加中央政府处理信息的控制范围的成本也很小。
数据来源
实证策略
计量模型设定
为了估计重庆升级为直辖市的处理效应，可以把四川省的乡镇作为重庆市的乡镇的反事情，基于两个地区的全样本进行DID估计。然而，作者发现，从全样本来看，重庆乡镇和四川乡镇的一些前变量特征（pre-treatment）并不平衡（imbalance），即1996年重庆和四川的灯光亮度、海拔、少数民族人口占比以及乡镇的数量均呈现出显著差异（见原文附表A2）。因此，简单的DID估计结果可能是有偏的。此外，作者进一步利用1992年到2013年的数据基于时间研究法进行估计，并不满足平行趋势假设（见原文附表A4）。在说明了，DID估计方法不可行之后，作者提出了空间断点回归模型：
其中， 为乡镇 的灯光亮度增长率，考虑到有些乡镇的亮度为0，因此定义 为：
为0-1虚拟变量，如果乡镇位于重庆，取值为1，若位于四川，取值为0。 为控制地理位置平滑功能的多项式，依据 Dell（2010），地理位置用乡镇质心的经度和维度表示。在这个模型中，捕获了重庆升级为直辖市的经济增长效应。回归过程中，文章聚类到县级层面。
断点回归的常用方法有两类：非参数局部线性回归和全局多项式回归。作者首先以距离四川重庆边界线30km为带宽进行了局部线性回归，并进而考虑不同带宽（20-100km）的稳健性（如图2）。此外，作者进一步报告了全局多项式回归的结果。
处理前的平衡性检验
断点回归的关键假设是平滑性假设，即除了处理效应，其他相关的变量，如分配变量和协变量，在四川-重庆的边界处都应该是连续的（不存在断点）。这一假设可以保证四川边界的乡镇是重庆边界乡镇的理想反事实对象。由于在1997年之前，重庆辖区隶属于重庆市，因此可以预期，边界处的前特征变量（pre-treatment characteristics）是平滑的。表2给出了比较的结果。根据表2，若将边界30km范围内的样本进行比较，重庆和四川的数据在1996年（处理之前）并没有显示出显著的差异性。图3以边界线为0点，绘制了距离边界线100KM范围以内的灯光亮度，结果进一步佐证了在1996年，重庆和四川边界的灯光亮度是没有断点的。
Stata 代码如下（相关数据可以在微信公众号 RStata 后台回复 0608 获取）：
* 表格 2 Panel A 代码
use data_si_chong, clear
reg mean1996 new if distance <= 30, cluster(county)
reg elevation new if distance <= 30, cluster(county)
reg slope new if distance <= 30, cluster(county)
use data_cencus2000.dta, clear
reg minorp new if distance <= 30, cluster(county)
use data_si_chong, clear
gen x_chengdu = 104.071208
gen y_chengdu = 30.665205
gen x_chongqing = 106.557078
gen y_chongqing = 29.569389
geodist ycenter xcenter y_chengdu x_chengdu, gen(dist_chengdu)
geodist ycenter xcenter y_chongq x_chongq, gen(dist_chongqing)
reg dist_chengdu new if distance <= 30, cluster(county)
reg dist_chongqing new if distance <= 30, cluster(county)
图 3 的 Stata 代码为：
use data_si_chong, clear
gen dist_g = distancereplace dist_g = -distance if new == 0
keep if distance <= 100 gen group = -100forvalues i = 1/20{ replace group = -100+(`i'*10) if dist_g >= -100 + (`i' * 10) & dist_g <- 90 + (`i' * 10) }
egen distance_mean = mean(dist_g), by(group)egen mean1996_mean = mean(mean1996), by(group)egen elevation_mean = mean(elevation), by(group)egen slope_mean = mean(slope), by(group)
reg mean1996 new dist_gpredict mean1996_ppredict se1,stdp gen mean1996_low = mean1996_p - 1.96 * se1gen mean1996_high = mean1996_p + 1.96 * se1by group, sort: egen float count = count(mean1996_p)
twoway ///scatter mean1996_mean distance_mean if /// distance <= 100, mcolor(black) || ///line mean1996_p mean1996_low mean1996_high dist_g /// if dist_g < 0 & dist_g >= -100, /// pstyle(p p3 p3) lcolor(black) lpattern(solid) || ///line mean1996_p mean1996_low mean1996_high dist_g /// if dist_g > 0 & dist_g <= 100, pstyle(p p3 p3) /// lcolor(black) lpattern(solid) /// ytitle(Light intensity of 1996) /// xtitle(Distance to cutoff) /// legend(off) xline(0,lpattern(dash)) /// text(1.7 50 'Chongqing',size(medlarge)) /// text(1.7 -50 'Sichuan',size(medlarge))gr export fig3.png, replace width(1200)接下来，作者进一步检验了四川和重庆边界的灯光强度的增长率是否在处理前存在平行趋势。为此，作者将（1）式中的因变量替换为：
并利用1992-2013年的数据，逐年进行回归，结果如图4所示。图4显示靠近边界处的处理和控制组在1997年以前的趋势基本是平衡的。
Stata 绘制代码如下：
use dynamic_1992-2013,clear
twoway ///
scatter beta t || ///
line beta t, lcolor(black) || ///
line beta_low t, lpattern(dash) lcolor(gray) || ///
line beta_high t, lpattern(dash) lcolor(gray) ///
xlabel(1992(2)2013,angle(45)) xtitle(Year) ///
yline(0, lpattern(dash)) ///
xline(1996, lpattern(dash)) ///
legend(off) ylabel(-1(0.5)2)
gr export fig4.png, replace width(1200)
主要结果
基本结论
作者首先依据（1）式进行了基准回归。回归的因变量为1996-2013年的灯光亮度增长率，回归方法分别为局部线性方法、局部二次项方法以及全部多项式回归，并且分别在带宽为30KM和带宽50KM处进行了回归。结果如表3所示，估计结果显示，重庆乡镇的灯光增长率显著高于四川。接下来，作者进一步以边界线为0点，绘制了距离边界线100KM范围以内的灯光亮度从1996-2013增长率的结果，结果显示存在明显断点。
表 3 的 Stata 代码如下：
use data_si_chong,clear
reg ln_diff_2013_1996 new xcenter ycenter if distance <= 30, cluster(county)reg ln_diff_2013_1996 new xcenter ycenter if distance <= 50, cluster(county)
gen x2 = xcenter^2gen y2 = ycenter^2gen xy = xcenter*ycenterreg ln_diff_2013_1996 new xcenter ycenter x2 y2 xy if distance <= 30, cluster(county)reg ln_diff_2013_1996 new xcenter ycenter x2 y2 xy if distance <= 50, cluster(county)
gen x3 = xcenter^3gen y3 = ycenter^3gen x2y = x2*ycentergen xy2 = xcenter*y2reg ln_diff_2013_1996 new x3 y3 x2y xy2 x2 y2 xy xcenter ycenter, cluster(county)
gen x4 = xcenter^4 gen x3y = xcenter^3*ycentergen x2y2 = xcenter^2*ycenter^2gen xy3 = xcenter*ycenter^3gen y4 = ycenter^4 reg ln_diff_2013_1996 new x4 x3y x2y2 xy3 y4 x3 y3 x2y xy2 x2 y2 xy xcenter ycenter, cluster(county)图 5 的 Stata 代码为：
use data_si_chong,clear
gen dist_g = distance
replace dist_g = -distance if new == 0
gen group = -100 if distance <= 100
forvalues i = 1(1)20{
replace group = -100 + (`i'*10) if dist_g >= -100 + (`i'*10) & dist_g < -90 + (`i'*10)
}
gen growthrate = ln(mean2013 + 0.01) - ln(mean1996 + 0.01)
egen growthrate_mean = mean(growthrate), by(group)
egen distance_mean = mean(dist_g), by(group)
qui reg growthrate new dist_g if distance <= 100
predict growthrate_p
predict se, stdp
gen growthrate_low = growthrate_p - 1.96 * se
gen growthrate_high = growthrate_p + 1.96 * se
by group, sort: egen float count = count(growthrate_p)
twoway ///
scatter growthrate_mean distance_mean ///
if distance <= 100, mcolor(black) || ///
line growthrate_p growthrate_low growthrate_high dist_g ///
if dist_g < 0 & dist_g >= -100, pstyle(p p3 p3) ///
lcolor(black) lpattern(solid) || ///
line growthrate_p growthrate_low growthrate_high dist_g if dist_g > 0 & dist_g <= 100, ///
pstyle(p p3 p3) lcolor(black) lpattern(solid) ///
ytitle(Night light growth rate 1996-2013) ///
xtitle(Distance to cutoff) legend(off) ///
xline(0, lpattern(dash)) ///
text(0.75 50 'Chongqing', size(medlarge)) ///
text(0.75 -50 'Sichuan', size(medlarge))
gr export fig5.png, replace width(1200)
稳健性检验
作者基于不同的带宽和模型设定进行了一系列的稳健性检验。图6中作者以10km的距离为增量，估计了20-100km带宽的结果。图A模型设定为线性；图B为二次多项式；图C为增加了一系列控制变量；图D为使用矩形核函数作为权重。结果均显示较为稳健。
Stata 代码为：
* 图 6 代码 panel A Linearuse data_si_chong,clear
matrix A = J(9, 2, 0)
forvalues i = 1/9{ qui reg ln_diff_2013_1996 new xcenter ycenter if distance <= (`i'+1)*10, cluster(county) matrix a`i' = e(b) matrix b`i' = e(V) scalar A`i'1 = a`i'[1,1] scalar A`i'2 = sqrt(b`i'[1,1]) matrix A[`i',1] = A`i'1 matrix A[`i',2] = A`i'2}
mat2txt, matrix(A) saving(matrixa) replaceinsheet using matrixa.txt,clearrename c1 coffrename c2 segen low = coff-1.96*segen high = coff+1.96*seegen v2 = group(v1)gen bw = (v2+1)*10twoway ///rcap low high bw, lcolor(cranberry) /// mcolor(midred) msymbol(O) mlabcolor(black) || ///scatter coff bw, mcolor(midblue) /// msymbol(O) mlabcolor(black) /// yline(0, lpattern(dash) lcolor(black)) /// legend(off) xtitle('Brandwidth') /// xlabel(20(10)100) ylabel(-0.5(0.5)2) /// ti('Panel A: Linear Polynomial') /// name(a, replace)
* 图 6 Panel B Quadraticuse data_si_chong,clear
gen x2 = xcenter^2gen y2 = ycenter^2gen xy = xcenter*ycenter
matrix B = J(9,2,0)
forvalues i = 1/9{ qui reg ln_diff_2013_1996 new xcenter ycenter x2 y2 xy if distance <= (`i'+1)*10, cluster(county) matrix a`i' = e(b) matrix b`i' = e(V) scalar B`i'1 = a`i'[1,1] scalar B`i'2 = sqrt(b`i'[1,1]) matrix B[`i',1] = B`i'1 matrix B[`i',2] = B`i'2}
mat2txt, matrix(B) saving(matrixb) replaceinsheet using matrixb.txt,clearrename c1 coffrename c2 segen low = coff-1.96*segen high = coff+1.96*seegen v2 = group(v1)gen bw = (v2+1)*10twoway ///rcap low high bw, lcolor(cranberry) /// mcolor(midred) msymbol(O) mlabcolor(black) || ///scatter coff bw, mcolor(midblue) /// msymbol(O) mlabcolor(black) /// yline(0, lpattern(dash) lcolor(black)) /// legend(off) xtitle('Brandwidth') /// xlabel(20(10)100) ylabel(-0.5(0.5)2) /// ti('Panel B: Quadratic Polynomial') /// name(b, replace)
* 图 6 Panel C Controluse data_si_chong, clear
matrix C = J(9,2,0)
forvalues i = 1/9{ qui reg ln_diff_2013_1996 new xcenter ycenter slope elevation seg1 seg2 seg3 seg4 if distance <= (`i'+1)*10, cluster(county) matrix a`i' = e(b) matrix b`i' = e(V) scalar C`i'1 = a`i'[1,1] scalar C`i'2 = sqrt(b`i'[1,1]) matrix C[`i',1] = C`i'1 matrix C[`i',2] = C`i'2}
mat2txt, matrix(C) saving(matrixc) replaceinsheet using matrixc.txt, clearrename c1 coffrename c2 segen low = coff-1.96*segen high = coff+1.96*seegen v2 = group(v1)gen bw = (v2+1)*10twoway ///rcap low high bw, lcolor(cranberry) /// mcolor(midred) msymbol(O) mlabcolor(black) || ///scatter coff bw, mcolor(midblue) msymbol(O) /// mlabcolor(black) yline(0, lpattern(dash) /// lcolor(black)) legend(off) /// xtitle('Brandwidth') xlabel(20(10)100) /// ylabel(-0.5(0.5)2) /// ti('Panel C: Additional Controls') /// name(c, replace)gr export fig6C.png, replace width(1200)
* 图 6 Panel D Weightuse data_si_chong,clearmatrix D = J(9,2,0)forvalue i = 1/9{ gen weight`i' = 1 - distance/(`i'+1)/10 if distance <= (`i'+1)*10 }
forvalues i = 1/9{ qui reg ln_diff_2013_1996 new xcenter ycenter if distance <= (`i'+1)*10 [weight = weight`i'],cluster(county) matrix a`i' = e(b) matrix b`i' = e(V) scalar D`i'1 = a`i'[1,1] scalar D`i'2 = sqrt(b`i'[1,1]) matrix D[`i',1] = D`i'1 matrix D[`i',2] = D`i'2}
mat2txt, matrix(D) saving(matrixd) replaceinsheet using matrixd.txt, clearrename c1 coffrename c2 segen low = coff-1.96*segen high = coff+1.96*seegen v2 = group(v1)gen bw = (v2+1)*10twoway rcap low high bw, lcolor(cranberry) /// mcolor(midred) msymbol(O) mlabcolor(black) || ///scatter coff bw, mcolor(midblue) msymbol(O) /// mlabcolor(black) yline(0, lpattern(dash) /// lcolor(black)) legend(off) /// xtitle('Brandwidth') xlabel(20(10)100) /// ylabel(-0.5(0.5)2) /// ti('Panel D: Triangular Kernel Weights') /// name(d, replace)
gr combine a b c d, row(2) gr export fig6.png, replace width(1200)由于内容过多，剩下的内容我们下次再更新～