. lxhuse smart_city, clear . merge m:1 id using https://file.lianxh.cn/data/s/smart_pt.dta // 匹配政策变量，生成政策时间变量 . drop _merge . replace pt=0 if pt==. . g did_1=0 . replace did_1=1 if pt==2012&year>=2011 // 政策时点提前1期 . g did_2=0 . replace did_2=1 if pt==2012&year>=2010 // 政策时点提前2期 . g did1=0 . replace did1=1 if pt==2012&year>=2013 // 政策时点滞后1期 . g did2=0 . replace did2=1 if pt==2012&year>=2014 // 政策时点滞后2期
. xtset id year . g lntgdp=lnrgdp*lnrgdp . gen lnww=log(ww) . gen lnso=log(so) . g lnrww=log(ww/pop) . xtreg lnrso did_1 lnrgdp lntgdp lninno lnurb lnopen lnss,fe . est store m1 . xtreg lnrso did_2 lnrgdp lntgdp lninno lnurb lnopen lnss,fe . est store m2 . xtreg lnrso did1 lnrgdp lntgdp lninno lnurb lnopen lnss,fe . est store m3 . xtreg lnrso did2 lnrgdp lntgdp lninno lnurb lnopen lnss,fe . est store m4
mat b = J(1000,1,0) // 系数矩阵 mat se = J(1000,1,0) // 标准误矩阵 mat p = J(1000,1,0) // P值矩阵
* 循环1000次 lxhget cao_treat.dta, replace lxhget cao_control.dta, replace forvalues i=1/1000{ use cao_treat.dta, clear qui xtset id year gen DID1=0 replace DID1=1 if pt==2005&year<2005 replace DID1=1 if pt==2009&year<2009 replace DID1=1 if pt==2012&year<2012 replace DID1=1 if pt==2013&year<2013 replace DID1=1 if pt==2014&year<2014 replace DID1=1 if pt==2015&year<2015 replace DID1=1 if pt==2016&year<2016 // 选取政策前年份 keep if DID1==1 sample 1, count by(id) // 根据id分组，每组随机抽取一个年份 keep id year // 得到所抽取样本的id编号 rename year policy_year save match_id_year.dta, replace // 另存id编号数据 merge 1:m id using cao_treat.dta // 与原数据匹配 qui xtset id year gen treatment = (_merge == 3) // 生成政策分组虚拟变量 gen period = (year >= policy_year) //生成政策时间虚拟变量 gen did0 = treatment*period append using cao_control.dta xtreg gdpr did0 invest consume export gov second agg innov i.year, fe cluster(id) * 将回归结果赋值到对应矩阵的对应位置 mat b[`i',1] = _b[did] // 系数矩阵 mat se[`i',1] = _se[did] // 标准误矩阵 mat p[`i',1] = 2*ttail(e(df_r), abs(_b[did]/_se[did])) // 计算P值并赋值于矩阵 } * 矩阵转化为向量 svmat b, names(coef) svmat se, names(se) svmat p, names(pvalue) * 删除空值并添加标签 drop if pvalue1 == . label var pvalue1 p值 label var coef1 估计系数 keep coef1 se1 pvalue1 gen tvalue = coef1/se1 // 计算t值 save placebo.dta, replace //关于p值，估计系数的文件，要用作画图 kdensity coef1, normal scheme(qleanmono) //绘制核密度分布图
mat b = J(1000,1,0) // 系数矩阵 mat se = J(1000,1,0) // 标准误矩阵 mat p = J(1000,1,0) // P值矩阵 lxhget new_area.dta, replace use new_area, clear xtset id year forvalues i=1/100 { use new_area, clear drop if treat==1 // 删除原来的处理组 drop treat gen random=1 sample 1, count gen treat=1 // 将新随机样本treat设为1 gen t=2005 // 将新随机样本政策时间设为2005年 save random_control_2005,replace
use new_area, clear drop if treat==1 // 删除原来的处理组 drop treat gen random=1 sample 1,count gen treat=1 // 将新随机样本treat设为1 gen t=2009 // 将新随机样本政策时间设为2009年 save random_control_2009,replace
use new_area, clear drop if treat==1 // 删除原来的处理组 drop treat gen random=1 sample 2,count gen treat=1 // 将新随机样本treat设为1 gen t=2012 // 将新随机样本政策时间设为2012年 save random_control_2012,replace
use new_area, clear drop if treat==1 // 删除原来的处理组 drop treat gen random=1 sample 2,count gen treat=1 // 将新随机样本treat设为1 gen t=2013 // 将新随机样本政策时间设为2013年 save random_control_2013,replace
use new_area, clear drop if treat==1 // 删除原来的处理组 drop treat gen random=1 sample 4,count gen treat=1 // 将新随机样本treat设为1 gen t=2014 // 将新随机样本政策时间设为2014年 save random_control_2014,replace
use new_area, clear drop if treat==1 // 删除原来的处理组 drop treat gen random=1 sample 5,count gen treat=1 // 将新随机样本treat设为1 gen t=2015 // 将新随机样本政策时间设为2015年 save random_control_2015,replace
use new_area, clear drop if treat==1 // 删除原来的处理组 drop treat gen random=1 sample 2,count gen treat=1 // 将新随机样本treat设为1 gen t=2016 // 将新随机样本政策时间设为2016年 save random_control_2016,replace
use random_control_2005,clear append using random_control_2009 append using random_control_2012 append using random_control_2013 append using random_control_2014 append using random_control_2015 append using random_control_2016 duplicates drop id year,force save random_control_all.dta,replace use new_area, clear drop treat gen random=0 gen start=0 append using random_control_all by id year,sort: gen test=_N drop if test==2 & random==0 by id,sort: egen t2=max(t) by id,sort: egen treat2=max(treat) gen DID0=1 if year>=t2&treat2==1 replace DID0=0 if DID0==. xtreg gdpr DID0 invest consume export gov second agg innov i.year, fe cluster(id) * 将回归结果赋值到对应矩阵的对应位置 mat b[`i',1] = _b[DID0] // 系数矩阵 mat se[`i',1] = _se[DID0] // 标准误矩阵 * 计算P值并赋值于矩阵 mat p[`i',1] = 2*ttail(e(df_r), abs(_b[DID0]/_se[DID0])) } * 矩阵转化为向量 svmat b, names(coef) svmat se, names(se) svmat p, names(pvalue) * 删除空值并添加标签 drop if pvalue1 == . label var pvalue1 p值 label var coef1 估计系数 keep coef1 se1 pvalue1 gen tvalue = coef1/se1 // 计算t值 save placebo.dta,replace // 关于p值，估计系数的文件，要用作画图 sum coef1 // 计算系数均值 kdensity coef1,normal scheme(qleanmono)
估计结果显示，此次变量 did 系数的均值为 -0.34，远小于基准回归结果 1.51，这表明国家级新区的政策效应表现出了明显的区位导向性，对设立国家级新区城市经济增长的带动效应最显著。
3. 同时随机化政策时间与处理组样本
以白俊红等 (2022) 为例，从所有样本中随机选取与原处理组相同数量的样本，并随机生成政策实施时间，构建城市与政策时间均随机的新处理组。基于此，重新估计基准回归模型，并随机重复 500 次实验，操作 Stata 代码如下：
clear set matsize 5000 mat b = J(500,1,0) mat se = J(500,1,0) mat p = J(500,1,0) lxhget inno_policy.dta, replace forvalues i=1/500{ use "inno_policy.dta" , clear xtset city year keep if year==2004 sample 71, count keep city save "atchcity.dta",replace merge 1:m city using "inno_policy.dta" gen groupnew=(_merge==3) // 生成伪处理组的虚拟变量 save "matchcity`i'.dta",replace
* 伪政策虚拟变量 use "inno_policy.dta",clear bsample 1, strata(city) keep year save "matchyear.dta", replace mkmat year, matrix(sampleyear) use "matchcity`i'.dta",replace xtset city year gen time = 0 foreach j of numlist 1/280 { replace time = 1 if (city == `j' & year >= sampleyear[`j',1]) } gen did=time*groupnew global xlist "lnagdp indust_stru finance ainternet market " xtreg entre_activation did $xlist i.year, fe robust
mat b[`i',1] = _b[did] mat se[`i',1] = _se[did] scalar df_r = e(N) - e(df_m) -1 mat p[`i',1] = 2*ttail(df_r,abs(_b[did]/_se[did])) } svmat b, names(coef) svmat se, names(se) svmat p, names(pvalue)
drop if pvalue1 == . label var pvalue1 p值 label var coef1 估计系数 save placebo.dta,replace // 保存数据
