use happy eduyrs using ESS1-6e01_0_F1.dta, clear
// Pull in outliers
qui sum eduyrs, detail
replace eduyrs = r(p99) if eduyrs >= r(p99) & !missing(eduyrs)
// Calculate density functions at different years of education
kdensity happy if eduyrs == 4, n(5000) k(gauss) bw(.8) gen(h0 d0) nograph
kdensity happy if eduyrs == 8, n(5000) k(gauss) bw(.8) gen(h1 d1) nograph
kdensity happy if eduyrs == 12, n(5000) k(gauss) bw(.8) gen(h2 d2) nograph
kdensity happy if eduyrs == 16, n(5000) k(gauss) bw(.8) gen(h3 d3) nograph
// Scale and place density functions
replace d0 = (d0*18) + 4
replace d1 = (d1*18) + 8
replace d2 = (d2*18) + 12
replace d3 = (d3*18) + 16
// Plot
twoway (line h0 d0 if inrange(h0,3,10)) ///
(line h1 d1 if inrange(h1,3,10)) ///
(line h2 d2 if inrange(h2,3,10)) ///
(line h3 d3 if inrange(h3,3,10)) ///
(lfit happy eduyrs) ///
, legend(off) xline(4 8 12 16) xlabel(0 (1) 22) ///
xtitle("Years of education") ytitle("Happiness") ///
yscale(range(3 10)) ylabel(3(1)10)
