********************************************************************************
/*						本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_appendix, 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: 省总人口
	
	
*/

	
	


********************************************************************************
********************** 图I1.政府对粮食生产的关注度 ******************************
********************************************************************************

	use $working_data/word_freq.dta ,clear

	twoway (connected word_freq_grain year if year>=2000&year<=2023, fcolor(none) lcolor(black) lpattern(solid) lwidth(thin) msymbol(o) msize(medium) mcolor(black) ) ||, xtitle("年份",size(small) alignment(top) margin(small)) xlabel(. 2000 2003 2016 2023,labsize(small)) xline(2003,lcolor(black) lstyle(refline) lwidth(vthin) lpattern(dash)) xline(2016,lcolor(black) lstyle(refline) lwidth(vthin) lpattern(dash)) ytitle("粮食生产相关词频(‰)", size(small) alignment(top) margin(medsmall)) ylabel(,labsize(small)) scheme(s1mono) legend(label(1 "粮食词频") size(small))  subtitle("省级政府工作报告中粮食生产词频",size(medium))
	
	graph export "$logfiles/figure_I1.png",replace width(1600) height(1200)



********************************************************************************
**************************** 图I2 .粮食生产在县层面的分化  ***********************
********************************************************************************


	use $working_data\county_panel.dta ,clear

	
			        	***--- 变量处理 ---***	
	
	bysort year: egen mean_grain = mean(Grain)
	gen Grain_dmean = Grain - mean_grain
	
						***--- 作图 ---***	
	
	twoway kdensity Grain_dmean if year==2003,lpattern(longdash) lwidth(medium) lcolor(gs5) || kdensity Grain_dmean if year==2016,lpattern(solid) lwidth(medium) lcolor(black)||, xtitle("均值中心化粮食产量(十万吨)",size(medlarge)) ytitle("核密度",size(medlarge) margin(medsmall)) xline(0,lcolor(black) lpattern(dash)) legend(label(1 "2003年") label(2 "2016年") ring(0) bplacement(neast) rows(1) size(small)) name(crpdist,replace) title("a.全国", size(Large)) scheme(s1mono) 
	
	
	twoway kdensity Grain_dmean if year==2003&is_zcq==1,lpattern(longdash) lwidth(medium) lcolor(gs5) || kdensity Grain_dmean if year==2016&is_zcq==1,lpattern(solid) lwidth(medium) lcolor(black)||, xtitle("均值中心化粮食产量(十万吨)",size(medlarge)) ytitle("核密度",size(medlarge) margin(medsmall)) xline(0,lcolor(black) lpattern(dash)) legend(label(1 "2003年") label(2 "2016年") ring(0) bplacement(neast) rows(1) size(small)) name(crpdist_zcq,replace) title("b.主产区", size(Large)) scheme(s1mono) ylabel(0(0.1)0.4)
	
	graph combine crpdist crpdist_zcq, cols(2) scheme(s1mono) ysize(4cm) xsize(10cm)
	graph export "$logfiles/figure_I2.png",replace width(2400) height(1100)





********************************************************************************
******************* 表I1 .2003-2016年分区域的粮食增产情况 ************************
********************************************************************************
	
	
	use $working_data\province_data.dta ,clear
	capture mkdir "$logfiles/table_I1"
	
	* 计算各省2003-2016的粮食增产
	gen is_yr2003 = 1 if year==2003
	gen is_yr2016 = 1 if year==2016
	
	bysort province_name: egen grain_03 = mean(grain*is_yr2003)
	bysort province_name: egen grain_16 = mean(grain*is_yr2016)
	
	gen d_grain = grain_16 - grain_03
	
	* 分别按照省份、产销区和地理区域统计粮食增产，输出三个文件，手动整理后形成最终展示的表格
	
	* a.分省
	preserve
	collapse d_grain if year==2016, by(province_name)
	eststo clear
	gsort - d_grain 
	gen d_grain_fmt = string(round(d_grain,0.001),"%20.3f")  //保留三位小数
	label variable d_grain_fmt "粮食增产(万吨)"
	label variable province_name "省份"
	export excel province_name d_grain_fmt using "$logfiles/table_I1/增产_按省.xlsx", firstrow(varlabels) replace 
	restore
	
	
	
	* b.分主产区、主销区和平衡区
	preserve
	collapse (sum) d_grain if year==2016, by(prod_region)
	eststo clear
	gsort - d_grain 
	gen d_grain_fmt = string(round(d_grain,0.001),"%20.3f")  //保留三位小数
	label variable d_grain_fmt "粮食增产(万吨)"
	label variable prod_region "产区"
	export excel prod_region d_grain_fmt using "$logfiles/table_I1/增产_按产区.xlsx", firstrow(varlabels) replace 
	restore
			
	
	
	* c.分地理区域
	preserve
	collapse (sum) d_grain if year==2016, by(geo_region)
	eststo clear
	gsort - d_grain 
	gen d_grain_fmt = string(round(d_grain,0.001),"%20.3f")  //保留三位小数
	label variable d_grain_fmt "粮食增产(万吨)"
	label variable geo_region "地理区域"
	export excel geo_region d_grain_fmt using "$logfiles/table_I1/增产_按地理区域.xlsx", firstrow(varlabels) replace 
	restore
	
	
	
	
*******************************************************************************
**************** 表A1 .按自然禀赋的增产和实际增产的相关性 ****************
*******************************************************************************

	use $working_data\county_panel.dta ,clear
	
						***--- 变量处理 ---***	
	
	* 生成变量的基期(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) cluster(city_id)
	gen baseline_sample = e(sample)
	
	
	* 用检验实际增产与潜在增产(用不同IV代表)的关系
	
	* (1). 基准IV
	reghdfe d_Grain dsim_Grain lnGDP_pc_03 if baseline_sample==1, absorb(city_id) vce(cluster city_id)
	
	outreg2 using "$logfiles/table_A1", replace excel nocons depvar dec(3) keep(dsim_Grain)
	
	* (2). 四种主粮
	
	reghdfe d_Grain dsim_mGrain lnGDP_pc_03 if baseline_sample==1, absorb(city_id) vce(cluster city_id)
	
	outreg2 using "$logfiles/table_A1", append excel nocons depvar dec(3) keep(dsim_mGrain)
	
	* (3). 三种主粮
	
	reghdfe d_Grain dsim_mmGrain lnGDP_pc_03 if baseline_sample==1, absorb(city_id) vce(cluster city_id)
	
	outreg2 using "$logfiles/table_A1", append excel nocons depvar dec(3) keep(dsim_mmGrain)
	
	* (4). 所有作物
	reghdfe d_Grain dsim_allcrop lnGDP_pc_03 if baseline_sample==1, absorb(city_id) vce(cluster city_id)
	
	outreg2 using "$logfiles/table_A1", append excel nocons depvar dec(3) keep(dsim_allcrop)
	
	
	
	
		
	
	
********************************************************************************
*********************** 表A2 .主要变量的描述性统计 *****************************
********************************************************************************

	use $working_data\county_panel.dta ,clear
	


	* 转换变量单位,便于展示
	* 比例变为百分比
	foreach v of varlist Share_PriGDP Share_SecGDP share_ov1k_new_firms{
		
		gen `v'_pct = `v'*100
		
	}
	
	* 产量由十万吨->万吨
	replace dsim_Grain = dsim_Grain*10
	replace Grain = Grain*10
	
	
	global vlist lnGDP_pc lnGDP_agr_pc lnGDP_ind_pc lnGDP_serv_pc lncstrct_land lnavgemp_old_firms lnpub_expd lnpub_expd_agr Grain dsim_Grain Share_PriGDP_pct Share_SecGDP_pct lnRainfall lnTempr
	
	* 样本期
	keep if year==2003|year==2016
	* 保持和回归样本一致
	drop if Grain==.|dsim_Grain==.|lnGDP_pc==.
	
	eststo clear
	estpost tabstat $vlist, statistics(mean sd n min max) columns(statistics)
	matrix sm_mean  = e(mean)
	matrix sm_sd = e(sd)
	matrix sm_n = e(count)
	matrix sm_min = e(min)
	matrix sm_max = e(max)
	
	matrix colnames sm_mean = $vlist
	matrix colnames sm_sd = $vlist
	matrix colnames sm_n = $vlist
	matrix colnames sm_min = $vlist
	matrix colnames sm_max = $vlist
	
	esttab using "$logfiles/table_A2.rtf" ,cells("mean(fmt(%20.3f) label(Mean)) sd(fmt(%20.3f) label(Std.dev.)) count(fmt(%g) label(Obs.)) min(fmt(%20.3f) label(Min)) max(fmt(%20.3f) label(Max))") noobs nostar nonumbers nomtitles replace 




	
********************************************************************************
*************************** 表A3.粮食增产与劳动力配置 ***************************
********************************************************************************


	use $working_data/census_data.dta ,clear

	
						***--- 变量处理 ---***	
	
	* 受限于普查数据可得性，以2000年作为基期
	gen is_yr2000 = 1 if year==2000
	foreach v of varlist Grain lnGDP_pc lnemply emply_shr*{
		
		bysort county_id (year): egen `v'_00 = mean(`v'*is_yr2000)
		
	}
	
	
	* 生成相对基期的差分
	foreach v of varlist Grain lnemply emply_shr_pri emply_shr_sec emply_shr_ter{
		
		gen d_`v' = `v' - `v'_00
		
	}
	
	

						***--- 逐年差分估计 ---***
		
	
	foreach y of varlist d_lnemply d_emply_shr_pri d_emply_shr_sec d_emply_shr_ter{
		
		ivreghdfe `y' lnGDP_pc_00 (d_Grain = dsim_Grain) if year==2010, absorb(city_id) vce(cluster city_id)
		outreg2 using "$logfiles/table_A3", ctitle(`y') addtext(Controls, Y, City FE, Y) append excel nocons nor2 depvar addstat(K-P Wald F stat, e(rkf)) keep(d_Grain) adec(3) dec(3)
		
		
	}
	

********************************************************************************
********************* 表A4. 已实现的转移支付数额*********************
********************************************************************************


	use $working_data/transfers.dta ,clear
	
	* 平减产粮大县奖励(transfer_daxian)和农业补贴(agr_subsidy)
	foreach v of varlist transfer_daxian agr_subsidy{
		
		gen `v'_df = `v'/cpi
		
	}
	
	* 平减后总额
	foreach v of varlist transfer_daxian agr_subsidy{
		
		egen tl_`v' = total(`v'_df)
		
	}
	
	
	
	* 每年度平减后的总转移支付
	egen tl_transfer = rowtotal(transfer_daxian_df agr_subsidy_df)
	
	
	* 所有年度总转移支付
	egen total_transfers = total(tl_transfer)
	

	* 输出表格
	* 在最后一行补充总和
	gl tl_transfer_daxian = tl_transfer_daxian[1]
	gl tl_agr_subsidy = tl_agr_subsidy[1]
	gl total_transfers = total_transfers[1]
	
	set obs 14
	replace transfer_daxian_df = $tl_transfer_daxian in 14
	replace agr_subsidy_df = $tl_agr_subsidy in 14
	replace tl_transfer = $total_transfers in 14
	
	drop cpi tl_transfer_daxian tl_agr_subsidy total_transfers
	
	
	tostring year,replace
	replace year = "总和" in 14
	
	
	* 保留小数点后两位
	foreach v of varlist transfer_daxian transfer_daxian_df agr_subsidy agr_subsidy_df tl_transfer{
		
		gen `v'_2d = string(round(`v',0.01),"%20.2f")
		
	}
	
	keep year *_2d
	renvarlab *_2d,subst(_2d )
	order year transfer_daxian transfer_daxian_df agr_subsidy agr_subsidy_df tl_transfer, first
	
	
	label variable year "年份"
	label variable transfer_daxian "产粮大县奖励"
	label variable transfer_daxian_df "产粮大县奖励   （2000年价）"
	label variable agr_subsidy "农业“四项”补贴"
	label variable agr_subsidy_df "农业“四项”补贴（2000年价）"
	label variable tl_transfer "加总（2000年价）"

	
	
	export excel using $logfiles/table_A4.xlsx,firstrow(varlabels) replace
	
	
	
	
	
********************************************************************************
**************** 表 A5.利益补偿机制设计和补偿效果：产销区不同标准 ****************
********************************************************************************

	use $working_data/county_panel.dta ,clear
	capture mkdir $logfiles/table_A5
	
	
	
	** (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)
	
	
	* 计算样本期间各县增产量，平均人口、人均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

	
	
	** (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.按人口补偿
	* 各产销区的总人口(万人)
	 bysort cqtype: egen cq_pop_16 = total(pop_16)
	* 各产销区的总成本(万元)
	 bysort cqtype: egen cq_gdp_cost = total((-1)*gdp_cost)
	 * 各产销区的补偿标准(元/人)
	 bysort cqtype: gen stad_bypop_cq = cq_gdp_cost/cq_pop_16
	 
	 * 补偿后反事实GDP	 
	 gen cfgdp_bypop_cq = GDP_16+stad_bypop_cq*pop_16
	 * 补偿后反事实人均GDP
	 gen cfgdppc_bypop_cq = cfgdp_bypop_cq/pop_16
	
	
	* 计算反事实后的区域发展差距
	* 基尼系数
	 ineqdeco cfgdppc_bypop_cq [fweight = pop_16_int]
	 gen gini_cfpop_cq = round(r(gini),0.001)
	 
	 
	* 人均GDP，高产组:低产组
	//组内均值(人口加权)
	bysort high_dgrain: egen mgdppc_by_dgrain_cfpop_cq = total(cfgdppc_bypop_cq*pop_16/ttlpop_by_dgrain)
	
	sum mgdppc_by_dgrain_cfpop_cq if high_dgrain==1,meanonly
	gen gdppc_high_cfpop_cq = r(mean)
	
	sum mgdppc_by_dgrain_cfpop_cq if high_dgrain==0,meanonly
	gen gdppc_low_cfpop_cq = r(mean)
	
	gen ratio_gdppc_cfpop_cq = round(gdppc_high_cfpop_cq/gdppc_low_cfpop_cq ,0.001)
	
	
					
	* B.按粮食增产补偿
	* 各产区总增产
	bysort cqtype: egen cq_dgrain_0316 = total(dgrain_0316)
	* 各产区补偿标准
	bysort cqtype: gen stad_bygrn_cq = cq_gdp_cost/cq_dgrain_0316
	* 换算单位元/吨
	gen stad_bygrn_cq_ypt = stad_bygrn_cq/10
	* 补偿后的反事实GDP
	gen cfgdp_bygrn_cq = GDP_16+stad_bygrn_cq*dgrain_0316
	* 补偿后反事实人均GDP
	gen cfgdppc_bygrn_cq = cfgdp_bygrn_cq/pop_16
	
	
	* 计算反事实后的区域发展差距
	* 基尼系数
	 ineqdeco cfgdppc_bygrn_cq [fweight = pop_16_int]
	 gen gini_cfgrn_cq = round(r(gini),0.001)
	
	* 人均GDP，高产组:低产组
	//组内均值(人口加权)
	bysort high_dgrain: egen mgdppc_by_dgrain_cfgrn_cq = total(cfgdppc_bygrn_cq*pop_16/ttlpop_by_dgrain)
	
	sum mgdppc_by_dgrain_cfgrn_cq if high_dgrain==1,meanonly
	gen gdppc_high_cfgrn_cq = r(mean)
	
	sum mgdppc_by_dgrain_cfgrn_cq if high_dgrain==0,meanonly
	gen gdppc_low_cfgrn_cq = r(mean)
	
	gen ratio_gdppc_cfgrn_cq = round(gdppc_high_cfgrn_cq/gdppc_low_cfgrn_cq ,0.001)
		
		

	* C.按耕地面积补偿
	* 各产销区总耕地面积
	bysort cqtype: egen cq_fmld_16 = total(farmland_area_16)
	* 各产销区补偿标准
	bysort year cqtype: gen stad_byland_cq = cq_gdp_cost/cq_fmld_16
	* 单位换算成元/亩
	gen stad_byland_cq_ypm = stad_byland_cq/1500
	
	* 补偿后反事实GDP
	gen cfgdp_byland_cq = GDP_16+stad_byland_cq*farmland_area_16
	* 补偿后反事实人均GDP	
	gen cfgdppc_byland_cq = cfgdp_byland_cq/pop_16
	
	
	* 计算反事实后的区域发展差距
	* 基尼系数
	 ineqdeco cfgdppc_byland_cq [fweight = pop_16_int]
	 gen gini_cfland_cq = round(r(gini),0.001)
	
	* 人均GDP，高产组:低产组
	//组内均值(人口加权)
	bysort high_dgrain: egen mgdppc_by_dgrain_cfland_cq = total(cfgdppc_byland_cq*pop_16/ttlpop_by_dgrain)
	
	sum mgdppc_by_dgrain_cfland_cq if high_dgrain==1,meanonly
	gen gdppc_high_cfland_cq = r(mean)
	
	sum mgdppc_by_dgrain_cfland_cq if high_dgrain==0,meanonly
	gen gdppc_low_cfland_cq = r(mean)
	
	gen ratio_gdppc_cfland_cq = round(gdppc_high_cfland_cq/gdppc_low_cfland_cq ,0.001)

	
	
	
	** (4).输出结果
	
	preserve
	collapse (mean) standard = stad_bypop_cq gini_real gini_cf = gini_cfpop_cq ratio_gdppc_real ratio_gdppc_cf = ratio_gdppc_cfpop_cq,by(cqtype)
	gen row = 1
	gen rownames = "按人口平分"
	save $logfiles/table_A5/compen_cq_1.dta,replace
	restore
	

	preserve
	collapse (mean) standard = stad_bygrn_cq_ypt gini_real gini_cf = gini_cfgrn_cq ratio_gdppc_real ratio_gdppc_cf = ratio_gdppc_cfgrn_cq,by(cqtype)
	gen row = 2
	gen rownames = "按粮食增产"
	save $logfiles/table_A5/compen_cq_2.dta,replace
	restore
	
	
	preserve
	collapse (mean) standard = stad_byland_cq_ypm gini_real gini_cf = gini_cfland_cq ratio_gdppc_real ratio_gdppc_cf = ratio_gdppc_cfland_cq,by(cqtype)
	gen row = 3
	gen rownames = "按耕地存量"
	save $logfiles/table_A5/compen_cq_3.dta,replace
	restore
	
	
	* 合并输出
	use $logfiles/table_A5/compen_cq_1.dta ,clear
	append using $logfiles/table_A5/compen_cq_2.dta
	append using $logfiles/table_A5/compen_cq_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 rownames *_str cqtype
	renvarlab *_str,subst(_str )
	
	gen cq = "主产区" if cqtype==1
	replace cq = "平衡区" if cqtype==2
	replace cq = "主销区" if cqtype==3
	drop cqtype
	
	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 rownames cq standard gini_real gini_cf ratio_gdppc_real ratio_gdppc_cf
	
	export excel using $logfiles/table_A5/table_A5.xls, firstrow(varlabels) replace
	
	
	
	
********************************************************************************
******************* 图A1.主产区和非主产区县域增产分布  **************************
********************************************************************************

	use $working_data\county_panel.dta ,clear
	
						***--- 变量处理 ---***	
	
	* 生成变量的基期(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)
	keep if e(sample)==1
	* 计算分地区的均值以便去中心化
	//主产区
	sum d_Grain if year==2016&is_zcq==1
	local mean_dgrain_zcq = `r(mean)'
	//非主产区
	sum d_Grain if year==2016&is_zcq==0
	local mean_dgrain_nzcq = `r(mean)'
	
	twoway kdensity d_Grain if year==2016&is_zcq==1,lpattern(solid) ||kdensity d_Grain if year==2016&is_zcq==0,lpattern(longdash) ||, xline(`mean_dgrain_zcq',lpattern(solid) lwidth(small) lstyle(refline)) xline(`mean_dgrain_nzcq',lpattern(longdash) lwidth(small) lstyle(refline) lcolor(gs10)) xlabel(-5(5)15) ytitle(核密度) xtitle({stSerif:{&Delta}}{it:{stSerif:Grain}} ({stSerif:十万吨})) legend(label(1 "主产区") label(2 "非主产区") ring(0) position(1) size(medsmall)) scheme(s1mono)
	graph export "$logfiles/figure_A1.png",replace width(1600) height(1200)

	
	
		
		
		
		
********************************************************************************
********************* 图A2.粮食增产量与期初粮食产量的相关性 *********************
********************************************************************************


	use $working_data\county_panel.dta ,clear

	
						***--- 变量处理 ---***	
	
	* 生成变量的基期(2003)值
	gen is_yr2003 = 1 if year==2003
	foreach v of varlist Grain{
		
		bysort county_id (year): egen `v'_03 = mean(`v'*is_yr2003)
		
	}
	
	
	* 生成相对基期的差分
	foreach v of varlist Grain{
		
		gen d_`v' = `v' - `v'_03
		
	}

	
	reg d_Grain Grain_03 if year==2016, vce(robust)
	local beta : display %6.2f _b[Grain_03]
	local se : display %6.2f _se[Grain_03]	
	local t_stat : display %6.2f _b[Grain_03]/_se[Grain_03]

	twoway scatter d_Grain Grain_03 if year==2016,mcolor(black) msize(small) || lfit d_Grain Grain_03 if year==2016, lcolor(black)||, scheme(s1mono) xtitle("{it:{stSerif:Grain{subscript:i,2003}}} ({stSerif:十万吨})") ytitle("{stSerif:{&Delta}}{it:{stSerif:Grain{subscript:i,2016}}} ({stSerif:十万吨})") legend(off) note("coef =`beta', SE =`se', t-stat =`t_stat'", size(small) width(5.7cm) ring(0) position(3) box bcolor(bg) lcolor(black) margin())
		
	graph export "$logfiles/figure_A2.png",replace width(1400) height(1200)
	
	
	
	
********************************************************************************
******************** 图A3.工具变量与内生变量的相关性 ****************************
********************************************************************************

	use $working_data\county_panel.dta ,clear

	
						***--- 变量处理 ---***	
	
	* 生成变量的基期(2003)值
	gen is_yr2003 = 1 if year==2003
	foreach v of varlist Grain{
		
		bysort county_id (year): egen `v'_03 = mean(`v'*is_yr2003)
		
	}
	
	
	* 生成相对基期的差分
	foreach v of varlist Grain{
		
		gen d_`v' = `v' - `v'_03
		
	}

						***--- 作图 ---***
						
		reg d_Grain dsim_Grain if year==2016, vce(robust)
		local beta : display %6.2f _b[dsim_Grain]
		local se : display %6.2f _se[dsim_Grain]	
		local t_stat : display %6.2f _b[dsim_Grain]/_se[dsim_Grain]
		
		twoway scatter d_Grain dsim_Grain if year==2016, mcolor(black) msize(small)|| lfit d_Grain dsim_Grain  if year==2016, lcolor(black) ||, scheme(s1mono) legend(off) ytitle("{stSerif:{&Delta}}{it:{stSerif:Grain}} ({stSerif:十万吨})") xtitle("{stSerif:{&Delta}}{it:{stSerif:Grain}}{it:{stSerif:{superscript:sim}}} ({stSerif:十万吨})") xlabel(,labsize(small)) ylabel(,labsize(small)) note("coef =`beta', SE =`se', t =`t_stat'",size(small) width(5.5cm) ring(0) position(3) box bcolor(bg) lcolor(black) margin())
		graph export "$logfiles/figure_A3.png",replace width(2000) height(1800)
	
	
	
	
	
	
*******************************************************************************
**************** 图A4.粮食增产对土地配置的影响 *********************************
*******************************************************************************



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

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


	
						***--- 逐年差分估计 ---***
						
	keep if year>=2003&year<=2016
	
	foreach v of varlist d_lncstrct_land d_lncstrct_land_new d_lnfarmland_dcrs d_lnslope_cstrct_land{
		
		statsby _b[d_Grain] _se[d_Grain],by(year) saving("$logfiles/figure_A4/est_`v'",replace): ivreghdfe `v' lnGDP_pc_03 (d_Grain = dsim_Grain), absorb(city_id) cluster(city_id)
		
		
	}
	
	

						***--- 作图输出 ---***
			
	* a.建设用地总面积
	use "$logfiles/figure_A4/est_d_lncstrct_land",clear
	rename _stat_1 est_coef
	rename _stat_2 est_se
	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   //末期=基期,设为0
	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(2003(1)2016,labsize(vsmall)) ylabel(-0.1(0.1)0.1) ytitle("估计系数", margin(small) size(medsmall)) subtitle("a.建设用地总面积",size(medlarge)) name(figure_A4_a,replace)
	
	
	
	* b.新增建设用地面积
	use "$logfiles/figure_A4/est_d_lncstrct_land_new",clear
	rename _stat_1 est_coef
	rename _stat_2 est_se
	replace est_coef = 0 if year==2003   //末期=基期,设为0
	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(2003(1)2016,labsize(vsmall)) ylabel(-0.2(0.1)0.2) ytitle("估计系数", margin(small) size(medsmall) ) subtitle("b.新开发的建设用地面积",size(medlarge)) name(figure_A4_b,replace)
	
	
	
	* c.被转换的耕地面积
	use "$logfiles/figure_A4/est_d_lnfarmland_dcrs",clear
	rename _stat_1 est_coef
	rename _stat_2 est_se
	replace est_coef = 0 if year==2003   //末期=基期,设为0
	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(2003(1)2016,labsize(vsmall)) ylabel(-0.4(0.1)0.4) ytitle("估计系数", margin(small) size(medsmall) ) subtitle("c.被转换用途的耕地面积",size(medlarge)) name(figure_A4_c,replace)
	
	
	
	* d.建设用地平均坡度
	use "$logfiles/figure_A4/est_d_lnslope_cstrct_land",clear
	rename _stat_1 est_coef
	rename _stat_2 est_se
	replace est_coef = 0 if year==2003   //末期=基期,设为0
	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(2003(1)2016,labsize(vsmall)) ylabel(-0.8(0.4)0.8) ytitle("估计系数", margin(small) size(medsmall)) subtitle("d.建设用地平均坡度",size(medlarge)) name(figure_A4_d,replace)
		
	
	* 合并
	graph combine figure_A4_a figure_A4_b figure_A4_c figure_A4_d, cols(2) scheme(s1mono) xsize(100) ysize(66)
	graph export "$logfiles/figure_A4/figure_A4.png",replace width(4000) height(2720)
	



******************************************************************************
******************** 图A5.粮食增产对企业规模与效率的影响 **********************
******************************************************************************

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

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


	
						***--- 逐年差分估计 ---***
						
	keep if year>=2003&year<=2016
	/*2003年的位置由于基期=末期，
	差分的自变量和因变量都为0，回归无意义，会出现红叉，
	后续会将其设为0，不影响结果。*/
	foreach v of varlist d_lnnum_new_firms d_share_ov1k_new_firms d_lnavgemp_old_firms d_tfp_firms{
		preserve
		//去掉数据缺失的年份
		drop if `v'==.
		statsby _b[d_Grain] _se[d_Grain],by(year) saving("$logfiles/figure_A5/est_`v'",replace): ivreghdfe `v' lnGDP_pc_03 (d_Grain = dsim_Grain), absorb(city_id) cluster(city_id)
		restore
		
	}


						***--- 作图输出 ---***
			
	* a.新注册企业数
	use "$logfiles/figure_A5/est_d_lnnum_new_firms.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(2003(1)2016,labsize(vsmall)) ylabel(-0.2(0.2)0.2) ytitle("估计系数", margin(small) size(medsmall) ) subtitle("a.进入企业数量",size(medlarge)) name(figure_A5_a,replace)
		
	
	
	* b.新注册企业中大企业占比
	use "$logfiles/figure_A5/est_d_share_ov1k_new_firms.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(2003(1)2016,labsize(vsmall)) ylabel(-0.1(0.05)0.05) ytitle("估计系数", margin(small) size(medsmall)) subtitle("b.新注册企业中大企业占比",size(medlarge)) name(figure_A5_b,replace)

	
		
	* c.在位工业企业平均就业人数
	use "$logfiles/figure_A5/est_d_lnavgemp_old_firms.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
	drop if est_coef==.  //去掉数据缺失的年份
	
	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(2003(1)2013,labsize(vsmall)) ylabel(-0.3(0.1)0.3) ytitle("估计系数", margin(small) size(medsmall) ) subtitle("c.在位工业企业平均就业人数",size(medlarge)) name(figure_A5_c,replace)
		
		
		
	* d.工业企业全要素生产率
	use "$logfiles/figure_A5/est_d_tfp_firms",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
	drop if est_coef==.  //去掉数据缺失的年份
	
	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(2003(1)2013,labsize(vsmall)) ylabel(-0.6(0.3)0.6) ytitle("估计系数", margin(small) size(medsmall) ) subtitle("d.工业企业全要素生产率",size(medlarge)) name(figure_A5_d,replace)
			
			
	* 合并
	graph combine figure_A5_a figure_A5_b figure_A5_c figure_A5_d, cols(2) scheme(s1mono) xsize(100) ysize(66)
	graph export "$logfiles/figure_A5/figure_A5.png",replace width(4000) height(2720)
		


	
	
*******************************************************************************
******************** 图A6.粮食增产对财政收入与支出的影响  **********************
*******************************************************************************
	
	use $working_data\county_panel.dta ,clear
	capture mkdir "$logfiles/figure_A6"

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


	
						***--- 逐年差分估计 ---***
						
	keep if year>=2003&year<=2016
	/*2003年的位置由于基期=末期，
	差分的自变量和因变量都为0，回归无意义，会出现红叉，
	后续会将其设为0，不影响结果。*/
	foreach v of varlist d_lntrsf_gen d_lntrsf_spec d_lnpub_expd d_lnpub_expd_agr{
		
		preserve
		drop if `v'==. //去掉数据缺失的年份
		statsby _b[d_Grain] _se[d_Grain],by(year) saving("$logfiles/figure_A6/est_`v'",replace): ivreghdfe `v' lnGDP_pc_03 (d_Grain = dsim_Grain), absorb(city_id) cluster(city_id)
		restore
		
	}
	
	

						***--- 作图输出 ---***
		
	* a.一般性转移支付收入
	use "$logfiles/figure_A6/est_d_lntrsf_gen.dta",clear
	rename _stat_1 est_coef
	rename _stat_2 est_se
	replace est_coef = 0 if year==2003  //末期=基期,设为0
	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
	drop if est_coef==.  //去掉数据缺失的年份
	
	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(2003(1)2007,labsize(vsmall)) ylabel(-0.4(0.2)0.4) ytitle("估计系数", margin(small) size(medsmall) ) subtitle("a.一般性转移支付收入",size(medlarge)) name(figure_A6_a,replace)

	
	
	
	* b.专项转移支付收入
	use "$logfiles/figure_A6/est_d_lntrsf_spec.dta",clear
	rename _stat_1 est_coef
	rename _stat_2 est_se
	replace est_coef = 0 if year==2003  //末期=基期,设为0
	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
	drop if est_coef==.  //去掉数据缺失的年份
	
	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(2003(1)2007,labsize(vsmall)) ylabel(-0.3(0.1)0.3) ytitle("估计系数", margin(small) size(medsmall) ) subtitle("b.专项转移支付收入",size(medlarge)) name(figure_A6_b,replace)

	
	
	
	
	* c.一般公共预算支出
	use "$logfiles/figure_A6/est_d_lnpub_expd.dta",clear
	rename _stat_1 est_coef
	rename _stat_2 est_se
	replace est_coef = 0 if year==2003  //末期=基期,设为0
	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
	drop if est_coef==.  //去掉数据缺失的年份
	
	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(2003(1)2007,labsize(vsmall)) ylabel(-0.1(0.1)0.1) ytitle("估计系数", margin(small) size(medsmall) ) subtitle("c.一般公共预算支出",size(medlarge)) name(figure_A6_c,replace)

	
	
	* d.农业相关公共支出
	use "$logfiles/figure_A6/est_d_lnpub_expd_agr.dta",clear
	rename _stat_1 est_coef
	rename _stat_2 est_se
	replace est_coef = 0 if year==2003  //末期=基期,设为0
	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
	drop if est_coef==.  //去掉数据缺失的年份
	
	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(2003(1)2006,labsize(vsmall)) ylabel(-0.2(0.2)0.4) ytitle("估计系数", margin(small) size(medsmall) ) subtitle("d.农业相关公共支出",size(medlarge)) name(figure_A6_d,replace)

	
			
	* 图片合并
	graph combine figure_A6_a figure_A6_b figure_A6_c figure_A6_d, cols(2) scheme(s1mono) xsize(100) ysize(66)
	graph export "$logfiles/figure_A6/figure_A6.png",replace width(2000) height(1360)
	
	



	
*******************************************************************************
************************* 表III1.工具变量排他性假设的检验 ************************
*******************************************************************************


	use $working_data\county_panel.dta ,clear

	
						***--- 变量处理 ---***	
	
	* 生成变量的基期(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-2002)
	gen is_yr1998 = 1 if year==1998
	gen is_yr2002 = 1 if year==2002
	
	bysort county_id (year): egen lnGDP_pc_98 = mean(lnGDP_pc*is_yr1998)
	bysort county_id (year): egen lnGDP_pc_02 = mean(lnGDP_pc*is_yr2002)
	
	gen d_lnGDP_pc_9802 = lnGDP_pc_02 - lnGDP_pc_98

	gl control lnGDP_pc_03 
	
	
	
	
			 ***--- 检验作物潜在产出份额对事前趋势的影响 ---***	
	
	* (1). 水稻
	
	reghdfe d_lnGDP_pc_9802 shr_potprod_rice $control if year==2016, absorb(city_id) vce(cluster city_id)
	
	outreg2 using "$logfiles/table_III1", ctitle("rice") addtext(Controls, Y, City FE, Y) keep(shr_potprod_rice) replace excel nocons depvar dec(3)
	
	
	
	* (2). 小麦
	
	reghdfe d_lnGDP_pc_9802 shr_potprod_wheat $control if year==2016, absorb(city_id) vce(cluster city_id)
	
	outreg2 using "$logfiles/table_III1", ctitle("wheat") addtext(Controls, Y, City FE, Y) keep(shr_potprod_wheat) append excel nocons depvar dec(3)
	
	
	
	* (3). 玉米
	
	reghdfe d_lnGDP_pc_9802 shr_potprod_maize $control if year==2016, absorb(city_id) vce(cluster city_id)
	
	outreg2 using "$logfiles/table_III1", ctitle("maize") addtext(Controls, Y, City FE, Y) keep(shr_potprod_maize) append excel nocons depvar dec(3)
	
	
	* (4). 大豆
	
	reghdfe d_lnGDP_pc_9802 shr_potprod_soybean $control if year==2016, absorb(city_id) vce(cluster city_id)
	
	outreg2 using "$logfiles/table_III1", ctitle("soybean") addtext(Controls, Y, City FE, Y) keep(shr_potprod_soybean) append excel nocons depvar dec(3)
	
	
	* (5). 谷子
	
	reghdfe d_lnGDP_pc_9802 shr_potprod_millet $control if year==2016, absorb(city_id) vce(cluster city_id)
	
	outreg2 using "$logfiles/table_III1", ctitle("millet") addtext(Controls, Y, City FE, Y) keep(shr_potprod_millet) append excel nocons depvar dec(3)
	
	
	* (6). 高粱
	
	reghdfe d_lnGDP_pc_9802 shr_potprod_sorghum $control if year==2016, absorb(city_id) vce(cluster city_id)
	
	outreg2 using "$logfiles/table_III1", ctitle("sorghum") addtext(Controls, Y, City FE, Y) keep(shr_potprod_sorghum) append excel nocons depvar dec(3)
	
	
	
			  ***--- 用其他IV估计 ---***	
	
	* (7). 构建IV时去除小麦和谷子
	
	ivreghdfe d_lnGDP_pc $control (d_Grain = dsim_Grain_dwm) if year==2016, absorb(city_id) cluster(city_id)
	
	outreg2 using "$logfiles/table_III1", ctitle("drop wheat&millet") addtext(Controls, Y, City FE, Y) keep(d_Grain) append excel nocons nor2 depvar dec(3) 
	
	
	
	* (8). 用所有作物的潜在产出份额作IV，进行过度识别检验
	
	global ivset = "(d_Grain = shr_potprod_rice shr_potprod_wheat shr_potprod_maize shr_potprod_soybean shr_potprod_millet shr_potprod_sorghum)"
	
	ivreghdfe d_lnGDP_pc $control $ivset if year==2016, absorb(city_id) cluster(city_id)

	outreg2 using "$logfiles/table_III1", ctitle("overID 2SLS") addtext(Controls, Y, City FE, Y) addstat(OverID p, e(jp)) adec(3) keep(d_Grain) append excel nocons nor2 depvar dec(3)

	
	
	
*******************************************************************************
*********** 表IV.粮食增产与县域经济发展：产区和地理区域的异质性 ****************
*******************************************************************************

	use $working_data\county_panel.dta ,clear

	
						***--- 变量处理 ---***	
	
	* 生成变量的基期(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
		
	}

		
			  ***--- 分产销区异质性 ---***
	
	* (1).主产区
	
	ivreghdfe d_lnGDP_pc lnGDP_pc_03 (d_Grain = dsim_Grain) if year==2016&is_zcq==1, absorb(city_id) cluster(city_id)
	
	outreg2 using "$logfiles/table_IV1", ctitle("zcq") addtext(Controls, Y, City FE, Y) replace excel nocons nor2 depvar addstat(K-P Wald F, e(rkf)) keep(d_Grain) adec(3) dec(3)
	
	
	* (2).主销区
	
	ivreghdfe d_lnGDP_pc lnGDP_pc_03 (d_Grain = dsim_Grain) if year==2016&is_zxq==1, absorb(city_id) cluster(city_id)
	
	outreg2 using "$logfiles/table_IV1", ctitle("zxq") addtext(Controls, Y, City FE, Y) append excel nocons nor2 depvar addstat(K-P Wald F, e(rkf)) keep(d_Grain) adec(3) dec(3)
	
	
	* (3).平衡区
	
	ivreghdfe d_lnGDP_pc lnGDP_pc_03 (d_Grain = dsim_Grain) if year==2016&is_phq==1, absorb(city_id) cluster(city_id)
	
	outreg2 using "$logfiles/table_IV1", ctitle("phq") addtext(Controls, Y, City FE, Y) append excel nocons nor2 depvar addstat(K-P Wald F, e(rkf)) keep(d_Grain) adec(3) dec(3)
	
	
	
			  ***--- 分地理区域异质性 ---***
			  
	* (4).东部

	ivreghdfe d_lnGDP_pc lnGDP_pc_03 (d_Grain = dsim_Grain) if year==2016&region_east==1, absorb(city_id) cluster(city_id)
	
	outreg2 using "$logfiles/table_IV1", ctitle("East") addtext(Controls, Y, City FE, Y) append excel nocons nor2 depvar addstat(K-P Wald F, e(rkf)) keep(d_Grain) adec(3) dec(3)
	
	
	* (5).中部
	
	ivreghdfe d_lnGDP_pc lnGDP_pc_03 (d_Grain = dsim_Grain) if year==2016&region_mid==1, absorb(city_id) cluster(city_id)
	
	outreg2 using "$logfiles/table_IV1", ctitle("Middle") addtext(Controls, Y, City FE, Y) append excel nocons nor2 depvar addstat(K-P Wald F, e(rkf)) keep(d_Grain) adec(3) dec(3)
	
	
	* (6).西部
	
	ivreghdfe d_lnGDP_pc lnGDP_pc_03 (d_Grain = dsim_Grain) if year==2016&region_west==1, absorb(city_id) cluster(city_id)
	
	outreg2 using "$logfiles/table_IV1", ctitle("West") addtext(Controls, Y, City FE, Y) append excel nocons nor2 depvar addstat(K-P Wald F, e(rkf)) keep(d_Grain) adec(3) dec(3)
	
	
	* (7) 东北
	
	ivreghdfe d_lnGDP_pc lnGDP_pc_03 (d_Grain = dsim_Grain) if year==2016&region_noreast==1, absorb(city_id) cluster(city_id)
	
	outreg2 using "$logfiles/table_IV1", ctitle("Northeast") addtext(Controls, Y, City FE, Y) append excel nocons nor2 depvar addstat(K-P Wald F, e(rkf)) keep(d_Grain) adec(3) dec(3)
		



	
********************************************************************************
******************** 表V1、图V1 .粮食增产与县域经济发展：稳健性检验I ***********
********************************************************************************

	use $working_data\county_panel.dta ,clear

	
	
						***--- 变量处理 ---***	
	
	* 生成匹配变量的基期(2003)值
	gen is_yr2003 = 1 if year==2003

	foreach v of varlist Grain lnGDP_pc lnRainfall lnTempr Share_PriGDP Share_SecGDP Farmland_share{
		
		bysort county_id (year): egen `v'_03 = mean(`v'*is_yr2003)
		
	}
	
	foreach v of varlist Grain lnGDP_pc{
		
		gen d_`v' = `v' - `v'_03
		
	}


						***--- 回归结果 ---***	
		
	
	* (1). 年份平均作为基期末期
	
	//基期: 2001-2003
	gen is_yr0103 = 1 if year>=2001&year<=2003
	bysort county_id (year): egen Grain_a0103 = mean(Grain*is_yr0103)
	bysort county_id (year): egen lnGDP_pc_a0103 = mean(lnGDP_pc*is_yr0103)
	bysort county_id (year): egen sim_Grain_a0103 = mean(sim_Grain*is_yr0103)
	
	
	//末期: 2014-2016
	gen is_yr1416 = 1 if year>=2014&year<=2016
	bysort county_id (year): egen Grain_a1416 = mean(Grain*is_yr1416)
	bysort county_id (year): egen lnGDP_pc_a1416 = mean(lnGDP_pc*is_yr1416)
	bysort county_id (year): egen sim_Grain_a1416 = mean(sim_Grain*is_yr1416)
	
	
	//年份平均间的差分
	gen d_Grain_avg = (Grain_a1416 - Grain_a0103)
	gen d_lnGDP_pc_avg = lnGDP_pc_a1416 - lnGDP_pc_a0103
	gen dsim_Grain_avg = sim_Grain_a1416 - sim_Grain_a0103
	
	
	//重新估计
	ivreghdfe d_lnGDP_pc_avg lnGDP_pc_a0103 (d_Grain_avg = dsim_Grain_avg) if year==2016, absorb(city_id) cluster(city_id)
	outreg2 using "$logfiles/table_V1", ctitle("average base") keep(d_Grain_avg) replace excel nocons nor2 depvar addstat(K-P Wald F stat, e(rkf)) dec(3) adec(3) addtext(LagGDP, N, Controls, Y, County FE, N, City FE, Y, City-Year FE, N) 
	
		
	
	* (2). 安慰剂
	// 生成结果变量的改革前变化(1998-2002)
	gen is_yr1998 = 1 if year==1998
	gen is_yr2002 = 1 if year==2002
	
	bysort county_id (year): egen lnGDP_pc_98 = mean(lnGDP_pc*is_yr1998)
	bysort county_id (year): egen lnGDP_pc_02 = mean(lnGDP_pc*is_yr2002)
	
	gen d_lnGDP_pc_9802 = lnGDP_pc_02 - lnGDP_pc_98
	
	ivreghdfe d_lnGDP_pc_9802 lnGDP_pc_03 (d_Grain = dsim_Grain) if year==2016, absorb(city_id) cluster(city_id)
	
	outreg2 using "$logfiles/table_V1", ctitle("placebo") keep(d_Grain) append excel nocons nor2 depvar addstat(K-P Wald F stat, e(rkf)) dec(3) adec(3) addtext(LagGDP, N, Controls, Y, County FE, N, City FE, Y, City-Year FE, N) 
	
	
	
	* (3). 控制趋势滞后项

	ivreghdfe d_lnGDP_pc lnGDP_pc_03 d_lnGDP_pc_9802 (d_Grain = dsim_Grain) if year==2016, absorb(city_id) cluster(city_id)
	
	outreg2 using "$logfiles/table_V1", ctitle("pretrends") keep(d_Grain) append excel nocons nor2 depvar addstat(K-P Wald F stat, e(rkf)) dec(3) adec(3) addtext(LagGDP, Y, Controls, Y, County FE, N, City FE, Y, City-Year FE, N) 
	
	
	
	* (4). 基准回归(保持样本相同)

	ivreghdfe d_lnGDP_pc lnGDP_pc_03 (d_Grain = dsim_Grain) if e(sample)==1, absorb(city_id) cluster(city_id)
	
	outreg2 using "$logfiles/table_V1", ctitle("baseline") keep(d_Grain) append excel nocons nor2 depvar addstat(K-P Wald F stat, e(rkf)) dec(3) adec(3) addtext(LagGDP, N, Controls, Y, County FE, N, City FE, Y, City-Year FE, N) 
	
	
	
	
	* (5). 倾向得分匹配
	* 按增产中位数分组
	egen grp_d_Grain = xtile(d_Grain) if year==2016, nquantiles(2) 
	gen is_yr2016 = 1 if year==2016
	bysort county_id: egen high_d_Grain = mean(grp_d_Grain*is_yr2016)
	replace high_d_Grain = high_d_Grain - 1
	
	* 匹配
	ivreghdfe d_lnGDP_pc lnGDP_pc_03 (d_Grain = dsim_Grain) if year==2016, absorb(city_id) vce(cluster city_id)
	
	capture drop matchv
	teffects psmatch (d_lnGDP_pc) (high_d_Grain lnGDP_pc_03 lnRainfall_03 lnTempr_03 Share_PriGDP_03 Share_SecGDP_03 Farmland_share_03 Grain_03 Slope) if e(sample)==1,vce(robust) pstolerance(1e-6) generate(matchv)
	
	* 平衡性检验(生成图#)
	tebalance summarize
	mat M = r(table) 
	coefplot matrix(M[,2]), noci xline(0) xlabel(-0.3(0.1)0.3, labsize(medsmall)) xline(-0.25 0.25, lpattern(dash)) scheme(s1mono) rename(lnGDP_pc_03 = "{stSerif:ln{it: GDP_pc}{it:{subscript:i,2003}}}" lnRainfall_03 = "{stSerif:ln{it: Rainfall}{it:{subscript:i,2003}}}" lnTempr_03 = "{stSerif:ln{it: Tempr}{it:{subscript:i,2003}}}" Share_PriGDP_03 = "{stSerif:{it: Share_PriGDP}{it:{subscript:i,2003}}}" Share_SecGDP_03 = "{stSerif:{it: Share_SecGDP}{it:{subscript:i,2003}}}" Farmland_share_03 = "{stSerif:{it: Farmland_Shr}{it:{subscript:i,2003}}}" Grain_03 = "{stSerif:{it: Grain}{it:{subscript:i,2003}}}" Slope = "{stSerif:{it: Slope}{it:{subscript:i,2003}}}") xtitle("标准化差值") ylabel(,labsize(vmall))
	graph export "$logfiles/figure_V1.png", replace width(1600) height(1400)
	
	
	ivreghdfe d_lnGDP_pc lnGDP_pc_03 (d_Grain = dsim_Grain) if matchv1!=., absorb(city_id) cluster(city_id)
	outreg2 using "$logfiles/table_V1", ctitle("PSM") keep(d_Grain) append excel nocons  nor2 depvar addstat(K-P Wald F stat, e(rkf)) dec(3) adec(3) addtext(LagGDP, N, Controls, Y, County FE, N, City FE, Y, City-Year FE, N) 
	
	
	
	
	* (6). 面板固定效应
	
	ivreghdfe lnGDP_pc i.year#c.lnGDP_pc_03 (Grain = sim_Grain) if year<=2016&year>=2003, absorb(county_id i.city_id#i.year) cluster(city_id) partial(i.year#c.lnGDP_pc_03)
	
	outreg2 using "$logfiles/table_V1", keep(Grain) append ctitle("panel fixed effects")excel nocons nor2 depvar dec(3) addstat(K-P Wald F stat, e(rkf)) adec(3) addtext(LagGDP, N, Controls, Y, County FE, Y, City FE, N, City-Year FE, Y) 
	
	

********************************************************************************
****************** 表V2 .粮食增产与县域经济发展：双重差分法 **********************
********************************************************************************


	use $working_data\county_panel.dta ,clear
	
	
						***--- 变量处理 ---***	
	
	
	* 生成分组变量基期(2003)值
	gen is_yr2003 = 1 if year==2003

	foreach v of varlist Grain Farmland_share{
		
		bysort county_id (year): egen `v'_03 = mean(`v'*is_yr2003)
		
	}
	
	* 以中位数分组
	//按期初粮食产量
	egen grain_2grp = xtile(Grain) if year==2003, nquantiles(2)
	bysort county_id: egen high_grain = mean(grain_2grp*is_yr2003)
	replace high_grain = high_grain - 1  //高组为1，低组为0
	
	//按期初耕地占比
	egen farmlandsr_2grp = xtile(Farmland_share) if year==2003, nquantiles(2)
	bysort county_id: egen high_farmlandshr = mean(farmlandsr_2grp*is_yr2003)
	replace high_farmlandshr = high_farmlandshr - 1  //高组为1，低组为0
	
	
	* post-treatment indicator
	gen post_2004 = (year>=2004)
	
	
	
						***--- 回归结果 ---***	

						
	* (1). 期初粮食产量
	
	reghdfe lnGDP_pc c.Grain_03#c.post_2004 if year>=1998&year<=2016,absorb(county_id i.city_id#i.year) vce(cluster city_id)
	
	outreg2 using "$logfiles/table_V2", replace excel nocons depvar dec(3) addtext(County FE, Y, City-Year FE, Y) 
	
	
	* (2). 期初粮食产量中位数分组(1=高,0=低)
	
	reghdfe lnGDP_pc c.high_grain#c.post_2004 if year>=1998&year<=2016,absorb(county_id i.city_id#i.year) vce(cluster city_id)
	
	outreg2 using "$logfiles/table_V2", append excel nocons depvar dec(3) addtext(County FE, Y, City-Year FE, Y) 
	
	
	* (3). 期初耕地占比
	
	reghdfe lnGDP_pc c.Farmland_share_03#c.post_2004 if year>=1998&year<=2016,absorb(county_id i.city_id#i.year) vce(cluster city_id)
	
	outreg2 using "$logfiles/table_V2", append excel nocons depvar dec(3) addtext(County FE, Y, City-Year FE, Y) 

	
	* (4). 期初耕地占比中位数分组(1=高,0=低)
	
	reghdfe lnGDP_pc c.high_farmlandshr#c.post_2004 if year>=1998&year<=2016,absorb(county_id i.city_id#i.year) vce(cluster city_id)
	
	outreg2 using "$logfiles/table_V2", append excel nocons depvar dec(3) addtext(County FE, Y, City-Year FE, Y)

	
	

*******************************************************************************
******************* 表 V3.粮食增产与县域经济发展：稳健性检验II ******************
*******************************************************************************

	use $working_data\county_panel.dta ,clear
	gen lnpop = ln(pop)
	
						***--- 变量处理 ---***	
	
	* 生成变量基期(2003)值
	gen is_yr2003 = 1 if year==2003

	foreach v of varlist Grain lnGDP_pc lnGDP lnLight lnPDen Sown_area lnpop{
		
		bysort county_id (year): egen `v'_03 = mean(`v'*is_yr2003)
		
	}
	
	foreach v of varlist Grain lnGDP_pc lnGDP lnLight lnPDen Sown_area lnpop{
		
		gen d_`v' = `v' - `v'_03
		
	}



						***--- 回归结果 ---***	

	* (1). 总量GDP
	
	ivreghdfe d_lnGDP lnGDP_pc_03 (d_Grain = dsim_Grain) if year==2016, absorb(city_id) cluster(city_id)
	
	outreg2 using "$logfiles/table_V3", ctitle("GDP") addtext(Controls, Y, City FE, Y)  replace excel nocons nor2 depvar dec(3) keep(d_Grain) addstat(K-P Wald F, e(rkf)) adec(3)
	
	
	* (2). 夜间灯光
	
	ivreghdfe d_lnLight lnLight_03 (d_Grain = dsim_Grain) if year==2016, absorb(city_id) cluster(city_id)
	
	outreg2 using "$logfiles/table_V3", ctitle("light") addtext(Controls, Y, City FE, Y) append excel nocons  nor2 depvar addstat(K-P Wald F, e(rkf)) keep(d_Grain) dec(3) adec(3)
	
	
	* (3). 人口密度
	
	ivreghdfe d_lnPDen lnPDen_03 (d_Grain = dsim_Grain) if year==2016, absorb(city_id) cluster(city_id)
	
	outreg2 using "$logfiles/table_V3", ctitle("population density") addtext(Controls, Y, City FE, Y) append excel nocons nor2 depvar addstat(K-P Wald F, e(rkf)) keep(d_Grain) dec(3) adec(3)
	 
	 
	* (4). 播种面积变化作为自变量
	
	ivreghdfe d_lnGDP_pc lnGDP_pc_03 (d_Sown_area = dsim_Grain) if year==2016, absorb(city_id) cluster(city_id)
	
	outreg2 using "$logfiles/table_V3", ctitle("sownarea") addtext(Controls, Y, City FE, Y) append excel nocons nor2 depvar addstat(K-P Wald F, e(rkf)) keep(d_Sown_area) dec(3) adec(3)
	
	
	* (5). 增加控制人口总量变化
	ivreghdfe d_lnGDP lnGDP_pc_03 d_lnpop (d_Grain = dsim_Grain) if year==2016, absorb(city_id) cluster(city_id)
	
	outreg2 using "$logfiles/table_V3", ctitle("Pop") addtext(Controls, Y, City FE, Y)  append excel nocons nor2 depvar dec(3) keep(d_Grain) addstat(K-P Wald F, e(rkf)) adec(3)
	
	
		

*******************************************************************************
********************** 表V4.对增产影响的替代性假说 *****************************
*******************************************************************************

	use $working_data\county_panel.dta ,clear
	
	
						***--- 变量处理 ---***	
	
	* 生成基期(2003)值
	gen is_yr2003 = 1 if year==2003

	foreach v of varlist Grain lnGDP_pc lnGDP Meat Oil confucious_temple soe_empshare{
		
		bysort county_id (year): egen `v'_03 = mean(`v'*is_yr2003)
		
	}
	
	foreach v of varlist Grain lnGDP_pc lnGDP Meat Oil{
		
		gen d_`v' = `v' - `v'_03
		
	}

	
		
						***--- 回归结果 ---***	
		
		
	* (1).2001年县平均关税
	
	ivreghdfe d_lnGDP_pc tarif_shock lnGDP_pc_03 (d_Grain = dsim_Grain) if year==2016, absorb(city_id) cluster(city_id)
	
	outreg2 using "$logfiles/table_V4", ctitle("trade:initial tariff") addtext(Controls, Y, City FE, Y) replace excel nocons nor2 depvar addstat(K-P Wald F, e(rkf)) adec(3) dec(3) keep(d_Grain)
	
	* (2). IV构建去除大豆
	
	ivreghdfe d_lnGDP_pc lnGDP_pc_03 (d_Grain = dsim_Grain_nosy) if year==2016, absorb(city_id) cluster(city_id)
	
	outreg2 using "$logfiles/table_V4", ctitle("trade:no soybeans") addtext(Controls, Y, City FE, Y) append excel nocons nor2 depvar addstat(K-P Wald F, e(rkf)) adec(3) dec(3) keep(d_Grain)

	
	* (3). 控制油料和肉类生产
	
	ivreghdfe d_lnGDP_pc d_Meat d_Oil lnGDP_pc_03 (d_Grain = dsim_Grain) if year==2016, absorb(city_id) cluster(city_id)
	
	outreg2 using "$logfiles/table_V4", ctitle("other agr. products") addtext(Controls, Y, City FE, Y) append excel nocons nor2 depvar addstat(K-P Wald F, e(rkf)) adec(3) dec(3) keep(d_Grain)
	
	
	* (4). 稻作文化: GAEZ的水稻潜在单产
	
	ivreghdfe d_lnGDP_pc potyield_rice lnGDP_pc_03 (d_Grain = dsim_Grain) if year==2016,absorb(city_id) cluster(city_id)
	
	outreg2 using "$logfiles/table_V4", ctitle("Rice-farming") addtext(Controls, Y, City FE, Y) append excel nocons  nor2 depvar addstat(K-P Wald F, e(rkf)) adec(3) dec(3) keep(d_Grain)
	
	
	* (5). 儒家文化
	
	ivreghdfe d_lnGDP_pc confucious_temple_03 lnGDP_pc_03 (d_Grain = dsim_Grain) if year==2016,absorb(city_id) cluster(city_id)
	
	outreg2 using "$logfiles/table_V4", ctitle("Confucius") addtext(Controls, Y, City FE, Y) append excel nocons  nor2 depvar addstat(K-P Wald F, e(rkf)) adec(3) dec(3) keep(d_Grain)
	
	
	* (6). 市场化：控制期初国企占比
	ivreghdfe d_lnGDP_pc soe_empshare_03 lnGDP_pc_03 (d_Grain = dsim_Grain) if year==2016,absorb(city_id) cluster(city_id)
	outreg2 using "$logfiles/table_V4", ctitle("SOE share") addtext(Controls, Y, City FE, Y) append excel nocons  nor2 depvar addstat(K-P Wald F, e(rkf)) adec(3) dec(3) keep(d_Grain)			
		
				
	* (7): 去掉东北地区县
	
	ivreghdfe d_lnGDP_pc lnGDP_pc_03 (d_Grain = dsim_Grain) if year==2016&region_noreast==0, absorb(city_id) cluster(city_id)
	
	outreg2 using "$logfiles/table_V4", ctitle("no northeastern") addtext(Controls, Y, City FE, Y) append excel nocons nor2 depvar addstat(K-P Wald F, e(rkf)) adec(3) dec(3) keep(d_Grain)
		
	
	* (8). 去掉沿海县
	
	ivreghdfe d_lnGDP_pc lnGDP_pc_03 (d_Grain = dsim_Grain) if year==2016&is_coastal==0, absorb(city_id) cluster(city_id)
	
	outreg2 using "$logfiles/table_V4", ctitle("no seaside county") addtext(Controls, Y, City FE, Y) append excel nocons nor2 depvar addstat(K-P Wald F, e(rkf)) adec(3) dec(3) keep(d_Grain)


	
	
	
*******************************************************************************
************************* 图VI1.粮食增产对土地财政的影响 ************************
*******************************************************************************

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

	
						***--- 变量处理 ---***	
	
	* 生成变量的基期(2003)值
	gen is_yr2003 = 1 if year==2003
	foreach v of varlist Grain lnGDP_pc lnlandinc lnlandexp{
		
		bysort county_id (year): egen `v'_03 = mean(`v'*is_yr2003)
		
	}
	
	
	* 生成相对基期的差分
	foreach v of varlist Grain lnlandinc lnlandexp{
		
		gen d_`v' = `v' - `v'_03
		
	}


						***--- 逐年差分估计 ---***
	/*2003年的位置由于基期=末期，
	差分的自变量和因变量都为0，回归无意义，会出现红叉，
	后续会将其设为0，不影响结果。*/
	foreach v of varlist d_lnlandinc d_lnlandexp{
		preserve
		drop if `v'==. //去掉数据缺失的年份
		statsby _b[d_Grain] _se[d_Grain],by(year) saving("$logfiles/figure_VI1/est_`v'",replace): ivreghdfe `v' lnGDP_pc_03 (d_Grain = dsim_Grain) if year>=2003&year<=2016, absorb(city_id) cluster(city_id)
	    restore
	}
		
		
						***--- 作图输出 ---***
						
		* a. 土地相关财政收入
		use "$logfiles/figure_VI1/est_d_lnlandinc",clear
		rename _stat_1 est_coef
		rename _stat_2 est_se
		replace est_coef = 0 if year==2003      //末期=基期，设为0
		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
		drop if 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(2003(1)2007,labsize(vsmall)) ylabel(-0.8(0.4)0.8) ytitle("估计系数", margin(small) size(medsmall)) subtitle("a.土地相关财政收入",size(large))  name(figure_VI1_a,replace)

		
		
		* b. 土地相关财政支出
		use "$logfiles/figure_VI1/est_d_lnlandexp",clear
		rename _stat_1 est_coef
		rename _stat_2 est_se
		replace est_coef = 0 if year==2003      //末期=基期，设为0
		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
		drop if 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(2003(1)2007,labsize(vsmall)) ylabel(-0.8(0.4)0.8) ytitle("估计系数", margin(small) size(medsmall) ) subtitle("b.土地相关财政支出",size(large)) name(figure_VI1_b,replace)

		
	graph combine figure_VI1_a figure_VI1_b, cols(2) scheme(s1mono) xsize(100) ysize(50)
	graph export "$logfiles/figure_VI1/figure_VI1.png",replace width(2500) height(1200)
		

		
		
		
		
		log close

	