/*** Data Analysis for Labor Supply and Income Inequality Paper ***/

/*** Analysis of ORG Data
1.  Graph of share of men and women who work >=49 or <=35 loghours per week.
2.  Graph of mean weekly work loghours for men and women.
3.  Table of Share men work >=49 or <=35, mean loghours, unemployment by MSA inequality quantile
4.  Table of Share men work >=49 or <=35, mean loghours, unemployment, worksd by (educationd,aged) cells
5.  Correlation commincsd, worksd, occsd, uhourse, communemp, workunemp, occunemp
***/

clear
capture log close
cd ~/Documents/dissertation/labor_supply
set mem 750m
log using org_analysis_ls, text replace
use org_ls, clear

collapse (mean) uhourse work49 work35, by(year sex)

tsset sex year, yearly

twoway (tsline work49 if sex==1, lpattern(solid) lcolor(black) ttick(1981(4)2005) tlabel(1981(4)2005))(tsline work35 if sex==1, lpattern(dash) lcolor(black) ttick(1981(4)2005) tlabel(1981(4)2005))(tsline work49 if sex==2, lpattern(shortdash_dot)  lcolor(black) ttick(1981(4)2005) tlabel(1981(4)2005))(tsline work35 if sex==2, lpattern(dot) lcolor(black) ttick(1981(4)2005) tlabel(1981(4)2005)), legend(position(9) label(1 "Men>=49") label(2 "Men<=35") label(3 "Women>=49") label(4 "Women<=35") cols(1) symxsize(5)) xtitle("Year") ytitle("Share") scheme(s1mono)

graph export props_all.png, as(png) replace

twoway (tsline uhourse if sex==1, lpattern(solid)  lcolor(black) ttick(1981(4)2005) tlabel(1981(4)2005))(tsline uhourse if sex==2, lpattern(dash) lcolor(black) ttick(1981(4)2005) tlabel(1981(4)2005)), legend(position(9) label(1 "Men") label(2 "Women") cols(1) symxsize(5)) xtitle("Year") ytitle("Mean loghours") scheme(s1mono)

graph export mean_loghours.png, as(png) replace


use org_ls, clear
collapse (mean) work49 work35 uhourse communemp if sex==1, by(commincsdq)
outsheet using communsdq.csv, comma replace

use org_ls, clear
collapse (mean) work49 work35 uhourse workunemp worksd if sex==1, by(educationd aged)
outsheet using work.csv, comma replace

use org_ls, clear
quantiles occsd, gen(occsdq) n(10)
collapse (mean) work49 work35 uhourse occunemp occsd if sex==1, by(occsdq)
outsheet using occ.csv, comma replace

use org_ls, clear
pwcorr uhourse wage commincsd worksd occsd communemp workunemp occunemp

use org_ls, clear

collapse (p90) p90=wage (p75) p75=wage (p50) p50=wage (p25) p25=wage (p10) p10=wage, by(year)
gen ninety_ten = p90/p10
gen ninety_fifty = p90/p50
gen fifty_ten = p50/p10
gen seventy_twenty = p75/p25
tsset year, yearly

#delimit ;
twoway (tsline ninety_ten, lpattern(solid) lcolor(black) ttick(1981(4)2005) tlabel(1981(4)2005))
(tsline ninety_fifty, lpattern(dash) lcolor(black) ttick(1981(4)2005) tlabel(1981(4)2005))
(tsline fifty_ten, lpattern(dash_dot) lcolor(black) ttick(1981(4)2005) tlabel(1981(4)2005))
, legend(position(9) label(1 "90/10") label(2 "90/50") label (3 "50/10") cols(1) symxsize(5)) xtitle("Year") ytitle("Ratio") scheme(s1mono);

#delimit cr

use org_ls, clear
gen grade=.
replace grade=gradeat if gradeat!=.
replace grade=grade92 if grade92!=.
collapse (mean) age grade rwage uhourse (sd) sdage=age sdeduc=grade sdrwage=rwage sdloghours=uhourse if sex==1 & uhourse>0 & uhourse!=., by(occ8090)
pwcorr age grade sdage sdeduc sdloghours sdrwage 

use org_ls, clear
gen grade=.
replace grade=gradeat if gradeat!=.
replace grade=grade92 if grade92!=.
collapse (mean) age grade rwage (sd) sdage=age sdeduc=grade sdrwage=rwage sdloghours=uhourse if sex==1 & uhourse>0 & uhourse!=., by(workgroup)
pwcorr age grade sdage sdeduc sdloghours sdrwage

log close


/*** PSID Analysis
1.  Means and Standard Deviations of all regression variables
2.  Percent of wage changes that entail a change of occupation
3.  Regressions
***/

capture log close
clear
set mem 500m
use psid_ls, clear
log using psid_analysis_ls, text replace

collapse (mean) loghours wage commincsd communemp worksd workunemp occsd occunemp
outsheet using desc_stats_psid.csv, comma replace

use psid_ls, clear
collapse (mean) loghours wage commincsd communemp worksd workunemp occsd occunemp, by(year)
outsheet using desc_stats_yearly_psid.csv, comma replace

use psid_ls, clear

sort panelid year

count if grwage05==1
count if grwage10==1
count if grwage15==1
count if grwage05==1 & occ8090!=occ8090[_n-1]
count if grwage10==1 & occ8090!=occ8090[_n-1]
count if grwage15==1 & occ8090!=occ8090[_n-1]


/*** Regressions
1.  All employed heads, no unemployment, community and workgroup dummies 
2.  All employed heads, no unemployment, community dummies, occupation inequality
3.  All employed heads, unemployment, community and workgroup dummies
4.  All fulltime heads, unemployment, community and workgroup dummies
5.  All salaried heads, unemployment, community and workgroup dummies
6.  All hourly heads, unemployment, community and workgroup dummies
7.  Marginal effects for all main effect plus interaction terms.
***/

foreach var of varlist loghours logwage logage occsd commincsd occsdxwage commincsdxwage worksd worksdxwage communemp workunemp occunemp {
	bys panelid: egen mean`var'=mean(`var')
	bys panelid:gen d`var'=`var'-mean`var'
}

local vars1 "logwage logage worksd commincsd worksdxwage commincsdxwage"
local vars2 "logwage logage occsd commincsd occsdxwage commincsdxwage"
local vars3 "logwage logage worksd commincsd worksdxwage commincsdxwage workunemp communemp"
local vars4 "logwage logage occsd commincsd occsdxwage commincsdxwage occunemp communemp"

/*** Logs ***/


egen occid=group(occ) if loghours!=. & logwage!=. & logage!=. & occsd!=. & commincsd!=. & communemp!=.
egen msaid=group(msa) if loghours!=. & logwage!=. & logage!=. & occsd!=. & commincsd!=. & communemp!=.
egen occxmsa=group(occid msaid)
egen workxmsa=group(workgroup msaid)

/*** 1 ***/

reg loghours logwage logage worksd commincsd worksdxwage commincsdxwage, cluster(msaid)
estout using regressions_ls.csv, style(tab) replace cells(b(star fmt(%9.6g)) se(par(`"="("' `")""') fmt(%9.6g))) stats(N, fmt(%9.0g) ) starlevels(* 0.10 ** 0.05 *** 0.01) varlabels(_cons Constant)  title("Pooled, nounemp, work") stard
parmest, saving(est1a, replace) rename(stderr stderr1) 

reg loghours logwage logage worksd commincsd worksdxwage commincsdxwage, cluster(workgroup)
parmest, saving(est1b, replace) rename(stderr stderr2)

reg loghours logwage logage worksd commincsd worksdxwage commincsdxwage, cluster(workxmsa)
parmest, saving(est1c, replace) rename(stderr stderr3)


/*** 2 ***/
reg loghours logwage logage occsd commincsd occsdxwage commincsdxwage, cluster(msaid) 
estout using regressions_ls.csv, style(tab) append cells(b(star fmt(%9.6g)) se(par(`"="("' `")""') fmt(%9.6g))) starlevels(* 0.10 ** 0.05 *** 0.01) stats(N, fmt(%9.0g)) varlabels(_cons Constant)  title("Pooled, nounemp, occ") stard
parmest, saving(est2a, replace) rename(stderr stderr1)

reg loghours logwage logage occsd commincsd occsdxwage commincsdxwage, cluster(occid)
parmest, saving(est2b, replace) rename(stderr stderr2)

reg loghours logwage logage occsd commincsd occsdxwage commincsdxwage, cluster(occxmsa)
parmest, saving(est2c, replace) rename(stderr stderr3)

/*** 3 ***/
areg loghours `vars1', cluster(msaid) absorb(panelid) 
estout using regressions_ls.csv, style(tab) append cells(b(star fmt(%9.6g)) se(par(`"="("' `")""') fmt(%9.6g))) starlevels(* 0.10 ** 0.05 *** 0.01) stats(N, fmt(%9.0g)) varlabels(_cons Constant)  title("Cluster, nounemployment") stard 
parmest, saving(est3a, replace) rename(stderr stderr1)

areg loghours `vars1', cluster(workgroup) absorb(panelid)
parmest, saving(est3b, replace) rename(stderr stderr2)

areg loghours `vars1', cluster(workxmsa) absorb(panelid)
parmest, saving(est3c, replace) rename(stderr stderr3)

/*** 4 ***/
areg loghours `vars2', cluster(msaid) absorb(panelid)  
estout using regressions_ls.csv, style(tab) append cells(b(star fmt(%9.6g)) se(par(`"="("' `")""') fmt(%9.6g))) starlevels(* 0.10 ** 0.05 *** 0.01) stats(N, fmt(%9.0g)) varlabels(_cons Constant)  title("Cluster, nounemployment, occ") stard
parmest, saving(est4a, replace) rename(stderr stderr1)

areg loghours `vars2', cluster(occid) absorb(panelid)
parmest, saving(est4b, replace) rename(stderr stderr2)

areg loghours `vars2', cluster(occxmsa) absorb(panelid)
parmest, saving(est4c, replace) rename(stderr stderr3)

/*** 5 ***/
areg loghours `vars3', cluster(msaid)  absorb(panelid)   
estout using regressions_ls.csv, style(tab) append cells(b(star fmt(%9.6g)) se(par(`"="("' `")""') fmt(%9.6g))) starlevels(* 0.10 ** 0.05 *** 0.01) stats(N, fmt(%9.0g)) varlabels(_cons Constant)  title("Base, unemployment") stard 
parmest, saving(est5a, replace) rename(stderr stderr1)

areg loghours `vars3', cluster(workgroup) absorb(panelid)
parmest, saving(est5b, replace) rename(stderr stderr2)

areg loghours `vars3', cluster(workxmsa) absorb(panelid)
parmest, saving(est5c, replace) rename(stderr stderr3)

/*** 6 ***/
areg loghours `vars4', cluster(msaid) absorb(panelid)  
estout using regressions_ls.csv, style(tab) append cells(b(star fmt(%9.6g)) se(par(`"="("' `")""') fmt(%9.6g))) starlevels(* 0.10 ** 0.05 *** 0.01) stats(N, fmt(%9.0g)) varlabels(_cons Constant)  title("Base, unemployment, occ") stard 
parmest, saving(est6a, replace) rename(stderr stderr1)

areg loghours `vars4', cluster(occid) absorb(panelid)
parmest, saving(est6b, replace) rename(stderr stderr2)

areg loghours `vars4', cluster(occxmsa) absorb(panelid)
parmest, saving(est6c, replace) rename(stderr stderr3)

/*** 7 ***/
areg loghours `vars3'  if hours>=30, cluster(msaid) absorb(panelid)  
estout using regressions_ls.csv, style(tab) append cells(b(star fmt(%9.6g)) se(par(`"="("' `")""') fmt(%9.6g))) starlevels(* 0.10 ** 0.05 *** 0.01) stats(N, fmt(%9.0g)) varlabels(_cons Constant)  title("Fulltime, unemployment") stard 
parmest, saving(est7a, replace) rename(stderr stderr1)

areg loghours `vars3'  if hours>=30, cluster(workgroup) absorb(panelid)
parmest, saving(est7b, replace) rename(stderr stderr2)

areg loghours `vars3'  if hours>=30, cluster(workxmsa) absorb(panelid)
parmest, saving(est7c, replace) rename(stderr stderr3)

/*** 8 ***/
areg loghours `vars4' if hours>=30, cluster(msaid) absorb(panelid)  
estout using regressions_ls.csv, style(tab) append cells(b(star fmt(%9.6g)) se(par(`"="("' `")""') fmt(%9.6g))) starlevels(* 0.10 ** 0.05 *** 0.01) stats(N, fmt(%9.0g)) varlabels(_cons Constant)  title("fulltime, unemployment, occ") stard 
parmest, saving(est8a, replace) rename(stderr stderr1)

areg loghours `vars4' if hours>=30, cluster(occid) absorb(panelid)
parmest, saving(est8b, replace) rename(stderr stderr2)

areg loghours `vars4' if hours>=30, cluster(occxmsa) absorb(panelid)
parmest, saving(est8c, replace) rename(stderr stderr3)

/*** 9 ***/
areg loghours `vars3' if slry_hrly==1, cluster(msaid) absorb(panelid)  
estout using regressions_ls.csv, style(tab) append cells(b(star fmt(%9.6g)) se(par(`"="("' `")""') fmt(%9.6g))) starlevels(* 0.10 ** 0.05 *** 0.01) stats(N, fmt(%9.0g)) varlabels(_cons Constant)  title("Salaried, unemployment") stard 
parmest, saving(est9a, replace) rename(stderr stderr1)

areg loghours `vars3' if slry_hrly==1, cluster(workgroup) absorb(panelid)
parmest, saving(est9b, replace) rename(stderr stderr2)

areg loghours `vars3' if slry_hrly==1, cluster(workxmsa) absorb(panelid)
parmest, saving(est9c, replace) rename(stderr stderr3)

/*** 10 ***/
areg loghours `vars4' if slry_hrly==1 , cluster(msaid) absorb(panelid)
estout using regressions_ls.csv, style(tab) append cells(b(star fmt(%9.6g)) se(par(`"="("' `")""') fmt(%9.6g))) starlevels(* 0.10 ** 0.05 *** 0.01) stats(N, fmt(%9.0g)) varlabels(_cons Constant)  title("Salaried, unemployment, occ") stard 
parmest, saving(est10a, replace) rename(stderr stderr1)

areg loghours `vars4' if slry_hrly==1, cluster(occid) absorb(panelid)
parmest, saving(est10b, replace) rename(stderr stderr2)

areg loghours `vars4' if slry_hrly==1, cluster(occxmsa) absorb(panelid)
parmest, saving(est10c, replace) rename(stderr stderr3)

/*** 11 ***/
areg loghours `vars3' if slry_hrly==3, cluster(msaid) absorb(panelid)
estout using regressions_ls.csv, style(tab) append cells(b(star fmt(%9.6g)) se(par(`"="("' `")""') fmt(%9.6g))) starlevels(* 0.10 ** 0.05 *** 0.01) stats(N, fmt(%9.0g)) varlabels(_cons Constant)  title("Hourly, unemployment") stard 
parmest, saving(est11a, replace) rename(stderr stderr1)

areg loghours `vars3' if slry_hrly==3, cluster(workgroup) absorb(panelid)
parmest, saving(est11b, replace) rename(stderr stderr2)

areg loghours `vars3' if slry_hrly==3, cluster(workxmsa) absorb(panelid)
parmest, saving(est11c, replace) rename(stderr stderr3)

/***12 ***/
areg loghours `vars4' if slry_hrly==3, cluster(msaid) absorb(panelid) 
estout using regressions_ls.csv, style(tab) append cells(b(star fmt(%9.6g)) se(par(`"="("' `")""') fmt(%9.6g))) starlevels(* 0.10 ** 0.05 *** 0.01) stats(N, fmt(%9.0g)) varlabels(_cons Constant)  title("Hourly, unemployment, occ") stard 
parmest, saving(est12a, replace) rename(stderr stderr1)

areg loghours `vars4' if slry_hrly==3, cluster(occid) absorb(panelid)
parmest, saving(est12b, replace) rename(stderr stderr2)

areg loghours `vars4' if slry_hrly==3, cluster(occxmsa) absorb(panelid)
parmest, saving(est12c, replace) rename(stderr stderr3)

/***13 ***/

areg loghours logwage logage commincsd commincsdxwage communemp, cluster(msaid)  absorb(panelid)   
estout using regressions_ls.csv, style(tab) append cells(b(star fmt(%9.6g)) se(par(`"="("' `")""') fmt(%9.6g))) starlevels(* 0.10 ** 0.05 *** 0.01) stats(N, fmt(%9.0g)) varlabels(_cons Constant)  title("Base, unemployment") stard 

/*** Make data set of the regression results from the above regressions ***/
global letter "a b c"

forval i=1/12 {
	foreach l of global letter {
		use est`i'`l', clear
		gen n=_n
		save, replace
	}
}

forval i=1/12 {
	use est`i'a, clear
	local lett "b c"
	foreach l of local lett {
		merge parm using est`i'`l', sort
		drop _merge
	}
	keep parm estimate dof stderr1 stderr2 stderr3 n
	gen regno=`i'
	save est`i', replace
}

use est1, clear

forval i=2/12 {
	append using est`i'
}

save estimates, replace

forval i=1/12 {
	foreach l of global letter {
		erase est`i'`l'.dta
	}
}

forval i=1/12 {
	erase est`i'.dta
}

use estimates, clear
gen stderr=abs(stderr1+stderr2-stderr3)
gen p=normalden(estimate/stderr)
save, replace

keep parm estimate n regno
gen type=1
save coef, replace

use estimates, clear
keep parm stderr p n regno
gen type=2
save stderr, replace

use coef, clear
append using stderr
gen star=""
replace star="***" if p<=.1
replace star="**" if p<=.05
replace star="*" if p<=0.01
sort regno n type
save estimates, replace

replace estimate=stderr if estimate==.
keep parm estimate stderr star regno type
outsheet using cluster_est.csv, comma replace

erase coef.dta
erase stderr.dta

log close