********************************************************************************
/*						本dofile用于生成正文中的所有图表					*/
********************************************************************************

	clear all           //Refresh.
    cap log close       //Close any open log file.
    set more off        //Show the result all at once.
    set rmsg off
	
	
	** 设定工作路径
	gl root="E:\2024-00768_数据"
	cd $root
	gl dofiles="$root\2024-00768+程序代码"
	gl working_data="$root\data"
	gl logfiles="$root\2024-00768+日志文件"

	
	log using $logfiles/log_main, replace 
	log on 
	
	

********************************************************************************
/*									变量含义								  */
********************************************************************************
/*
	

	（1）县级数据中的变量，用于文章的主体分析。在实际分析中可能也会使用其差分、基期值，生成过程详见具体代码：
	Grain：县粮食产量（10万吨）。
	lnGDP_pc：按2000年价格指数平减的县实际人均GDP（元/人）的对数值。
	lnTempr：对数年均气温（摄氏度）。
	lnRainfall：对数年总降水量（毫米）。
	Share_PriGDP: 第一产业（农业）GDP占总GDP的百分比。
	Share_SecGDP：第二产业（工业）GDP占总GDP的百分比。
	lnGDP_agr_pc: 按2000年价格指数平减的县实际第一产业人均GDP（元/人）的对数值。
	lnGDP_ind_pc: 按2000年价格指数平减的县实际第二产业人均GDP（元/人）的对数值。
	lnGDP_serv_pc: 按2000年价格指数平减的县实际第三产业人均GDP（元/人）的对数值。
	lnGDP_nagr_pc: 按2000年价格指数平减的县实际非农产业（也即第二和第三产业加总）人均GDP（元/人）的对数值。
	lncstrct_land：县建设用地面积（平方公里）的对数值。
	lncstrct_land_new：县新增的建设用地面积（平方公里）的对数值。
	lnfarmland_dcrs：被转化为其他用途的耕地面积（平方公里）的对数值。
	lnslope_cstrct_land：县建设用地的平均坡度的对数值。
	lnnum_new_firms：县新注册企业的数量的对数值。
	share_ov1k_new_firms: 新注册企业中注册资本超过1000万元的比例。
	lnavgemp_old_firms：县域2003年以前注册的工业企业的平均就业人数。
	tfp_firms：工业企业的平均全要素生产率取对数。
	lntrsf_gen：县级政府所收到的一般性财政转移支付收入的对数值。
	lntrsf_spec：县级政府所收到的专项财政转移支付收入的对数值。
	lnpub_expd：县级政府的一般预算公共支出的对数值。
	lnpub_expd_agr：县级政府的农业相关公共支出的对数值。
	is_zcq：县是否属于粮食主产省的虚拟变量，是=1，否=0。
	is_phq：县是否属于粮食产销平衡省的虚拟变量，是=1，否=0
	is_zxq：县是否属于粮食主销省的虚拟变量，是=1，否=0。
	cqtype：县所属产销区的类型：1=主产区，2=平衡区，3=主销区。
	region_east: 是否属于东部地区，是=1，否=0。
	region_mid: 是否属于中部地区，是=1，否=0。
	region_west: 是否属于西部地区，是=1，否=0。
	region_noreast: 是否属于东北地区，是=1，否=0。
	lnemply: 就业总人口数的对数值。
	emply_shr_pri: 第一产业就业占总就业的百分比。
	emply_shr_sec: 第二产业就业占总就业的百分比。
	emply_shr_ter: 第三产业就业占总就业的百分比。
	GDP_pc: 县按2000年价格平减后的实际人均GDP（元/人）。
	pop: 县年末总人口（万人）。
	GDP: 县按2000年价格平减后的实际GDP（万元）。
	farmland_area：县域耕地总面积（平方公里）。
	Farmland_share：县域耕地面积占总面积的百分比。
	Sown_area：县农作物播种面积（万公顷）。
	shr_potprod_rice：县的潜在水稻总产出占全国水稻潜在总产出的比例。
	shr_potprod_wheat：县的潜在水稻总产出占全国小麦潜在总产出的比例。
	shr_potprod_maize：县的潜在水稻总产出占全国玉米潜在总产出的比例。
	shr_potprod_soybean：县的潜在水稻总产出占全国大豆潜在总产出的比例。
	shr_potprod_millet：县的潜在水稻总产出占全国谷子潜在总产出的比例。
	shr_potprod_sorghum：县的潜在水稻总产出占全国高粱潜在总产出的比例。
	dsim_Grain：按自然禀赋分配的潜在增产量（工具变量）。
	dsim_Grain_dwm：在构建时去除了小麦和谷子的按自然禀赋分配的潜在增产量（工具变量）。
	dsim_Grain_nosy：在构建时去除了大豆的按自然禀赋分配的潜在增产量（工具变量）。
	tarif_shock：2001年县层面按产业结构加总的平均进口关税。
	sim_Grain：按自然禀赋分配的潜在总产量。
	potyield_rice：县种植水稻的平均潜在可达单产。
	lnLight：县平均夜间灯光强度的对数值。
	lnPDen：县的人口密度对数值。
	confucious_temple：县孔庙的数量。
	soe_empshare：期初县域国有企业的就业占比。
	lnlandinc：县财政与土地出让相关的收入（国有土地使用权出让金收入）的对数值。
	lnlandexp：县财政与土地出让相关的支出（国有土地使用权出让金支出）的对数值。

	（2）省级数据中的变量（用于区域增产的描述性统计图1、表I1、核算总体增产成本表4）：
	grain: 省粮食产量，单位为万吨。
	gdppc: 省实际人均GDP
	mpop_per_county: 本省县级行政区的平均人口
	pop: 省总人口
	
	
*/

	

	
********************************************************************************
********************** 图1.粮食生产与主产区贡献 ********************************
********************************************************************************

	use $working_data\province_data.dta, clear
	
	
	* 变量处理
	//粮食产量单位从万吨换算为亿吨
	replace grain = grain/10000
	//计算全国总产出、主产区产出和占比
	bysort year: egen ttl_grain = total(grain)
	bysort year: egen zcq_grain = total(grain*is_zcq)
	gen grain_shr_zcq = (zcq_grain/ttl_grain)*100
	
 
	twoway (connected ttl_grain year if year>=1998&year<=2023, fcolor(none) lcolor(black) lwidth(medium) msymbol(triangle) msize(small) mcolor(black) lpattern(dash) lwidth(thin)) || (connected zcq_grain year if year>=1998&year<=2023, fcolor(none) lcolor(black) lwidth(medium) msymbol(circle) msize(small) mcolor(black) lpattern(solid) lwidth(thin)) || (connected grain_shr_zcq year if year>=1998&year<=2023, yaxis(2) fcolor(none) lcolor(black) lwidth(medium) msymbol(+) msize(small) mcolor(black) lpattern(dash_dot) lwidth(medium)) ||, xtitle("年份",size(medlarge) alignment(top) margin(small)) xlabel(1998 2003 2016 2023,labsize(small)) xline(2003,lcolor(black) lstyle(refline) lpattern(dash)) xline(2016,lcolor(black) lpattern(dash)) ytitle("粮食产量(亿吨)",axis(1) size(medlarge) alignment(top) margin(medsmall)) ytitle("粮食产量贡献(%)",axis(2) size(medlarge) alignment(top) margin(medsmall)) ylabel(,labsize(small) axis(1)) ylabel(70(2)80,labsize(small) axis(2)) scheme(s1mono) legend(label(1 "全国粮食产量") label(2 "主产区粮食产量") label(3 "主产区粮食产量贡献") size(medsmall) ring(0) rows(3) position(11)) 
	graph export "$logfiles/figure_1.png",replace width(2400) height(1800)
	
	

	
	

********************************************************************************
********************* 图2.县域粮食产量与人均GDP的关系  *************************
********************************************************************************

	use $working_data\county_panel.dta ,clear
	

	* 人均单位转换：元->万元
	replace GDP_pc = GDP_pc/10000
	
	* 删去极端值
	winsor2 Grain GDP_pc, suffix(_tr) trim cuts(1 99) by(year)

	
	
	//2000
	reg GDP_pc_tr Grain_tr if year==2000, vce(robust)
	local beta : display %6.3f _b[Grain_tr]
	local se : display %6.3f _se[Grain_tr]	
	local p_val : display %6.3f 2*ttail(e(df_r), abs(_b[Grain_tr]/_se[Grain_tr]))
	twoway scatter GDP_pc_tr Grain_tr if year==2000, msize(vsmall) mcolor(gs10) || lfit GDP_pc_tr Grain_tr if year==2000,lcolor(black) ||, scheme(s1mono) xtitle("粮食产量(十万吨)") ytitle("人均GDP(万元)",margin(medsmall)) legend(off) note("coef = `beta', SE = `se', p = `p_val'",ring(0) position(4) box bcolor(bg) lcolor(black)) name(gdpgrain_00,replace) subtitle("2000年")
	
	
	
	//2003
	reg GDP_pc_tr Grain_tr if year==2003, vce(robust)
	local beta : display %06.3f _b[Grain_tr]
	local se : display %6.3f _se[Grain_tr]	
	local p_val : display %6.3f 2*ttail(e(df_r), abs(_b[Grain_tr]/_se[Grain_tr]))
	twoway scatter GDP_pc_tr Grain_tr if year==2003, msize(vsmall) mcolor(gs10) || lfit GDP_pc_tr Grain_tr if year==2003,lcolor(black) ||, scheme(s1mono) xtitle("粮食产量(十万吨)") ytitle("人均GDP(万元)",margin(medsmall)) legend(off) note("coef = `beta', SE = `se', p = `p_val'",ring(0) position(4) box bcolor(bg) lcolor(black)) name(gdpgrain_03,replace) subtitle("2003年")
	
	
	
	//2016
	reg GDP_pc_tr Grain_tr if year==2016, vce(robust)
	local beta : display %06.3f _b[Grain_tr]
	local se : display %6.3f _se[Grain_tr]	
	local p_val : display %6.3f 2*ttail(e(df_r), abs(_b[Grain_tr]/_se[Grain_tr]))
	twoway scatter GDP_pc_tr Grain_tr if year==2016, msize(vsmall) mcolor(gs10) || lfit GDP_pc_tr Grain_tr if year==2016,lcolor(black) ||, scheme(s1mono) xtitle("粮食产量(十万吨)") ytitle("人均GDP(万元)",margin(medsmall)) legend(off) note("coef = `beta', SE = `se', p = `p_val'",ring(0) position(4) box bcolor(bg) lcolor(black)) name(gdpgrain_16,replace) subtitle("2016年")

	
	
	//2020
	reg GDP_pc_tr Grain_tr if year==2020, vce(robust)
	local beta : display %06.3f _b[Grain_tr]
	local se : display %6.3f _se[Grain_tr]	
	local p_val : display %6.3f 2*ttail(e(df_r), abs(_b[Grain_tr]/_se[Grain_tr]))
	twoway scatter GDP_pc_tr Grain_tr if year==2020, msize(vsmall) mcolor(gs10) || lfit GDP_pc_tr Grain_tr if year==2020,lcolor(black) ||, scheme(s1mono) xtitle("粮食产量(十万吨)") ytitle("人均GDP(万元)",margin(medsmall)) legend(off) note("coef = `beta', SE = `se', p = `p_val'",ring(0) position(4) box bcolor(bg) lcolor(black)) name(gdpgrain_20,replace) subtitle("2020年")
	

	* 合并输出
	graph combine gdpgrain_00 gdpgrain_03 gdpgrain_16 gdpgrain_20, scheme(s1mono) xsize(30) ysize(25)
	graph export "$logfiles/figure_2.png", replace height(2500) width(3000)


	
	
********************************************************************************
****************** 表1 .粮食增产与县域经济发展：基准回归 ************************
********************************************************************************
	
	
	use $working_data\county_panel.dta ,clear
	
	capture mkdir "$logfiles/table_1"  //创建单独文件夹输出

						***--- 变量处理 ---***
	
	* 生成变量的基期(2003)值
	gen is_yr2003 = 1 if year==2003
	foreach v of varlist Grain lnGDP_pc lnTempr lnRainfall Share_PriGDP Share_SecGDP{
		
		bysort county_id (year): egen `v'_03 = mean(`v'*is_yr2003)
		
	}
	
	* 生成相对基期的差分
	foreach v of varlist Grain lnGDP_pc{
		
		gen d_`v' = `v' - `v'_03
		
	}
	
	
	
						***--- 工具变量估计 ---***
			//注:本部分将输出两个结果文件，"table1_2sls"为表1的2SLS回归结果,"table1_fststage"为表1的第一阶段回归结果

	est clear
	*(1). 无控制变量
	** 2SLS
	ivreghdfe d_lnGDP_pc (d_Grain = dsim_Grain) if year==2016, absorb(city_id) vce(cluster city_id)
	outreg2 using "$logfiles/table_1/table1_2sls", title("2SLS Baseline") append nocons nor2 dec(3) excel addstat(K-P Wald F, e(rkf)) adec(3) addtext(City FE,Y)
	** First-stage
	reghdfe d_Grain dsim_Grain if e(sample)==1, absorb(city_id) vce(cluster city_id)
	outreg2 using "$logfiles/table_1/table1_fststage", title("First Stage") append nocons dec(3) keep(dsim_Grain) excel	
	
	
	*(2).增加控制基期人均GDP
	** 2SLS
	ivreghdfe d_lnGDP_pc lnGDP_pc_03 (d_Grain = dsim_Grain) if year==2016, absorb(city_id) vce(cluster city_id)
	outreg2 using "$logfiles/table_1/table1_2sls", title("2SLS Baseline") append nocons nor2 dec(3) excel addstat(K-P Wald F, e(rkf)) adec(3) addtext(City FE,Y)
	** First-stage
	reghdfe d_Grain dsim_Grain lnGDP_pc_03 if e(sample)==1, absorb(city_id) vce(cluster city_id)
	outreg2 using "$logfiles/table_1/table1_fststage", title("First Stage") append nocons dec(3) keep(dsim_Grain) excel	
	
	
	*(3).增加控制期初年均温
	** 2SLS
	ivreghdfe d_lnGDP_pc lnGDP_pc_03 lnTempr_03 (d_Grain = dsim_Grain) if year==2016, absorb(city_id) vce(cluster city_id)
	outreg2 using "$logfiles/table_1/table1_2sls", title("2SLS Baseline") append nocons nor2 dec(3) excel addstat(K-P Wald F, e(rkf)) adec(3) addtext(City FE,Y)
	** First-stage
	reghdfe d_Grain dsim_Grain lnGDP_pc_03 lnTempr_03 if e(sample)==1, absorb(city_id) vce(cluster city_id)
	outreg2 using "$logfiles/table_1/table1_fststage", title("First Stage") append nocons dec(3) keep(dsim_Grain) excel	
	
	
	*(4).增加控制期初年降水
	** 2SLS
	ivreghdfe d_lnGDP_pc lnGDP_pc_03 lnTempr_03 lnRainfall_03 (d_Grain = dsim_Grain) if year==2016, absorb(city_id) vce(cluster city_id)
	outreg2 using "$logfiles/table_1/table1_2sls", title("2SLS Baseline") append nocons nor2 dec(3) excel addstat(K-P Wald F, e(rkf)) adec(3) addtext(City FE,Y)
	** First-stage
	reghdfe d_Grain dsim_Grain lnGDP_pc_03 lnTempr_03 lnRainfall_03 if e(sample)==1, absorb(city_id) vce(cluster city_id)
	outreg2 using "$logfiles/table_1/table1_fststage", title("First Stage") append nocons dec(3) keep(dsim_Grain) excel	
	
	
	*(5).增加控制期初农业GDP占比
	** 2SLS
	ivreghdfe d_lnGDP_pc lnGDP_pc_03 lnTempr_03 lnRainfall_03 Share_PriGDP_03 (d_Grain = dsim_Grain) if year==2016, absorb(city_id) vce(cluster city_id)
	outreg2 using "$logfiles/table_1/table1_2sls", title("2SLS Baseline") append nocons nor2 dec(3) excel addstat(K-P Wald F, e(rkf)) adec(3) addtext(City FE,Y)
	** First-stage
	reghdfe d_Grain dsim_Grain lnGDP_pc_03 lnTempr_03 lnRainfall_03 Share_PriGDP_03 if e(sample)==1, absorb(city_id) vce(cluster city_id)
	outreg2 using "$logfiles/table_1/table1_fststage", title("First Stage") append nocons dec(3) keep(dsim_Grain) excel	
	
	
	*(6).增加控制期初农业GDP占比
	** 2SLS
	ivreghdfe d_lnGDP_pc lnGDP_pc_03 lnTempr_03 lnRainfall_03 Share_PriGDP_03 Share_SecGDP_03 (d_Grain = dsim_Grain) if year==2016, absorb(city_id) vce(cluster city_id)
	outreg2 using "$logfiles/table_1/table1_2sls", title("2SLS Baseline") append nocons nor2 dec(3) excel addstat(K-P Wald F, e(rkf)) adec(3) addtext(City FE,Y)
	** First-stage
	reghdfe d_Grain dsim_Grain lnGDP_pc_03 lnTempr_03 lnRainfall_03 Share_PriGDP_03 Share_SecGDP_03 if e(sample)==1, absorb(city_id) vce(cluster city_id)
	outreg2 using "$logfiles/table_1/table1_fststage", title("First Stage") append nocons dec(3) keep(dsim_Grain) excel	
	
	
	
********************************************************************************
******************** 表2 .粮食增产与县域经济发展：分部门影响  *******************
********************************************************************************
	
	use $working_data\county_panel.dta ,clear

						***--- 变量处理 ---***
	
	* 生成变量的基期(2003)值
	gen is_yr2003 = 1 if year==2003
	
	foreach v of varlist Grain lnGDP_pc lnGDP_agr_pc lnGDP_nagr_pc lnGDP_ind_pc lnGDP_serv_pc{
		
		bysort county_id (year): egen `v'_03 = mean(`v'*is_yr2003)
		
	}
	
	* 生成相对基期的差分
	foreach v of varlist Grain lnGDP_agr_pc lnGDP_nagr_pc lnGDP_ind_pc lnGDP_serv_pc{
		
		gen d_`v' = `v' - `v'_03
		
	}
	
	
	gl control lnGDP_pc_03
	
	est clear
	foreach y of varlist d_lnGDP_agr_pc d_lnGDP_nagr_pc d_lnGDP_ind_pc d_lnGDP_serv_pc{
		
		ivreghdfe `y' $control (d_Grain = dsim_Grain) if year==2016, absorb(city_id) vce(cluster city_id)
		outreg2 using "$logfiles/table_2", title("Effects by Sectors") append nocons nor2 dec(3) excel addstat(K-P Wald F, e(rkf)) adec(3) keep(d_Grain) addtext(Controls,Y,City FE,Y)
		
	}
		

		
		
	
*******************************************************************************
****************** 图3 .粮食增产对经济发展的各年份累积影响 *********************
*******************************************************************************


	use $working_data\county_panel.dta ,clear
	capture mkdir "$logfiles/figure_3"

						***--- 变量处理 ---***
	
	* 生成变量的基期(2003)值
	gen is_yr2003 = 1 if year==2003
	foreach v of varlist Grain lnGDP_pc{
		
		bysort county_id (year): egen `v'_03 = mean(`v'*is_yr2003)
		
	}
	
	* 生成相对基期的差分
	foreach v of varlist Grain lnGDP_pc{
		
		gen d_`v' = `v' - `v'_03
		
	}
	
	
	* 以1998-2016逐年作为差分末期分别估计
	/*2003年的位置由于基期=末期，
	差分的自变量和因变量都为0，回归无意义，会出现红叉，
	后续会将其设为0，不影响结果。*/
	keep if year>=1998&year<=2016
	gl control lnGDP_pc_03
	statsby _b[d_Grain] _se[d_Grain],by(year) saving("$logfiles/figure_3/est_fig3",replace): ivreghdfe d_lnGDP_pc $control (d_Grain = dsim_Grain), absorb(city_id) cluster(city_id)

	
	
	* 结果输出图片
	preserve
	use "$logfiles/figure_3/est_fig3.dta",clear
	rename _stat_1 est_coef
	rename _stat_2 est_se
	* 2003年基期=末期，设为0作基准
	replace est_coef = 0 if year==2003
	replace est_se = 0 if year==2003
	//95% CI
	gen lower_95ci = est_coef - 1.96*est_se
	gen upper_95ci = est_coef + 1.96*est_se
	
	twoway connected est_coef year,msize(small) || rline lower_95ci upper_95ci year, lpattern(dash) ||, scheme(s1mono) legend(off) yline(0,lpattern(refline)) xline(2003) xlabel(1998(1)2016,labsize(vsmall)) ylabel(-0.2(0.1)0.2)  ytitle("估计系数", margin(small) size(medium)) xtitle("年份",size(medium))
	graph export "$logfiles/figure_3/figure_3.png",replace width(2400) height(1800)
	
	restore
		
	
	
*******************************************************************************	
****************** 表3 .粮食增产与县域经济发展：机制分析 ***********************
*******************************************************************************


	use $working_data\county_panel.dta ,clear
	capture mkdir "$logfiles/table_3"

					***--- 变量处理 ---***
						
	gl y_mechanism lncstrct_land lncstrct_land_new lnfarmland_dcrs lnslope_cstrct_land lnnum_new_firms share_ov1k_new_firms lnavgemp_old_firms tfp_firms lntrsf_gen lntrsf_spec lnpub_expd lnpub_expd_agr
	
	* 生成变量的基期(2003)值
	gen is_yr2003 = 1 if year==2003
	foreach v of varlist Grain lnGDP_pc   $y_mechanism{
		
		bysort county_id (year): egen `v'_03 = mean(`v'*is_yr2003)
		
	}
	
	* 生成相对基期的差分
	foreach v of varlist Grain $y_mechanism{
		
		gen d_`v' = `v' - `v'_03
		
	}


	
						***--- 回归结果 ---***
						
	//注: 本部分将输出三个结果文件，分别对应表3的Panel A, B, C.
	gl control lnGDP_pc_03
	
	* Panel A. 土地配置
	* (1). 建设用地面积
	
	ivreghdfe d_lncstrct_land $control (d_Grain = dsim_Grain) if year==2016, absorb(city_id) cluster(city_id)
	
	outreg2 using "$logfiles/table_3/table_3_A", title("Mechanism_A: Land Allocation") ctitle(Construction Land Area) append nocons nor2 dec(3) excel addstat(K-P Wald F, e(rkf)) adec(3) keep(d_Grain) addtext(Controls,Y,City FE,Y)

	
	* (2).新增建设用地面积
	
	ivreghdfe d_lncstrct_land_new $control (d_Grain = dsim_Grain) if year==2016, absorb(city_id) cluster(city_id)
	
	outreg2 using "$logfiles/table_3/table_3_A", title("Mechanism_A: Land Allocation") ctitle(Increased Construction Land Area) append nocons nor2 dec(3) excel addstat(K-P Wald F, e(rkf)) adec(3) keep(d_Grain) addtext(Controls,Y,City FE,Y)

	
	* (3).被转化的耕地面积
	
	ivreghdfe d_lnfarmland_dcrs $control (d_Grain = dsim_Grain) if year==2016, absorb(city_id) cluster(city_id)
	
	outreg2 using "$logfiles/table_3/table_3_A", title("Mechanism_A: Land Allocation") ctitle(Decreased Farmland Area) append nocons nor2 dec(3) excel addstat(K-P Wald F, e(rkf)) adec(3) keep(d_Grain) addtext(Controls,Y,City FE,Y)

	
	* (4).建设用地平均坡度
	
	ivreghdfe d_lnslope_cstrct_land $control (d_Grain = dsim_Grain) if year==2016, absorb(city_id) cluster(city_id)
		
	outreg2 using "$logfiles/table_3/table_3_A", title("Mechanism_A: Land Allocation") ctitle(Slope of Construction Land) append nocons nor2 dec(3) excel addstat(K-P Wald F, e(rkf)) adec(3) keep(d_Grain) addtext(Controls,Y,City FE,Y)
	
	
	
	
	* Panel B. 企业表现
	* (1).新注册企业数量
	
	ivreghdfe d_lnnum_new_firms $control (d_Grain = dsim_Grain) if year==2016, absorb(city_id) cluster(city_id)
	
	outreg2 using "$logfiles/table_3/table_3_B", title("Mechanism_B: Firm Performance") ctitle(Number of New Registration) append nocons nor2 dec(3) excel addstat(K-P Wald F, e(rkf)) adec(3) keep(d_Grain) addtext(Controls,Y,City FE,Y)
	
	
	
	* (2).新注册企业中大企业占比(注册资本>1000万元)
	
	ivreghdfe d_share_ov1k_new_firms $control (d_Grain = dsim_Grain) if year==2016, absorb(city_id) cluster(city_id)
	
	outreg2 using "$logfiles/table_3/table_3_B", title("Mechanism_B: Firm Performance") ctitle(Share of Large Firms) append nocons nor2 dec(3) excel addstat(K-P Wald F, e(rkf)) adec(3) keep(d_Grain) addtext(Controls,Y,City FE,Y)
	
	
	* (3).在位企业平均就业人数
	//受限于数据可得性，工企相关结果差分末期为2013
	ivreghdfe d_lnavgemp_old_firms $control (d_Grain = dsim_Grain) if year==2013, absorb(city_id) cluster(city_id)
	
	outreg2 using "$logfiles/table_3/table_3_B", title("Mechanism_B: Firm Performance") ctitle(Employment of Incumbent Firms) append nocons nor2 dec(3) excel addstat(K-P Wald F, e(rkf)) adec(3) keep(d_Grain) addtext(Controls,Y,City FE,Y)
	
	
	
	* (4).工业企业平均全要素生产率
	//受限于数据可得性，工企相关结果差分末期为2013
	ivreghdfe d_tfp_firms $control (d_Grain = dsim_Grain) if year==2013, absorb(city_id) cluster(city_id)  
	
	outreg2 using "$logfiles/table_3/table_3_B", title("Mechanism_B: Firm Performance") ctitle(TFP) append nocons nor2 dec(3) excel addstat(K-P Wald F, e(rkf)) adec(3) keep(d_Grain) addtext(Controls,Y,City FE,Y)
	
	
	
	
	
	* Panel C.财政资源与公共支出
	//受限于数据可得性，农业支出差分末期为2006，其余财政相关结果差分末期为2007
	
	* (1).一般性转移支付收入
	
	ivreghdfe d_lntrsf_gen $control (d_Grain = dsim_Grain) if year==2007, absorb(city_id) cluster(city_id)
	
	outreg2 using "$logfiles/table_3/table_3_C", title("Mechanism_C: Fiscal Capacity and Public Expenditures") ctitle(General Transfers) append nocons nor2 dec(3) excel addstat(K-P Wald F, e(rkf)) adec(3) keep(d_Grain) addtext(Controls,Y,City FE,Y)
	
	
	* (2).专项转移支付收入
	
	ivreghdfe d_lntrsf_spec $control (d_Grain = dsim_Grain) if year==2007, absorb(city_id) cluster(city_id)
	
	outreg2 using "$logfiles/table_3/table_3_C", title("Mechanism_C: Fiscal Capacity and Public Expenditures") ctitle(Special Transfers) append nocons nor2 dec(3) excel addstat(K-P Wald F, e(rkf)) adec(3) keep(d_Grain) addtext(Controls,Y,City FE,Y)

	
	* (3). 公共预算支出
	
	ivreghdfe d_lnpub_expd $control (d_Grain = dsim_Grain) if year==2007, absorb(city_id) cluster(city_id)
	
	outreg2 using "$logfiles/table_3/table_3_C", title("Mechanism_C: Fiscal Capacity and Public Expenditures") ctitle(Public Expenditures) append nocons nor2 dec(3) excel addstat(K-P Wald F, e(rkf)) adec(3) keep(d_Grain) addtext(Controls,Y,City FE,Y)
	
	
	* (4). 农业相关支出
	
	ivreghdfe d_lnpub_expd_agr $control (d_Grain = dsim_Grain) if year==2006, absorb(city_id) cluster(city_id)
	
	outreg2 using "$logfiles/table_3/table_3_C", title("Mechanism_C: Fiscal Capacity and Public Expenditures") ctitle(Public Agricultural Expenditures) append nocons nor2 dec(3) excel addstat(K-P Wald F, e(rkf)) adec(3) keep(d_Grain) addtext(Controls,Y,City FE,Y)
	
	
	
	
********************************************************************************
************************ 表4 .粮食增产的经济机会成本估算 ************************
********************************************************************************
	
	
	
						*1.用县样本估算
	
	use $working_data/county_panel.dta ,clear
	capture mkdir $logfiles/table_4
	
	
	
	** (1).获取基准估计系数，以及分产区的估计系数
	
	* 生成变量的基期(2003)值
	gen is_yr2003 = 1 if year==2003
	foreach v of varlist Grain lnGDP_pc{
		
		bysort county_id (year): egen `v'_03 = mean(`v'*is_yr2003)
		
	}
	
	* 生成相对基期的差分
	foreach v of varlist Grain lnGDP_pc{
		
		gen d_`v' = `v' - `v'_03
		
	}

	* 基准估计 (保存系数保留小数点后三位，不会影响结论)
	ivreghdfe d_lnGDP_pc lnGDP_pc_03 (d_Grain = dsim_Grain) if year==2016, absorb(city_id) vce(cluster city_id)
	gen coef_base = round(_b[d_Grain],0.001)
	gl ebase = round(_b[d_Grain],0.001)
	
	* 分产区
	gen coef_cq = .
	ivreghdfe d_lnGDP_pc lnGDP_pc_03 (d_Grain = dsim_Grain) if year==2016&is_zcq==1, absorb(city_id) vce(cluster city_id)
	replace coef_cq = round(_b[d_Grain],0.001) if is_zcq==1
	gl ezcq = round(_b[d_Grain],0.001)
	
	ivreghdfe d_lnGDP_pc lnGDP_pc_03 (d_Grain = dsim_Grain) if year==2016&is_zxq==1, absorb(city_id) vce(cluster city_id)
	replace coef_cq = round(_b[d_Grain],0.001) if is_zxq==1
	gl ezxq = round(_b[d_Grain],0.001)

	ivreghdfe d_lnGDP_pc lnGDP_pc_03 (d_Grain = dsim_Grain) if year==2016&is_phq==1, absorb(city_id) vce(cluster city_id)
	replace coef_cq = round(_b[d_Grain],0.001) if is_phq==1
	gl ephq = round(_b[d_Grain],0.001)
	
	
	** (2).计算增产成本
	
	* 计算样本期间各县增产量，平均人口、人均GDP以及GDP
	gen is_yr2016 = 1 if year==2016
	foreach v of varlist GDP_pc pop{
		
		bysort county_id (year): egen `v'_03 = mean(`v'*is_yr2003)
		bysort county_id (year): egen `v'_16 = mean(`v'*is_yr2016)
		
	}
	
	
	bysort county_id (year): egen Grain_16 = mean(Grain*is_yr2016)
	
	
	gen dgrain_0316 = Grain_16 - Grain_03
	gen mgdppc_0316 = (GDP_pc_03 + GDP_pc_16)/2
	gen mpop_0316 = (pop_03 + pop_16)/2
	gen mgdp_0316 = mgdppc_0316*mpop_0316
	

	drop if dgrain_0316==.|mgdppc_0316==.|mpop_0316==.
	keep if year==2016
	
	
	
	* A.按照基准系数估计成本
	* i.各县增产成本(GDP)=产量变化*估计系数*人均GDP*人口 (见正文p12)
	gen gdp_cost = dgrain_0316*coef_base*mgdppc_0316*mpop_0316
	//转换单位：万元-->亿元
	replace gdp_cost = gdp_cost/10000 
	* ii. 样本县增产成本加总
	bysort year: egen ttl_gdp_cost = total(gdp_cost)
		
	
	* iii. 人均成本=产量变化*估计系数*人均GDP
	gen gdppc_cost = dgrain_0316*coef_base*mgdppc_0316
	
	* iv. 按人口加权平均得到全国人均成本
	* 计算样本总人口
	egen na_mpop_0316 = total(mpop_0316)
	
	egen na_gdppc_cost = total(gdppc_cost*mpop_0316/na_mpop_0316)

	
	* B.产销区按照不同系数估计成本
	* i.各县增产成本(GDP)=产量变化*估计系数*人均GDP*人口 (见正文p12)
	gen gdp_cost_cq = dgrain_0316*coef_cq*mgdppc_0316*mpop_0316
	//转换单位：亿元
	replace gdp_cost_cq = gdp_cost_cq/10000
	* ii. 样本县增产成本加总
	bysort year: egen ttl_gdp_cost_cq = total(gdp_cost_cq)
	
	
	* iii. 人均成本=产量变化*估计系数*人均GDP
	gen gdppc_cost_cq = dgrain_0316*coef_cq*mgdppc_0316
	* iv. 按人口加权平均得到全国人均成本
	* 计算样本总人口
	bysort year: egen na_gdppc_cost_cq = total(gdppc_cost_cq*mpop_0316/na_mpop_0316)
	
	
	
	** (3).输出县样本结果
	
	*A.按基准估计：总量
	preserve
	* 避免同一个变量在collapse生成不同统计量时冲突
	gen gdp_cost_formedian = gdp_cost
	
	collapse (mean) mean_cost = gdp_cost (median) median_cost = gdp_cost_formedian (mean) total_cost = ttl_gdp_cost
	gen row = "A1"
	save $logfiles/table_4/coun_cost_a1.dta,replace
	restore

	
	*A.按基准估计：人均
	preserve
	* 人口取整，便于作为权重
    gen mpop_0316_int = round(mpop_0316)
	collapse (mean) mean_cost = na_gdppc_cost (median) median_cost = gdppc_cost [fweight = mpop_0316_int]
	gen row = "A2"
	save $logfiles/table_4/coun_cost_a2.dta,replace
	
	restore
	
	
	
	*B.产销区不同系数估计：总量
	preserve
	* 避免同一个变量在collapse生成不同统计量时冲突
	gen gdp_cost_formedian_cq = gdp_cost_cq
	collapse (mean) mean_cost = gdp_cost_cq (median) median_cost = gdp_cost_formedian_cq (mean) total_cost = ttl_gdp_cost_cq
	gen row = "B1"
	save $logfiles/table_4/coun_cost_b1.dta,replace
	restore
	
	
	*B.产销区不同系数估计：人均
	preserve
	* 人口取整，便于作为权重
    gen mpop_0316_int = round(mpop_0316)
	collapse (mean) mean_cost = na_gdppc_cost_cq (median) median_cost = gdppc_cost_cq [fweight = mpop_0316_int]
	gen row = "B2"
	save $logfiles/table_4/coun_cost_b2.dta,replace
	
	restore
	
	
	
	*结果合并：
	use $logfiles/table_4/coun_cost_a1.dta,clear
	append using $logfiles/table_4/coun_cost_a2.dta
	append using $logfiles/table_4/coun_cost_b1.dta
	append using $logfiles/table_4/coun_cost_b2.dta
	save $logfiles/table_4/coun_cost.dta,replace

	
	
	
	
							*2. 用省样本估算
	
	use $working_data/province_data.dta,clear


	** (1).获取基准估计系数，以及分产区的估计系数
	
	* 获取县样本的估计系数
	gen coef_base = $ebase
	gen coef_cq = $ezcq if is_zcq==1
	replace coef_cq = $ezxq if is_zxq==1
	replace coef_cq = $ephq if is_phq==1
	

	** (2).计算增产成本
	* 样本期间平均GDP、总人口和县均人口和粮食产量变化
	gen is_yr2003 = 1 if year==2003
	gen is_yr2016 = 1 if year==2016
	
	foreach v of varlist gdppc mpop_per_county grain pop{
		
		bysort province_name: egen `v'_03 = mean(`v'*is_yr2003)
		bysort province_name: egen `v'_16 = mean(`v'*is_yr2016)
		
	}
	
	* 增产(单位转换：万吨->十万吨)
	gen dgrain_0316 = (grain_16 - grain_03)/10
	bysort province_name: gen mgdppc_0316 = (gdppc_03 + gdppc_16)/2
	bysort province_name: gen mpop_per_county_0316 = (mpop_per_county_03 + mpop_per_county_16)/2
	bysort province_name: gen mpop_0316 = (pop_03 + pop_16)/2
	
	drop if dgrain_0316==.|mgdppc_0316==.|mpop_0316==.
	keep if year==2016
	
	
	* A.按照基准系数估计成本
	* i.各省增产成本(GDP)=产量变化*估计系数*人均GDP*县平均人口 (见正文p12)
	gen gdp_cost = (-1)*dgrain_0316*0.059*mgdppc_0316*mpop_per_county_0316
	//单位转换为亿元
	replace gdp_cost = gdp_cost/10000
	* ii.各省成本加总
	egen ttl_gdp_cost = total(gdp_cost)
		
	
	* iii.各省人均成本=总成本/总人口数，再将单位转换为元/人
	gen gdppc_cost_0316 = (gdp_cost/mpop_0316)*10000
	* iv. 按人口加权平均得到全国人均成本
	* 计算样本总人口
	egen na_mpop_0316 = total(mpop_0316)
	egen na_gdppc_cost = total(gdppc_cost_0316*mpop_0316/na_mpop_0316)
	
		
	
	* B.按照产销区不同系数估计成本
	* i.各省增产成本(GDP)=产量变化*估计系数*人均GDP*县平均人口 (见正文p12)
	gen gdp_cost_cq = dgrain_0316*coef_cq*mgdppc_0316*mpop_per_county_0316
	//单位转换为亿元
	replace gdp_cost_cq = gdp_cost_cq/10000
	* ii.各省成本加总
	bysort year: egen ttl_gdp_cost_cq = total(gdp_cost_cq)
	
	
	* iii.各省人均成本=总成本/总人口数，再将单位转换为元/人
	gen gdppc_cost_cq = (gdp_cost_cq/mpop_0316)*10000
	* iv. 按人口加权平均得到全国人均成本
	* 计算样本总人口
	egen na_gdppc_cost_cq = total(gdppc_cost_cq*mpop_0316/na_mpop_0316)
	
	* 将总人口数保存到全局暂元中以便后续计算人均转移支付
	gl ttlpop = na_mpop_0316[1]
	
	
	** (3).输出省样本结果
	*A.按基准估计：总量
	preserve
	collapse (mean) cost_by_prov = ttl_gdp_cost
	gen row = "A1"
	save $logfiles/table_4/prov_cost_a1.dta,replace
	restore
	
	
	*A.按基准估计：人均
	preserve
	collapse (mean) cost_by_prov = na_gdppc_cost
	gen row = "A2"
	save $logfiles/table_4/prov_cost_a2.dta,replace
	restore
	
		
	*B.产销区不同系数估计：总量
	preserve
	collapse (mean) cost_by_prov = ttl_gdp_cost_cq
	gen row = "B1"
	save $logfiles/table_4/prov_cost_b1.dta,replace
	restore
	
	
	*B.按基准估计：人均
	preserve
	collapse (mean) cost_by_prov = na_gdppc_cost_cq
	gen row = "B2"
	save $logfiles/table_4/prov_cost_b2.dta,replace
	restore
		
		
	*结果合并：
	use $logfiles/table_4/prov_cost_a1.dta,clear
	append using $logfiles/table_4/prov_cost_a2.dta
	append using $logfiles/table_4/prov_cost_b1.dta
	append using $logfiles/table_4/prov_cost_b2.dta
	save $logfiles/table_4/prov_cost.dta,replace
		
	
	
							*3. 汇总结果
	
	use $logfiles/table_4/coun_cost.dta ,clear
	merge 1:1 row using $logfiles/table_4/prov_cost.dta
	drop _merge
	
	
	* (1).已实现的转移支付
	** 总量= 13215.05亿元 (见附表A4)
	gen transfer = 13215.05
	** 总人口数(之前已经保存到ttlpop暂元中)
	gen pop = $ttlpop
	** 人均：除以，并调整单位为元/人
	replace transfer = (transfer/pop)*10000 if row=="A2"|row=="B2"	
	
	
	* (2).真实总成本
	** 总量 = 还未覆盖的成本(省样本) + 转移支付
	gen realcost = cost_by_prov - transfer
	** 转移支付占总成本百分比(因为总成本用负值表示，转移支付用正值，这里求百分比时取绝对值)
	gen transfer_shr = abs((transfer/realcost)*100)
	
	drop pop
	
	
	** 输出表格
	gen rownames = "A.按基准估计：总量" if row=="A1"
	replace rownames = "A.按基准估计：人均" if row=="A2"
	replace rownames = "B.分产销区异质性：总量" if row=="B1"
	replace rownames = "B.分产销区异质性：人均" if row=="B2"
	
	//保留两位小数
	foreach v of varlist *cost cost_* transfer*{
		
		gen `v'_2dig = string(round(`v',0.01),"%20.2f") 
	}
	
	keep row rownames *_2dig 
	renvarlab *_2dig,subst(_2dig )
	
	label variable mean_cost "平均值"
	label variable median_cost "中位数"
	label variable total_cost "总和"
	label variable cost_by_prov "省样本估算"
	label variable transfer "已实现的转移支付"
	label variable transfer_shr "转移支付占比(%)"
	label variable realcost "真实总成本"
	
	order row,first
	order row rownames mean_cost median_cost total_cost cost_by_prov transfer transfer_shr realcost
	
	
	export excel using "$logfiles/table_4/table_4.xls" , firstrow(varlabels) replace
	
	
	
	
********************************************************************************
**************************** 表5.利益补偿机制的比较 ****************************
********************************************************************************


	use $working_data/county_panel.dta ,clear
	capture mkdir $logfiles/table_5
	
	
	
	** (1).计算样本县增产成本，以便后续计算补偿标准
	
	* 生成变量的基期(2003)值
	gen is_yr2003 = 1 if year==2003
	foreach v of varlist Grain lnGDP_pc{
		
		bysort county_id (year): egen `v'_03 = mean(`v'*is_yr2003)
		
	}
	
	* 生成相对基期的差分
	foreach v of varlist Grain lnGDP_pc{
		
		gen d_`v' = `v' - `v'_03
		
	}

	* 基准 (保留小数点后三位，不会影响结论)
	ivreghdfe d_lnGDP_pc lnGDP_pc_03 (d_Grain = dsim_Grain) if year==2016, absorb(city_id) vce(cluster city_id)
	gen coef_base = round(_b[d_Grain],0.001)
	gl ebase = round(_b[d_Grain],0.001)
	
	* 分产区
	gen coef_cq = .
	ivreghdfe d_lnGDP_pc lnGDP_pc_03 (d_Grain = dsim_Grain) if year==2016&is_zcq==1, absorb(city_id) vce(cluster city_id)
	replace coef_cq = round(_b[d_Grain],0.001) if is_zcq==1
	gl ezcq = round(_b[d_Grain],0.001)
	
	ivreghdfe d_lnGDP_pc lnGDP_pc_03 (d_Grain = dsim_Grain) if year==2016&is_zxq==1, absorb(city_id) vce(cluster city_id)
	replace coef_cq = round(_b[d_Grain],0.001) if is_zxq==1
	gl ezxq = round(_b[d_Grain],0.001)

	ivreghdfe d_lnGDP_pc lnGDP_pc_03 (d_Grain = dsim_Grain) if year==2016&is_phq==1, absorb(city_id) vce(cluster city_id)
	replace coef_cq = round(_b[d_Grain],0.001) if is_phq==1
	gl ephq = round(_b[d_Grain],0.001)
	
		
	* 计算样本期间各县增产量，平均人口、人均GDP以及GDP
	gen is_yr2016 = 1 if year==2016
	foreach v of varlist GDP_pc pop GDP farmland_area{
		
		bysort county_id (year): egen `v'_03 = mean(`v'*is_yr2003)
		bysort county_id (year): egen `v'_16 = mean(`v'*is_yr2016)
		
	}
	
	bysort county_id (year): egen Grain_16 = mean(Grain*is_yr2016)
	gen dgrain_0316 = Grain_16 - Grain_03
	gen mgdppc_0316 = (GDP_pc_03 + GDP_pc_16)/2
	gen mpop_0316 = (pop_03 + pop_16)/2
	gen mgdp_0316 = mgdppc_0316*mpop_0316
	
	drop if dgrain_0316==.|mgdppc_0316==.|mpop_0316==.
	keep if year==2016
	
	* 按照基准系数估计成本
	* i.各县增产成本(GDP)=产量变化*估计系数*人均GDP*人口 (见正文p12)
	gen gdp_cost = dgrain_0316*coef_base*mgdppc_0316*mpop_0316
	//转换单位：万元-->亿元
	replace gdp_cost = gdp_cost/10000 
	* ii. 样本县增产成本加总
	bysort year: egen ttl_gdp_cost = total(gdp_cost)

	* iii. 人均成本=产量变化*估计系数*人均GDP
	gen gdppc_cost = dgrain_0316*coef_base*mgdppc_0316
	
	* iv. 按人口加权平均得到全国人均成本
	egen na_mpop_0316 = total(mpop_0316)
	egen na_gdppc_cost = total(gdppc_cost*mpop_0316/na_mpop_0316)
	

	
	** (2). 2016年事实的区域发展差距
	
	* 基尼系数
	* 16年人口(万人->人并取整)作为权重
	gen pop_16_int = round(pop*10000)
	ineqdeco GDP_pc_16 [fweight = pop_16_int]
	gen gini_real = round(r(gini),0.001)
	
	* 人均GDP，高产组:低产组
	egen high_dgrain = xtile(dgrain_0316),nquantiles(2)
	//高组=1，低组=0
	replace high_dgrain = high_dgrain - 1
	
	//组内均值(人口加权)
	bysort high_dgrain: egen ttlpop_by_dgrain = total(pop_16)
	bysort high_dgrain: egen mgdppc_by_dgrain = total(GDP_pc_16*pop_16/ttlpop_by_dgrain)
	
	sum mgdppc_by_dgrain if high_dgrain==1,meanonly
	gen gdppc_high_real = r(mean)
	
	sum mgdppc_by_dgrain if high_dgrain==0,meanonly
	gen gdppc_low_real = r(mean)
	
	gen ratio_gdppc_real = round(gdppc_high_real/gdppc_low_real,0.001)
	
	
	
	** (3).按不同标准补偿并计算补偿后的发展差距
	

	* A.按人口补偿
	* 16年样本县总人口(万人)
	egen napop_16 = total(pop_16)
    * 计算补偿标准并换算成元/人
	gen stad_bypop = abs((ttl_gdp_cost/napop_16)*10000)
	* 补偿后各县反事实GDP(万元)
	gen cfgdp_bypop = GDP_16+stad_bypop*pop_16
	* 补偿后各县反事实人均GDP(元/人)
	gen cfgdppc_bypop = cfgdp_bypop/pop_16
	
	* 计算反事实后的区域发展差距
	* 基尼系数
	ineqdeco cfgdppc_bypop [fweight = pop_16_int]
	gen gini_cfpop = round(r(gini),0.001)
	
	* 人均GDP，高产组:低产组
	//组内均值(人口加权)
	bysort high_dgrain: egen mgdppc_by_dgrain_cfpop = total(cfgdppc_bypop*pop_16/ttlpop_by_dgrain)
	
	sum mgdppc_by_dgrain_cfpop if high_dgrain==1,meanonly
	gen gdppc_high_cfpop = r(mean)
	
	sum mgdppc_by_dgrain_cfpop if high_dgrain==0,meanonly
	gen gdppc_low_cfpop = r(mean)
	
	gen ratio_gdppc_cfpop = round(gdppc_high_cfpop/gdppc_low_cfpop ,0.001)
	
	
	
	
		
	* B.按照粮食增产补偿
	* 样本县的总增产量(十万吨)
	egen na_dgrain_0316 = total(dgrain_0316)
	* 计算补偿标准(亿元/十万吨)
	gen stad_bydgrain = abs(ttl_gdp_cost/na_dgrain_0316)
	* 换算成元/吨
	gen stad_bydgrain_ypt = stad_bydgrain*1000
	
	* 补偿后反事实GDP(万元)
	gen cfgdp_bydgrain = GDP_16+stad_bydgrain*dgrain_0316*10000
	* 补偿后反事实人均GDP(元/人)
	gen cfgdppc_bydgrain = cfgdp_bydgrain/pop_16
	
	
	* 计算反事实后的区域发展差距
	* 基尼系数
	ineqdeco cfgdppc_bydgrain [fweight = pop_16_int]
	gen gini_cfgrn = round(r(gini),0.001)
	* 人均GDP，高产组:低产组
	//组内均值(人口加权)
	bysort high_dgrain: egen mgdppc_by_dgrain_cfgrn = total(cfgdppc_bydgrain*pop_16/ttlpop_by_dgrain)
	
	sum mgdppc_by_dgrain_cfgrn if high_dgrain==1,meanonly
	gen gdppc_high_cfgrn = r(mean)
	
	sum mgdppc_by_dgrain_cfgrn if high_dgrain==0,meanonly
	gen gdppc_low_cfgrn = r(mean)
	
	gen ratio_gdppc_cfgrn = round(gdppc_high_cfgrn/gdppc_low_cfgrn ,0.001)
	
	
	
	
	* C.按耕地面积
	* 样本县总耕地面积
	egen nacrld_16 = total(farmland_area_16)
	* 按耕地面积补偿标准
	gen stad_byfmld = abs(ttl_gdp_cost/nacrld_16)
	* 换算成元/亩
	gen stad_byfmld_ypm = stad_byfmld*10000/1500

	* 补偿后的反事实GDP(万元)
	gen cfgdp_byfmld = GDP_16+stad_byfmld*farmland_area_16*10000
	* 补偿后反事实人均GDP(元/人)
	gen cfgdppc_byfmld = cfgdp_byfmld/pop_16
	
	
	
	* 计算反事实后的区域发展差距
	* 基尼系数
	ineqdeco cfgdppc_byfmld [fweight = pop_16_int]
	gen gini_cfland = round(r(gini),0.001)
	
	* 人均GDP，高产组:低产组
	//组内均值(人口加权)
	bysort high_dgrain: egen mgdppc_by_dgrain_cfland = total(cfgdppc_byfmld*pop_16/ttlpop_by_dgrain)
	
	sum mgdppc_by_dgrain_cfland if high_dgrain==1,meanonly
	gen gdppc_high_cfland = r(mean)
	
	sum mgdppc_by_dgrain_cfland if high_dgrain==0,meanonly
	gen gdppc_low_cfland = r(mean)
	
	gen ratio_gdppc_cfland = round(gdppc_high_cfland/gdppc_low_cfland ,0.001)

	
	
	
	** (4).输出结果
	
	preserve
	* 展示补偿标准的绝对值(由负)
	collapse (mean) standard = stad_bypop gini_real gini_cf = gini_cfpop ratio_gdppc_real ratio_gdppc_cf = ratio_gdppc_cfpop
	gen row = 1
	gen rownames = "按人口平分"
	save $logfiles/table_5/compen_1.dta,replace
	restore
	

	preserve
	collapse (mean) standard = stad_bydgrain_ypt gini_real gini_cf = gini_cfgrn ratio_gdppc_real ratio_gdppc_cf = ratio_gdppc_cfgrn
	gen row = 2
	gen rownames = "按粮食增产"
	save $logfiles/table_5/compen_2.dta,replace
	restore
	
	
	preserve
	collapse (mean) standard = stad_byfmld_ypm gini_real gini_cf = gini_cfland ratio_gdppc_real ratio_gdppc_cf = ratio_gdppc_cfland
	gen row = 3
	gen rownames = "按耕地存量"
	save $logfiles/table_5/compen_3.dta,replace
	restore
	
	
	* 合并输出
	use $logfiles/table_5/compen_1.dta ,clear
	append using $logfiles/table_5/compen_2.dta
	append using $logfiles/table_5/compen_3.dta
	replace standard = round(standard,0.01)
	
	//保留小数位数
	gen standard_str = string(standard,"%20.2f")
	foreach v of varlist gini_* ratio_*{
		
		gen `v'_str = string(`v',"%20.3f")
	}
	
	keep row rownames *_str
	renvarlab *_str,subst(_str )
	
	
	label variable standard "补偿标准"
	label variable gini_real "基尼系数：2016事实"
	label variable gini_cf "基尼系数：一次性补偿后"
	label variable ratio_gdppc_real "人均GDP之比：2016事实"
	label variable ratio_gdppc_cf "人均GDP之比：一次性补偿后"
	order row rownames standard gini_real gini_cf ratio_gdppc_real ratio_gdppc_cf
	
	export excel using $logfiles/table_5/table_5.xls, firstrow(varlabels) replace
	
	
	
	
	
	
	
	
	log close
