*Analysis file of ACS data for NYBS DOH Wage Study


use "${clean_data}\acs2017_clean.dta", clear

keep if nurse==1


*Gen annual hours
	gen ann_hrs = wkswork1*uhrswork
	label var ann_hrs "Annual Hours"
	*Breakdown of Annual Hours
	forvalues i=1(1)3 {
		if `i'==1  local z "Registered Nurses"
		if `i'==2  local z "Licensed Practical Nurses"
		if `i'==3  local z "Nursing Assistants and Other Aides"
		histogram ann_hr if xn1==`i' & nys==1, title("Distribution of Annual Hours") subtitle("Among `z'")
		gr export  "${output}\histo_xn`i'_annhr.pdf", replace
		}
	*There are some individuals reporting extremely high hours/wk. I am going to cut those with hr/wk in the 95th percentile. Cut point (50 hr/wk) is same in NYS and full pop. 
	sum uhrswork if nys==1, d
	sum uhrswork, d
	*drop if uhrswork>50
	winsor2 hrwage, replace cuts(5 95) trim
	winsor2 ann_hr, replace cuts(5 95) trim





*Prepare final analysis variables
gen ln_ann_hrs = ln(ann_hrs)
gen work = uhrswork>0
gen outside_inc = (ftotinc - incwage)/1000
replace outside_inc = 0 if incwage>=999998
replace outside_inc = 0  if ftotinc>=999998
gen female = sex==2
gen married = marst==1 | marst==2
tab empstat, gen(xemp_)
gen otherind = xind==3
gen pnchlt5 = nchlt5>0
gen pnchgt5 = nchgt5>0
gen pt = uhrs<35
label var ln_ann_hrs "Log Annual Hours"
label var xemp_1 "Employed"
label var xemp_2 "Unemployed"
label var married "Married"
label var female "Female"
label var outside_inc "Household Income Less Own Earnings"
label var work "Positive Hours"
label var age "Age"
label var black "Black"
label var other "Other"
label var hisp "Hispanic"
label var white "White"
label var nchlt5 "Children under 5"
label var nchgt5 "Children over 5"
label var pnchlt5 "Has children under 5"
label var pnchgt5 "Has children over 5"
label var metro "Metropolitan"
label var hosp "Hospitals"
label var skf "Skilled Nursing Facilities"
label var otherind "Other"
label var edu_lths "Less than High School" 
label var edu_hsch "High School" 
label var edu_scoll "Some College"
label var edu_coll "College"
label var pt "Part-Time"




*Breakdown of Hourly Wage for 3 nurse types, states, and industry
preserve
	keep if work==1
	drop if xn1==.
	collapse (mean) hrwage [pweight=perwt], by(xn1 xind xst)
	export excel using "${output}\summ_hr_wage_states.xls", firstrow(varlabels) replace
	restore
	
*Breakdown of Hourly Wage for 3 nurse types,  industry WITHIN NYS
preserve
	keep if work==1
	drop if xn1==.
	keep if nys==1
	collapse (mean) hrwage [pweight=perwt], by(xn1 xind ny_edr)
	export excel using "${output}\summ_hr_wage_nys.xls", firstrow(varlabels) replace
	restore
	
	
	

	
*******************************************************************************
*******************************************************************************
***		Remainder of analysis only looks within NYS							***
*******************************************************************************
*******************************************************************************	
keep if nys==1



*Summary of variables
local var work xemp_1 xemp_2 ann_hrs hosp skf otherind 
eststo sumstatsA1: estpost tabstat `var'   if xn1==1, statistics(mean sd) columns(statistics) 
eststo sumstatsA2: estpost tabstat `var'   if xn1==2, statistics(mean sd) columns(statistics) 
eststo sumstatsA3: estpost tabstat `var'   if xn1==3, statistics(mean sd) columns(statistics) 
esttab sumstatsA1 sumstatsA2 sumstatsA3 using "${output}\summ_t1.csv" , main(mean) aux(sd) label nostar unstack nogaps nonotes mtitles("Registered Nurses" "Licensed Practical Nurses" "Nursing Assistants and Other Aides")  replace


local var edu_lths edu_hsch edu_scoll edu_coll age female black other hisp white married pnchlt5 pnchgt5 metro
eststo sumstatsA1: estpost tabstat `var'   if xn1==1, statistics(mean sd) columns(statistics) 
eststo sumstatsA2: estpost tabstat `var'   if xn1==2, statistics(mean sd) columns(statistics) 
eststo sumstatsA3: estpost tabstat `var'   if xn1==3, statistics(mean sd) columns(statistics) 
esttab sumstatsA1 sumstatsA2 sumstatsA3 using "${output}\summ_t2.csv", main(mean) aux(sd) label nostar unstack nogaps nonotes mtitles("Registered Nurses" "Licensed Practical Nurses" "Nursing Assistants and Other Aides") replace



	
*Participation Analysis for each nurse type 
local grab
forvalues i=1(1)3 {
	reg work ln_hrwage i.xind c.ln_hrwage#i.xind outside_inc i.sex age c.age#c.age i.xedu black other hisp  nchlt5 nchgt5 i.marst i.metro i.ny_edr i.multyear if xn1==`i' , robust
	eststo m`i'_a: margins, dydx(ln_hrwage) at(xind=(1 2 3))  atmeans post
	probit work ln_hrwage i.xind c.ln_hrwage#i.xind outside_inc i.sex age c.age#c.age i.xedu black other hisp  nchlt5 nchgt5 i.marst i.metro i.ny_edr i.multyear if xn1==`i' , robust
	eststo m`i'_b: margins, dydx(ln_hrwage) at(xind=(1 2 3))  atmeans post
	local grab `grab' m`i'_a m`i'_b
	}
esttab `grab' using "${output}\participation_e.csv",  replace

local grab
forvalues i=1(1)3 {
	reg work ln_hrwage i.xind c.ln_hrwage#i.xind outside_inc i.sex age c.age#c.age i.xedu black other hisp  nchlt5 nchgt5 i.marst i.metro i.ny_edr i.multyear if xn1==`i' , robust
	eststo m`i'_a: margins, dydx(ln_hrwage)   atmeans post
	probit work ln_hrwage i.xind c.ln_hrwage#i.xind outside_inc i.sex age c.age#c.age i.xedu black other hisp  nchlt5 nchgt5 i.marst i.metro i.ny_edr i.multyear if xn1==`i' , robust
	eststo m`i'_b: margins, dydx(ln_hrwage)  atmeans post
	local grab `grab' m`i'_a m`i'_b
	}
esttab `grab' using "${output}\participation_total_e.csv",  replace


*Hours Elasticity
forvalues nr=1(1)3 {
	reg ln_ann_hrs ln_hrwage i.xind c.ln_hrwage#i.xind outside_inc i.sex age c.age#c.age i.xedu black other hisp  nchlt5 nchgt5 i.marst i.metro i.ny_edr i.multyear if xn1==`nr', robust
	eststo h`nr': margins, dydx(ln_hrwage) at(xind=(1 2 3)) atmeans vsquish post
	}
esttab h1 h2 h3 using "${output}\hours_e.csv",  replace

forvalues nr=1(1)3 {
	reg ln_ann_hrs ln_hrwage i.xind c.ln_hrwage#i.xind outside_inc i.sex age c.age#c.age i.xedu black other hisp  nchlt5 nchgt5 i.marst i.metro i.ny_edr i.multyear if xn1==`nr', robust
	eststo h`nr': margins, dydx(ln_hrwage)  atmeans vsquish post
	}
esttab h1 h2 h3 using "${output}\hours_total_e.csv",  replace


/*
*Hours Elasticity by PT FT
forvalues nr=1(1)3 {
	reg ln_ann_hrs ln_hrwage i.xind c.ln_hrwage#i.xind outside_inc i.sex age c.age#c.age i.xedu black other hisp  nchlt5 nchgt5 i.marst i.metro i.ny_edr i.multyear if xn1==`nr', robust
	eststo h`nr': margins, dydx(ln_hrwage) at(xind=(1 2 3)) atmeans vsquish post
	}
esttab h1 h2 h3 using "${output}\hours_e.csv",  replace
*/

* DO I LOOK AT ELASTICITIES FOR THOSE WHO ARE WORKING LESS THAN FY-FT?
forvalues nr=1(1)3 {
	reg ln_ann_hrs ln_hrwage i.xind c.ln_hrwage#i.xind outside_inc i.sex age c.age#c.age i.xedu black other hisp  nchlt5 nchgt5 i.marst i.metro i.ny_edr i.multyear if xn1==`nr' & pt==1, robust
	eststo h`nr': margins, dydx(ln_hrwage) at(xind=(1 2 3)) atmeans vsquish post
	}
esttab h1 h2 h3 using "${output}\hours_e_pt.csv",  replace

forvalues nr=1(1)3 {
	reg ln_ann_hrs ln_hrwage i.xind c.ln_hrwage#i.xind outside_inc i.sex age c.age#c.age i.xedu black other hisp  nchlt5 nchgt5 i.marst i.metro i.ny_edr i.multyear if xn1==`nr' & pt==0, robust
	eststo h`nr': margins, dydx(ln_hrwage) at(xind=(1 2 3)) atmeans vsquish post
	}
esttab h1 h2 h3 using "${output}\hours_e_ft.csv",  replace