cd "C:\Users\Reimers\Desktop\Pre-semester Econometrics\Pre-semester Econometrics\2019"

clear all

set obs 10000
set seed 1 

gen s = runiform(0,1)
	replace s = 1 if s <= 0.1
	replace s = 0 if s > 0.1 & s != 1
	
********************************************************************************
***			 	Example 1: Linear model w/ one expl. variable  				   * 
********************************************************************************
gen x1 = (_n/_N +2)^2 - 4 + rnormal(0,1)*0.5				//generate explanatory variable
gen x1sq = x1^2
gen c = 5.0
gen u = rnormal(0,1)						//generate error term with mean 0, variance 1
gen y1 = c + 2*x1 + u						//generate dependent variable
gen y2 = c + 2*x1 - 0.5*x1^2 + u			//generate dependent variable

local options1 scheme(s2mono) graphregion(color(white) fcolor(white)) legend(region(lcolor(white))) 
graph twoway (scatter y1 x1 if s == 1, msym(x) mc(gs3) ytick(0(5)25)),  xtitle(x) `options1'
		graph export "input/data_ex1.pdf" , as(pdf) replace
reg y1 x1		
		
graph twoway (scatter y2 x1 if s == 1, msym(x) mc(gs3) ytick(0(1)10)), xtitle(x) `options1'
		graph export "input/data_ex1_2.pdf" , as(pdf) replace

		
		
********************************************************************************
***			 		Example 2: Coefficient Interpretation 					   * 
********************************************************************************
gen D = runiform(0,1) 			//An Dummy variable
	replace D = 0 if D < 0.5
	replace D = 1 if D >= 0.5
	
gen lnX3 = ln(rnormal(0,1)*exp(_n/_N))	

*	1) y is continuous
	gen y_c = c -0.5*x1 -0.5*D + 0.25*lnX3 + u
		 reg y_c x1 D lnX3 if s == 1
		 eststo Cont , title("Cont.")
		
*	2) y is log form
	gen lny = 3*c -2.14*x1 + 0.275*D + 0.51*lnX3 + u
		 reg lny x1 D lnX3 if s == 1
		 eststo Log, title("Log")
		
*	3) y is dummy
	gen t = -5 + 1.5*x1 +2.25*D - 0.6*lnX3 + u
		gen y_D = 1 if t >= 0
		replace y_D = 0 if t < 0
		  reg y_c x1 D lnX3 if s == 1
		  eststo Dummy, title("Dummy")
		
	
	local optest1 /*title("")*/ mtitle("y is Continuous" "y is in log-terms" "y is a Dummy" ) b(a1) se(a2) nonumber coeflab(x1 "Dep.: Continuous" D "Dep.:  Dummy" lnX3 "Dep.:  Log")
	esttab Cont  Log Dummy using "input/CoeffInt.tex", replace  `optest1'   //fragm

	
	
********************************************************************************
***			 		Example 3: Homo- & Heteroskedasticity 					   * 
********************************************************************************
gen u_het = rnormal(0,1)*0.75*x1			//generate heteroskedastic error term with mean 0, variance 1

gen y_het1 = c + 2*x1  + u_het			//generate dependent variable
gen y_het2 = c + 2*x1 - 0.5*x1^2 + u_het			//generate dependent variable

*Check y to see homo/heteroskedasticity
graph twoway (scatter y1 x1 if s == 1, msym(x) mc(gs3) ylab(0(5)30)) ///
			 (lfit y1 x1 if s == 1,  lc(cranberry) ylab(0(5)30)), ytitle(y) ylab(0(5)30) xtitle(x) `options1'
		graph export "input/data_ex1_lfit.pdf" , as(pdf) replace		
graph twoway (scatter y_het1 x1 if s == 1, msym(x) mc(gs3) ylab(0(5)30)) ///
		     (lfit y_het1 x1 if s == 1,  lc(cranberry) ylab(0(5)30) ), ytitle(y) ylab(0(5)30) xtitle(x) `options1'
		graph export "input/data_het.pdf" , as(pdf) replace

*Check the error term to see homo/heteroskedasticity
graph twoway (scatter u x1 if s == 1, msym(x) mc(gs3) ylab(-10(5)20)) ///
		     (lfit u x1 if s == 1,  lc(cranberry) ylab(-10(5)20)), ytitle(u (homoskedastic)) xtitle(x) `options1'
		graph export "input/data_u.pdf" , as(pdf) replace
graph twoway (scatter u_het x1 if s == 1, msym(x) mc(gs3) ) ///
			 (lfit u_het x1 if s == 1,  lc(cranberry) ), ytitle(u (heteroskedastic)) xtitle(x) ylab(-10(5)20) `options1'
		graph export "input/data_uhet.pdf" , as(pdf) replace		
		

*Compare standard errors in homo/hetero/het-robust regression setups		
*1) Homoskedasticity
reg y1 x1	if s == 1
estat hettest 
eststo hom, title("hom")

*2) Heteroskedasticity
reg y_het1 x1 	if s == 1
estat hettest 
eststo het, title("het")

*3) Heteroskedasticity robust 
reg y_het1 x1 if s == 1, rob	
eststo hetrob, title("hetrob")

*--> Adj. for heterosked. in standard errors_ Var(hatbeta) increases

local optest1 /*title("")*/ mtitle("Homoskedastic u" "Heteroskedastic u" "Het. u, robust Std. Errors" ) b(a1) se(a2) nonumber coeflab(x1 "$ \beta_1 $" _cons "$ \beta_0 $")
esttab hom  het hetrob using "input/HomHet.tex", replace  `optest1'   //fragm


reg y2 x1 x1sq	if s == 1
estat hettest 

reg y_het2 x1 	if s == 1
estat hettest 

reg y_het2 x1 x1sq if s == 1, rob	
	

********************************************************************************
***			 		Example 4: Monte Carlo Simulation	 					   * 
********************************************************************************

foreach s in 02 05 075 {
	gen mcb_`s' = .

	qui {
		forvalues n = 1/1000 {
			gen s_ = runiform(0,1)
				replace s_= 1 if s_ <= 0.`s'				//Sample size: ~200 Obs., 2% of pop.
				replace s_ = 0 if s_ > 0.`s' & s_ != 1
			
				reg y1 x1	if s_ == 1
				matrix A = e(b)
				local b = A[1,1]
					replace mcb_`s' = `b' in `n'
			
			drop s_
		}
	}
}

local options1  graphregion(color(white) fcolor(white)) legend(region(lcolor(white))) 

graph twoway	(kdensity mcb_02, lc(cranberry) legend(label(1 "2% Sample" )) ) /// 
				(kdensity mcb_05, lc(navy) legend(label(2 "5% Sample" )) ) ///
				(kdensity mcb_075, lc(black) legend(label(3 "7.5% Sample" )) ), `options1' ytitle("Density (in %)") xtitle("Sample estimates")  ///
				legend(order(1 2 3)) ytick(0(1)15) xtick(1.8(0.05)2.2)  //lab(1 "2% Sample" 2 "5% Sample" 3 "7.5% Sample") )
	
		graph export "input/mc_sim.pdf" , as(pdf) replace		
		
				
********************************************************************************
***			 		Example 5: Coefficient Interpretation 					   * 
********************************************************************************				
clear all
use "USA_all_std.dta" //--> CPS 2005 data, downloadable at online site of US Bureau of Labor Statistics 				
				
keep hid pid hr_wage age hhsize employed hrs male emp_wj edu ind wgt				
order pid hr_wage age hhsize employed hrs male emp_wj edu ind wgt	 			
	
bysort male: sum employed [aweight = wgt]
tab edu if employed == 1
sum hr_wage age hhsize male if employed == 1 [aweight = wgt] 			
				
graph twoway	(kdensity hr_wage if employed == 1 [aweight = wgt], lc(cranberry) )			//--> What a right skewed distribution!!		
	
gen ln_wage = log(hr_wage)
sum ln_wage if employed == 1 [aweight = wgt]
graph twoway (kdensity ln_wage if employed == 1 [aweight = wgt], lc(cranberry) )				//--> log-normalizing helps in reduced skewness	

*-----------------------------------------------------------------*	
*Example 1: depvar = wage (cont.), indepvar = age, age2, male, edu	
gen age2 = age^2		
				
reg hr_wage i.edu age age2 male if employed == 1 [aweight = wgt]
eststo contcont, title("Cont.-Cont.")

*1.1) Interprete!
*		edu: baseline is less than secondary education
*		age: Hump shape, display() command for critical age
*		male: Compare men/women of same age/educ - gap in wages!
*1.2) Explanatory power R2?
*1.3) Example for perfect multicollinearity
gen exper = age-16 if edu == 1				//less than college
	replace exper = age-18 if edu == 2		//secondary completed
	replace exper = age-25 if edu == 3		//more than secondary

reg hr_wage i.edu age age2 exper male if employed == 1 [aweight = wgt] //Surprise!

*Example 2: depvar = log(wage) (cont.), indepvar = age, age2, male, edu	
reg ln_wage i.edu age age2 male if employed == 1 [aweight = wgt]
eststo logcont, title("Log-Cont.")

*2.1)	Interprete coefficients (semi-elasticities)
*2.2)	Relate semi-elasticity to estimate in cont.-cont. model
*2.3)	Mention: Log-log regressions
*2.4)	IMPORTANT: Estimated model vs. implied relationship in levels!!!!!		

*Example 3: LPM. depvar = employment (cont.), indepvar = age, age2, male, edu	
reg employed i.edu age age2 male [aweight = wgt]
eststo LPM, title("LPM")

*3.1) 	Interprete
*3.2)	Show how LPM over/undershoots above 1/below0
predict hat_emp, xb
sum hat_emp, detail

sort hat_emp
gen id = _n
graph twoway scatter hat_emp  id  


*Example 4: Partial effect w/ interaction term
reg ln_wage i.edu age age2 male male#i.edu if employed == 1 [aweight = wgt]
eststo Int, title("Int.Action")

*4.1) Interprete. Given setup, best to interprete for being a men on top of being a women

local coeflab 
local optest1 /*title("")*/ mtitle("Cont.-Cont." "Log-Cont." "Linear Prob. Model" "incl. Interactions" ) b(a1) se(a2) nonumber coeflab(" `coeflab' ")
esttab contcont logcont LPM using "input/CoeffInt_new.tex", replace  `optest1'   //fragm


********************************************************************************
***			 			Example 6: 	Hypothesis Testing	 					   * 
********************************************************************************
gen edu1 = 0 
	replace edu1 = 1 if edu == 1
gen edu2 = 0 
	replace edu2 = 1 if edu == 2
gen edu3 = 0 
	replace edu3 = 1 if edu == 3	
	
reg ln_wage edu2 edu3 age age2 male  if employed == 1 [aweight = wgt]

*5.1) Test for equal effects of secondary/tertiary education
matrix A = e(b)
test  edu2 = edu3 = A[1,1]
testparm edu2 edu3

*------------------------------------------------------------------------------*
*				Hypothesis testing under collinearity						   *

reg ln_wage i.edu age age2 exper male if employed == 1 [aweight = wgt]

replace exper = exper + rnormal(0,1)		//Add a little bit of randomness
reg ln_wage i.edu exper  male if employed == 1 [aweight = wgt] 			//Experience alone
reg ln_wage i.edu exper age age2  male if employed == 1 [aweight = wgt] //incl. age --> Multicollinearity


********************************************************************************
***			 			Example 7: Heteroskedasticity	 					   * 
********************************************************************************
reg ln_wage i.edu age age2 male  if employed == 1 [aweight = wgt]
estat hettest 
//We have heteroskedasticity, clearly. Maybe including other variables can help?
//F.e.: Ability - unobserved, pos. correlated with educ and earnings. As educ rises,
//		ability might also be rising. It's in the error, so that u_i rises as x increases
//		--> Volatility not constant!

reg ln_wage i.edu age age2 male hhsize i.ind if employed == 1 [aweight = wgt]
estat hettest 
predict hat_wage, xb
gen hat_u = ln_wage-hat_wage				

*Check graphically for heteroskedasticity
*1) res & age
	graph twoway (scatter hat_u age if employed == 1 [aweight = wgt], ylab(-10(2)10) mc(gs0) msym(x) msize(tiny))			//--> What a right skewed distribution!!		
*	--> Rather(more or less constant

*2)	res & educ
	graph twoway (scatter hat_u edu if employed == 1 [aweight = wgt], xlab(1(1)3) ylab(-15(3)15) mc(gs0) msym(x) msize(tiny))			//--> What a right skewed distribution!!		
*	--> There some heteroskedasticity here!

*3) res & hhsize?
	graph twoway (scatter hat_u hhsize if employed == 1 [aweight = wgt], xlab(1(1)10) ylab(-15(3)15) mc(gs0) msym(x) msize(tiny))			//--> What a right skewed distribution!!		
*	--> Absolutely clear!
				
*Manual version of the heteroskedasticity test from White
*1) Run main reg, extract fitted residuals 
*	--> see above

*2) Run regression of squared fitted res on expl. vars/yhat
	gen hat_u2 = hat_u^2
	reg hat_u2 i.edu age age2 male hhsize i.ind if employed == 1 [aweight = wgt]
*	--> Check F stat!

*	Check a more encompassing model w/ interaction terms...
	reg hat_u2 i.edu age age2 male hhsize i.ind age#i.ind i.edu#male i.edu#hhsize if employed == 1 [aweight = wgt]
*	--> Check F stat!

*	Check a model that just relies on fitted dep. variable
	reg hat_u2 hat_wage if employed == 1 [aweight = wgt]
*	--> Check F stat!

*Baseline:  We have to adjust error terms for heteroskedasticity. To do this, in Stata you can add 
*			the "vce(robust)" option to your reg command
reg ln_wage i.edu age age2 male hhsize i.ind if employed == 1 [aweight = wgt], vce(robust)
			
*------------------------------------------------------------------------------*
*							Example: feasible GLS							   *
*1)	Run main reg, extract residuals
local regv edu2 edu3 age age2 male hhsize

reg ln_wage `regv' if employed == 1 [aweight = wgt]
predict fit if employed == 1, xb
predict res if employed == 1, residuals
*	--> See above

*2)	Compute transformed residuals, reg on x
	gen res2 = res^2			//Conditional variance of u given x
	gen ln_res2 = log(res2)		
	reg ln_res2 `regv' if employed == 1 [aweight = wgt]
	
*3) Compute predicted residuals hat g, transform
	predict hat_g if employed == 1, xb			//Fitted value for log(conditional variance of u given x)
	gen hat_h = exp(hat_g)		//retransform to levels

*4	Estimate transformed model
	local tregv
	foreach v in `regv' {
		gen t`v' = `v'/sqrt(hat_h)
		local tregv `tregv' t`v'	
	}
reg ln_wage `tregv' if employed == 1 [aweight = wgt]

reg ln_wage edu2 edu3 age age2 male hhsize if employed == 1 [aweight = wgt]







