# /* Examples of Regression Models
# The probability of labor participation using the Current Population Data
# Data  cps_msa2013 : 2013 Current population survey for 3 MSA
# 3 MSA : NY, LA, CHICAGO
# Created by Jin Man Lee
# Updated: 11/14/2022 */  
# Clear all memory and working space
rm(list=ls())

# Load all required libraries. If not installed, use "install.packages("foreign") 
library(foreign)
library(dplyr)
library(stargazer)
library(lmtest)
library(whitestrap)

# Use the following line if you use R in bigblue server
# cps <- read.dta("/var/www/html/jlee141/econdata/cps_data/cps_msa2013.dta")
cps <- read.dta("https://bigblue.depaul.edu/jlee141/econdata/cps_data/cps_msa2013.dta") 

# Create some variables and remove any invalid cases 
cps1 <- cps %>%  
  mutate(educ92a = as.numeric(educ92)) %>% 
  mutate(educyr = case_when( 
    educ92a %in%  1 ~ 1,
    educ92a %in%  2 ~ 4,
    educ92a %in%  3 ~ 6,
    educ92a %in%  4 ~ 8,
    educ92a %in%  5 ~ 9,
    educ92a %in%  6 ~ 10,
    educ92a %in%  7 ~ 11,
    educ92a %in%  8:9 ~ 12,
    educ92a %in%  10:12 ~ 14,
    educ92a %in%  13 ~ 16,
    educ92a %in%  14 ~ 18,
    educ92a %in%  15 ~ 20,
    educ92a %in%  16 ~ 21)) %>% 
  mutate(edusq = educyr^2,
         lhrwage = log(hrwage)) %>% 
  mutate(laborp = ifelse(as.numeric(lfstat) == 1 | as.numeric(lfstat) ==2,1, 
                         ifelse(is.na(lfstat),NA,0))) %>% 
  filter(hrwage > 0) %>% 
  dplyr::select(female,educyr,lhrwage,edusq,wbho) %>% 
  filter(complete.cases(.)) 

# Let's make sure the data has any issues
summary(cps1)
table(cps1$wbho)
table(cps1$wbho,cps1$female)

# Regression with dummy variables
ols1 <- lm(lhrwage ~ educyr + wbho,data=cps1)
summary(ols1)
cps1$yhat1 <- predict(ols1)
plot(x=cps1$educyr,y=cps1$yhat1,col=cps1$wbho)

ols2 <- lm(lhrwage ~ educyr + female,data=cps1)
summary(ols2)
cps1$yhat2 <- predict(ols2)
plot(x=cps1$educyr,y=cps1$yhat2)

ols3 <- lm(lhrwage ~ educyr + female + female*educyr,data=cps1)
summary(ols3)
cps1$yhat3 <- predict(ols3)
plot(x=cps1$educyr,y=cps1$yhat3)

# a Breusch-Pagan Test to determine if heteroscedasticity is present
bptest(ols3)

# White's test
white_test(ols3)

# Summary Table for the Regression Model
stargazer(ols1,ols2,ols3,type="text",title = "Comparison the OLS regressions")