Libraries

# install.packages('lmtest')

library(readr) # read csv
library(lmtest) # stat tests
library(sandwich) # cov matrixes
library(readxl) # csv and excel files
library(psych) # descriptive statistics
library(knitr) # for tables
library(xtable) # for latex tables
library(ggplot2) # for graphs
library(car) # for linear hypothesis testing
library(stargazer) # latex regression tables

Downloading the data

ern_data <- read_csv(
  'C:/Users/Popov/Documents/Studies/NES_studies/R/Econometrics_1/HA3/cps99_ps3.csv',
                      show_col_types = FALSE)

Descriptive statistics

Task 5

OLS: \(ln(ahe)\) ~ \(yrseduc + female + yrseduc*female\)

ln_ahe_1 <- lm(log(ahe) ~ yrseduc + female + yrseduc*female, data = ern_data)
rob_se_ln_ahe_1<- sqrt(diag(vcovHC(ln_ahe_1, type = "HC1"))) #for plots
print(coeftest(ln_ahe_1, vcov = vcovHC, type = 'HC1'))
## 
## t test of coefficients:
## 
##                  Estimate Std. Error  t value  Pr(>|t|)    
## (Intercept)     1.5875698  0.0174548  90.9533 < 2.2e-16 ***
## yrseduc         0.0826756  0.0012774  64.7199 < 2.2e-16 ***
## female         -0.4635077  0.0270460 -17.1377 < 2.2e-16 ***
## yrseduc:female  0.0166647  0.0019651   8.4803 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Regressions for men/women (Task 5.b)

# Define intercepts and slopes for men and women
intercept_men <- coef(ln_ahe_1)["(Intercept)"]
slope_men <- coef(ln_ahe_1)["yrseduc"]

intercept_women <- coef(ln_ahe_1)["(Intercept)"] + coef(ln_ahe_1)["female"]
slope_women <- coef(ln_ahe_1)["yrseduc"] + coef(ln_ahe_1)["yrseduc:female"]

# Create data for plotting
plot_data <- data.frame(yrseduc = seq(min(ern_data$yrseduc), max(ern_data$yrseduc), length.out = 100))

# Calculate predicted values for men and women
plot_data$pred_men <- intercept_men + slope_men * plot_data$yrseduc
plot_data$pred_women <- intercept_women + slope_women * plot_data$yrseduc

# Create data for legend
legend_data <- data.frame(
  group = c("Men", "Women"),
  intercept = c(intercept_men, intercept_women),
  slope = c(slope_men, slope_women)
)

# Plot both regression lines and add intercepts and slopes to the legend
ggplot(plot_data, aes(x = yrseduc)) +
  geom_line(aes(y = pred_men, color = "Men"), linetype = "solid") +
  geom_line(aes(y = pred_women, color = "Women"), linetype = "solid") +
  labs(x = "Years of Education", y = "log(Average Hourly Earnings)") +
  ggtitle("Regression Lines for Men and Women") +
  scale_color_manual(values = c("blue", "red"), labels = c("Men", "Women")) +
  theme_minimal() +
  theme(legend.title = element_blank(),
        legend.position = "top",
        legend.box.margin = margin(0, 0, 0, 0),
        legend.spacing.x = unit(0.1, "cm")) +
  scale_color_manual(values = c("blue", "red"), 
                     labels = c(paste("Men Intercept:",
                                      round(intercept_men,3),
                                      ", Slope:",round(slope_men, 3)),
                                paste("Women Intercept:",
                                      round(intercept_women, 3),
                                       ", Slope:", round(slope_women, 3))))

Alternative way

ggplot(data = ern_data) +
  geom_point(aes(x = yrseduc, y = log(ahe), color = factor(female)),
             size = 0.7) +
  scale_color_manual(values = c("#0000A5", "#D30085")) +
  geom_abline(
    aes(
      intercept = ln_ahe_1$coefficients["(Intercept)"], 
      slope = ln_ahe_1$coefficients["yrseduc"]), 
      color = "#007AFF", 
      size = 1.1) +
  geom_abline(
    aes(intercept = (
        ln_ahe_1$coefficients["(Intercept)"] + 
        ln_ahe_1$coefficients["female"]), 
      slope = (
        ln_ahe_1$coefficients["yrseduc"] + 
        ln_ahe_1$coefficients["yrseduc:female"])), 
      color = "#FF00A2", 
      size = 1.1) +
  labs(
    color = "Female",
    x = "Years of Education",
    y = "Average Hourly Earnings in Logs") +
  theme_light()

Test for the difference between regression lines for men/women (task 5.d)

f_H0 <- c("female = 0", #same intercepts
          "yrseduc:female = 0") # same slopes
linearHypothesis(ln_ahe_1, f_H0, vcov = vcovHC(ln_ahe_1, type = "HC1"))
## Linear hypothesis test
## 
## Hypothesis:
## female = 0
## yrseduc:female = 0
## 
## Model 1: restricted model
## Model 2: log(ahe) ~ yrseduc + female + yrseduc * female
## 
## Note: Coefficient covariance matrix supplied.
## 
##   Res.Df Df      F    Pr(>F)    
## 1  37808                        
## 2  37806  2 1213.2 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

So, H0 is rejected - regression lines are not the same for men/women.

Task 6

OLS: \(ln(ahe)\) ~ \(yrseduc + female + female:yrseduc + male:yrseduc\)

# Create the variable 'male'
ern_data$male <- ifelse(ern_data$female == 0, 1, 0)

# Estimate the regression
ln_ahe_2 <- lm(log(ahe) ~ yrseduc + female + female:yrseduc + male:yrseduc, data = ern_data)
rob_se_ln_ahe_2<- sqrt(diag(vcovHC(ln_ahe_2, type = "HC1"))) #for plots
print(coeftest(ln_ahe_2, vcov = vcovHC, type = 'HC1'))
## 
## t test of coefficients:
## 
##                  Estimate Std. Error  t value  Pr(>|t|)    
## (Intercept)     1.5875698  0.0174548  90.9533 < 2.2e-16 ***
## yrseduc         0.0826756  0.0012774  64.7199 < 2.2e-16 ***
## female         -0.4635077  0.0270460 -17.1377 < 2.2e-16 ***
## yrseduc:female  0.0166647  0.0019651   8.4803 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

male:yrseduc was dropped for causing multicollinearity

Task 7

OLS \(log(ahe)\) ~ yrseduc + female + female:yrseduc + age + I(age^2) + I(age^3)

ln_ahe_i <- lm(log(ahe) ~ yrseduc + female + female:yrseduc + age + I(age^2) + I(age^3), data = ern_data)
se_ln_ahe_i<- sqrt(diag(vcovHC(ln_ahe_i, type = "HC1"))) #for Stargayzer 


ln_ahe_ii <- lm(log(ahe) ~ yrseduc + female + female:yrseduc + age + I(age^2) + I(age^3) + female:age + female:I(age^2) + female:I(age^3), data = ern_data)
se_ln_ahe_ii<- sqrt(diag(vcovHC(ln_ahe_ii, type = "HC1"))) #for Stargayzer 


stargazer(ln_ahe_i, ln_ahe_ii, se = list(se_ln_ahe_i, se_ln_ahe_ii), type = "text", keep.stat = c("n"))
## 
## ===========================================
##                    Dependent variable:     
##                ----------------------------
##                          log(ahe)          
##                     (1)            (2)     
## -------------------------------------------
## yrseduc           0.081***      0.080***   
##                   (0.001)        (0.001)   
##                                            
## female           -0.519***        0.042    
##                   (0.027)        (0.350)   
##                                            
## age               0.103***      0.113***   
##                   (0.013)        (0.017)   
##                                            
## I(age2)          -0.002***      -0.002***  
##                   (0.0003)      (0.0004)   
##                                            
## I(age3)          0.00001***    0.00001***  
##                  (0.00000)      (0.00000)  
##                                            
## yrseduc:female    0.021***      0.020***   
##                   (0.002)        (0.002)   
##                                            
## female:age                       -0.025    
##                                  (0.026)   
##                                            
## female:I(age2)                   0.0003    
##                                  (0.001)   
##                                            
## female:I(age3)                  -0.00000   
##                                 (0.00000)  
##                                            
## Constant           -0.230        -0.452*   
##                   (0.174)        (0.233)   
##                                            
## -------------------------------------------
## Observations       37,810        37,810    
## ===========================================
## Note:           *p<0.1; **p<0.05; ***p<0.01

Gender gap (years = 25, educ = 16) (Task 7.b)

# Define the values for the variables
yrseduc_value <- 16
age_value <- 25

# Calculate the predicted earnings for males and females using regression (ii)
predicted_male <- predict(ln_ahe_ii, newdata = data.frame(yrseduc = yrseduc_value, female = 0, # for male 
age = age_value, `I(age^2)` = age_value^2, `I(age^3)` = age_value^3))


predicted_female <- predict(ln_ahe_ii, newdata = data.frame(yrseduc = yrseduc_value, female = 1, # for female
age = age_value, `I(age^2)` = age_value^2, `I(age^3)` = age_value^3))

# Calculate the gender gap
gender_gap <- predicted_male - predicted_female

# Print the gender gap
cat("Gender gap:\n", gender_gap)
## Gender gap:
##  0.08428115

Gender gap (years = 55, educ = 16) (Task 7.c)

# Define the values for the variables
age_value_2 <- 55

# Calculate the predicted earnings for males and females using regression (ii)
predicted_male_2 <- predict(ln_ahe_ii, newdata = data.frame(yrseduc = yrseduc_value, female = 0, # for male 
age = age_value_2, `I(age^2)` = age_value_2^2, `I(age^3)` = age_value_2^3))


predicted_female_2 <- predict(ln_ahe_ii, newdata = data.frame(yrseduc = yrseduc_value, female = 1, # for female
age = age_value_2, `I(age^2)` = age_value_2^2, `I(age^3)` = age_value_2^3))

# Calculate the gender gap
gender_gap_2 <- predicted_male_2 - predicted_female_2

# Print the gender gap
cat("Gender gap:\n", gender_gap_2)
## Gender gap:
##  0.2269328

Does gender gap depend on age? (Task 7.d)

f_H0 <- c("female:age = 0", 
          "female:I(age^2) = 0",
          "female:I(age^3) = 0") # if so, gender gap doesn't depend on age
linearHypothesis(ln_ahe_ii, f_H0, vcov = vcovHC(ln_ahe_ii, type = "HC1"))
## Linear hypothesis test
## 
## Hypothesis:
## female:age = 0
## female:I(age^2) = 0
## female:I(age^3) = 0
## 
## Model 1: restricted model
## Model 2: log(ahe) ~ yrseduc + female + female:yrseduc + age + I(age^2) + 
##     I(age^3) + female:age + female:I(age^2) + female:I(age^3)
## 
## Note: Coefficient covariance matrix supplied.
## 
##   Res.Df Df      F    Pr(>F)    
## 1  37803                        
## 2  37800  3 27.624 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The null hypothesis is declined, gender gap does depend on age.