clear all
sysdir
adopath
cd "D:\Data for paper"   

/*==================================================
              1. CPI
==================================================*/                  
** 1.1. Annual CPI by Province
import excel "D:\Data for paper\居民消费价格指数.xls", sheet("分省年度数据") firstrow clear
destring 时间, i(年) replace
gather 北京市 - 新疆维吾尔自治区, variable(province) value(index)
rename 时间 year
keep if inrange(year,2014,2020)
summ index
sort province year
gen CPI = 1 if year==2014
replace CPI = CPI[_n-1] * (index[_n]/100) if inrange(year,2015,2020)  
keep if inlist(year,2014,2016,2018,2020)
summ CPI

** 1.2. Province Name
replace province = usubstr(province,1,2)
replace province = "黑龙江" if province=="黑龙"
replace province = "内蒙古" if province=="内蒙"

** 1.3. Match province code
merge m:1 province using "省份编码对照"
drop _merge
rename province provname
rename pcode province_m
rename CPI CPI_m
keep year province_m CPI_m
save CPI_m.dta ,replace
rename province_m province_f
rename CPI_m CPI_f
save CPI_f.dta ,replace


/*==================================================
              2. Consolidate into panel data
==================================================*/    
**2.1. Vertical merger
use 2020mo_fa_child_2.dta,clear
append using 2018mo_fa_child_2.dta
append using 2016mo_fa_child_2.dta
append using 2014mo_fa_child_2.dta
save 2014-2020mo_fa_child_2.dta, replace

use 2014-2020mo_fa_child_2.dta,clear
drop num*
bys pid_m: egen num2 = count(pid_m) 
tab num2
xtset pid_m year

**2.2. Horizontal merger
tab province_m   
replace province_m=province2_m if province_m==-9 | province_m==.  

merge  m:1 year province_m using CPI_m.dta
list pid_m year province_m province2_m if _merge==1
keep if _merge!=2
rename _merge _merge_p_m

tab _merge_f      
tab province_f  
replace province_f=province2_f if province_f==-9 | province_f==.
merge  m:1 year province_f using CPI_f.dta  
keep if _merge!=2
rename _merge _merge_p_f

/*==================================================
              3. Variables
==================================================*/    
**3.1.Mother's or Father's Earing (adjusted for CPI)
gen inc_m=income_m
replace inc_m=emp_income_m if year==2020
tab inc_m
replace inc_m=. if inc_m<0
gen inc_cpi_m=inc_m/CPI_m
tab inc_cpi_m
gen lninc_cpi_m=ln(inc_cpi_m+1)
summ lninc_cpi_m

gen inc_f=income_f
replace inc_f=emp_income_f if year==2020
tab inc_f
replace inc_f=. if inc_f<0
gen inc_cpi_f=inc_f/CPI_f
tab inc_cpi_f
gen lninc_cpi_f=ln(inc_cpi_f+1)
summ lninc_cpi_f

**3.2. Child-Related Emotional Labor Context
gen emo_demand=e_5

**3.3. living with parents or parents-in-law
gen grandparents_live=0
replace grandparents_live=1 if par_f_reside_m==1 | par_m_reside_m==1 |par_f_reside_f==1|  par_m_reside_f==1
tab grandparents_live

**3.4. health status
tab health_m
recode health_m (1 2 3 = 1 "Healthy")(4 = 2 "Average")(5 = 3 "Unhealthy"), gen(healthy_m )
tab healthy_m
tab health_f
tabstat health_f,stat(mean n) by(health_f) 
recode health_f (1 2 3 = 1 "Healthy")(4 = 2 "Average")(5 = 3 "Unhealthy"), gen(healthy_f )
tab healthy_f

**3.5. region (using the western region as the reference group)
gen east_m=0
label variable east_m "Eastern region"
replace east_m=1 if (province_m==11 | province_m==12 | province_m==13 |province_m==21| province_m==31 | province_m==32 | province_m==33 | province_m==35 | province_m==37 | province_m==44 | province_m==46) 
tab east_m
gen middle_m=0
replace middle_m=1 if (province_m==14  | province_m==22 | province_m==23 | province_m==34 | province_m==36 | province_m==41 | province_m==42 | province_m==43)  
tab middle_m
label variable middle_m "Central region"

gen east_f=0
label variable east_f "Eastern region"
replace east_f=1 if (province_f==11 | province_f==12 | province_f==13 |province_f==21| province_f==31 | province_f==32 | province_f==33 | province_f==35 | province_f==37 | province_f==44 | province_f==46) 
tab east_f
gen middle_f=0
replace middle_f=1 if (province_f==14  | province_f==22 | province_f==23 | province_f==34 | province_f==36 | province_f==41 | province_f==42 | province_f==43)  
tab middle_f
label variable middle_f "Central region"

**3.6. emotional intelligence
gen emointel_m=popularity_m
gen emointel_f=popularity_f

/*==================================================
         4. Remove observations from the sample
==================================================*/    
**4.1. Exclude self-employed observations
keep if job_type_m>2  
tab lninc_cpi_m if job_type_m==.
replace job_type_m=4 if job_type_m==. 
tab job_type_m

**4.2. Exclude observations with missing values for the key variables
keep if lninc_cpi_m!=.  
tab year, sum(lninc_cpi_m)
egen lninc_cpi_m_mean=mean(lninc_cpi_m), by(year)
tab lninc_cpi_m_mean

keep if emo_demand !=.   

tab _merge_f  
keep if _merge_f==3  

**4.2. Exclude observations with missing values for the control variables
tab edu_m
replace edu_m=edu2_m if edu_m==.

tab urban_m
replace urban_m=urban2_m if urban_m==.

tab edu_f
replace edu_f=edu2_f if edu_f==.

egen miss1 = rowmiss(age_m edu_m childnumber2_m health_m internet_m urban_m c0_2_m c3_5_m east_m middle_m ) 
tab miss1
keep if miss1 == 0

egen miss3 = rowmiss(edu_f health_f)
tab miss3
keep if miss3 == 0
tab year, sum(lninc_cpi_m) 
tab year, sum(lninc_cpi_f) 

/*==================================================
         5. Descriptive Statistics
==================================================*/    
**5.1. Mother Sample 
tabstat lninc_cpi_m emo_demand age_m age2_m edu_m healthy_m internet_m urban_m childnumber2_m c0_2_m c3_5_m  edu_f healthy_f grandparents_live east_m middle_m housework_m whour_m admin_m contract_m emointel_m emo_m, stats(N mean sd ) col(s) format(%9.4f)  

**5.2. Father Sample 
tabstat lninc_cpi_f emo_demand age_f age2_f edu_f healthy_f internet_f urban_f childnumber2_f c0_2_f c3_5_f  edu_m healthy_m  grandparents_live east_f middle_f housework_f whour_f admin_f contract_f emointel_f emo_f if job_type_f>2 , stats(N mean sd ) col(s) format(%9.4f)  

/*==================================================
             6. Regression Analysis
==================================================*/
**6.1. Baseline  Regression
xtreg lninc_cpi_m emo_demand i.year,fe r 
est store m_1

xtreg lninc_cpi_f emo_demand i.year if job_type_f>2,fe r 
est store m_2

xtreg lninc_cpi_m emo_demand age_m age2_m edu_m healthy_m internet_m urban_m childnumber2_m c0_2_m c3_5_m edu_f healthy_f grandparents_live east_m middle_m i.year,fe r 
est store m_3

xtreg lninc_cpi_f emo_demand age_f age2_f edu_f healthy_f internet_f urban_f childnumber2_f c0_2_f c3_5_f edu_m healthy_m grandparents_live east_f middle_f i.year if job_type_f>2,fe r 
est store m_4

esttab  m_1 m_2 m_3 m_4 , scalar(r2 r2_a N) compress  ///
                star(* 0.1 ** 0.05 *** 0.01)  ///
				b(%6.4f) se(%6.4f)

**6.2. Mechanism Analysis
**6.2.1. Working Hours
replace housework_m=housework_original_m if year==2016 | year==2018
tab housework_m 
xtreg housework_m emo_demand age_m age2_m edu_m healthy_m internet_m urban_m childnumber2_m c0_2_m c3_5_m  edu_f healthy_f grandparents_live east_m middle_m  i.year,fe r  
est store m_5

xtreg  lninc_cpi_m housework_m emo_demand age_m age2_m edu_m childnumber2_m healthy_m internet_m urban_m c0_2_m c3_5_m  edu_f healthy_f grandparents_live east_m middle_m i.year,fe r  
est store m_6

**6.2.2. Working Hours
tab whour_m
xtreg whour_m emo_demand age_m age2_m edu_m childnumber2_m healthy_m internet_m urban_m c0_2_m c3_5_m  edu_f healthy_f  grandparents_live east_m middle_m  i.year,fe r 
est store m_7

xtreg lninc_cpi_m whour_m emo_demand age_m age2_m edu_m childnumber2_m healthy_m internet_m urban_m c0_2_m c3_5_m  edu_f healthy_f grandparents_live east_m middle_m i.year,fe r 
est store m_8

esttab  m_5 m_6 m_7 m_8, scalar(r2 r2_a N) compress  ///
                star(* 0.1 ** 0.05 *** 0.01)  ///
				b(%6.4f) se(%6.4f)
 
**6.3. Heterogeneity Analysis
**6.3.1.Emotional Intelligence
xtreg lninc_cpi_m emo_demand age_m age2_m edu_m healthy_m internet_m urban_m childnumber2_m c0_2_m c3_5_m  edu_f healthy_f grandparents_live east_m middle_m i.year if (emointel_m<6),fe r  
est store m_9

xtreg lninc_cpi_m emo_demand age_m age2_m edu_m healthy_m internet_m urban_m childnumber2_m c0_2_m c3_5_m  edu_f healthy_f grandparents_live east_m middle_m i.year if (emointel_m>5),fe r 
est store m_10

**6.3.2. Emotional Disposition
tab emo_m
xtreg lninc_cpi_m emo_demand age_m age2_m edu_m health_m internet_m urban_m childnumber2_m c0_2_m c3_5_m edu_f healthy_f grandparents_live east_m middle_m i.year  if (emo_m==1 | emo_m==2 ),fe r 
 est store m_11
 
xtreg lninc_cpi_m emo_demand age_m age2_m edu_m health_m internet_m urban_m childnumber2_m c0_2_m c3_5_m edu_f healthy_f grandparents_live east_m middle_m i.year  if (emo_m==4 | emo_m==3 ),fe r   
est store m_12

**6.3.3. Educational Attainment
tab edu_m
xtreg lninc_cpi_m emo_demand age_m age2_m healthy_m internet_m urban_m childnumber2_m c0_2_m c3_5_m  edu_f healthy_f grandparents_live east_m middle_m i.year if edu_m<4,fe r 
est store m_13

xtreg lninc_cpi_m emo_demand age_m age2_m healthy_m internet_m urban_m childnumber2_m c0_2_m c3_5_m  edu_f healthy_f grandparents_live east_m middle_m i.year if edu_m>3,fe r 
est store m_14

 esttab  m_9 m_10 m_11 m_12 m_13 m_14, scalar(r2 r2_a N) compress  ///
                star(* 0.1 ** 0.05 *** 0.01)  ///
				b(%6.4f) se(%6.4f)

				
**6.3.4. Employment Stability 
tab contract_m
xtreg lninc_cpi_m emo_demand age_m age2_m edu_m healthy_m internet_m urban_m childnumber2_m c0_2_m c3_5_m  edu_f healthy_f grandparents_live east_m middle_m i.year if contract_m==0,fe r  
est store m_15
xtreg lninc_cpi_m emo_demand age_m age2_m edu_m healthy_m internet_m urban_m childnumber2_m c0_2_m c3_5_m  edu_f healthy_f grandparents_live east_m middle_m i.year if contract_m==1,fe r 
est store m_16

**6.3.4. Occupational status 
tab admin_m
xtreg lninc_cpi_m emo_demand age_m age2_m edu_m healthy_m internet_m urban_m childnumber2_m c0_2_m c3_5_m  edu_f healthy_f grandparents_live east_m middle_m i.year if admin_m==0,fe r 
est store m_17
xtreg lninc_cpi_m emo_demand age_m age2_m edu_m healthy_m internet_m urban_m childnumber2_m c0_2_m c3_5_m  edu_f healthy_f grandparents_live east_m middle_m i.year if admin_m==1,fe r
est store m_18

 esttab  m_15 m_16 m_17 m_18, scalar(r2 r2_a N) compress  ///
                star(* 0.1 ** 0.05 *** 0.01)  ///
				b(%6.4f) se(%6.4f)