# Load required libraries
library(wooldridge)  
library(lmtest)      
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
library(car)         
## Loading required package: carData
# Load the data
data(wage2)

—- Part (i): Original Model —-

cat("\n=== Part (i): Original Model Analysis ===\n")
## 
## === Part (i): Original Model Analysis ===
model1 <- lm(log(wage) ~ educ + exper + tenure + married + 
             black + south + urban, data = wage2)

# Display results
summary_model1 <- summary(model1)
cat("\nOriginal Model Results:\n")
## 
## Original Model Results:
print(summary_model1)
## 
## Call:
## lm(formula = log(wage) ~ educ + exper + tenure + married + black + 
##     south + urban, data = wage2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.98069 -0.21996  0.00707  0.24288  1.22822 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  5.395497   0.113225  47.653  < 2e-16 ***
## educ         0.065431   0.006250  10.468  < 2e-16 ***
## exper        0.014043   0.003185   4.409 1.16e-05 ***
## tenure       0.011747   0.002453   4.789 1.95e-06 ***
## married      0.199417   0.039050   5.107 3.98e-07 ***
## black       -0.188350   0.037667  -5.000 6.84e-07 ***
## south       -0.090904   0.026249  -3.463 0.000558 ***
## urban        0.183912   0.026958   6.822 1.62e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3655 on 927 degrees of freedom
## Multiple R-squared:  0.2526, Adjusted R-squared:  0.2469 
## F-statistic: 44.75 on 7 and 927 DF,  p-value: < 2.2e-16
# Calculate the approximate percentage difference in wages for blacks
black_coef <- coef(model1)["black"]
black_effect <- (exp(black_coef) - 1) * 100

cat("\nWage Differential Analysis for Blacks vs Non-blacks:\n")
## 
## Wage Differential Analysis for Blacks vs Non-blacks:
cat("Coefficient:", round(black_coef, 4), "\n")
## Coefficient: -0.1883
cat("Standard Error:", round(summary_model1$coefficients["black", "Std. Error"], 4), "\n")
## Standard Error: 0.0377
cat("t-statistic:", round(summary_model1$coefficients["black", "t value"], 4), "\n")
## t-statistic: -5.0004
cat("p-value:", round(summary_model1$coefficients["black", "Pr(>|t|)"], 4), "\n")
## p-value: 0
cat("Percentage difference in wages:", round(black_effect, 2), "%\n")
## Percentage difference in wages: -17.17 %

—- Part (ii): Adding Squared Terms —-

cat("\n=== Part (ii): Analysis with Squared Terms ===\n")
## 
## === Part (ii): Analysis with Squared Terms ===
model2 <- lm(log(wage) ~ educ + exper + I(exper^2) + tenure + 
             I(tenure^2) + married + black + south + urban, 
             data = wage2)

# Test joint significance of exper^2 and tenure^2
squared_test <- linearHypothesis(model2, 
                                c("I(exper^2) = 0", "I(tenure^2) = 0"))

cat("\nJoint Test Results for Squared Terms:\n")
## 
## Joint Test Results for Squared Terms:
print(squared_test)
## 
## Linear hypothesis test:
## I(exper^2) = 0
## I(tenure^2) = 0
## 
## Model 1: restricted model
## Model 2: log(wage) ~ educ + exper + I(exper^2) + tenure + I(tenure^2) + 
##     married + black + south + urban
## 
##   Res.Df    RSS Df Sum of Sq      F Pr(>F)
## 1    927 123.82                           
## 2    925 123.42  2   0.39756 1.4898  0.226

—- Part (iii): Education Returns by Race —-

cat("\n=== Part (iii): Education Returns by Race ===\n")
## 
## === Part (iii): Education Returns by Race ===
model3 <- lm(log(wage) ~ educ + exper + tenure + married + 
             black + south + urban + educ:black, 
             data = wage2)

# Display results
summary_model3 <- summary(model3)
cat("\nModel with Race-Education Interaction:\n")
## 
## Model with Race-Education Interaction:
print(summary_model3)
## 
## Call:
## lm(formula = log(wage) ~ educ + exper + tenure + married + black + 
##     south + urban + educ:black, data = wage2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.97782 -0.21832  0.00475  0.24136  1.23226 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  5.374817   0.114703  46.859  < 2e-16 ***
## educ         0.067115   0.006428  10.442  < 2e-16 ***
## exper        0.013826   0.003191   4.333 1.63e-05 ***
## tenure       0.011787   0.002453   4.805 1.80e-06 ***
## married      0.198908   0.039047   5.094 4.25e-07 ***
## black        0.094809   0.255399   0.371 0.710561    
## south       -0.089450   0.026277  -3.404 0.000692 ***
## urban        0.183852   0.026955   6.821 1.63e-11 ***
## educ:black  -0.022624   0.020183  -1.121 0.262603    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3654 on 926 degrees of freedom
## Multiple R-squared:  0.2536, Adjusted R-squared:  0.2471 
## F-statistic: 39.32 on 8 and 926 DF,  p-value: < 2.2e-16
# Calculate return to education for blacks and non-blacks
educ_nonblack <- coef(model3)["educ"]
educ_interaction <- coef(model3)["educ:black"]
educ_black <- educ_nonblack + educ_interaction

cat("\nReturn to Education Analysis:\n")
## 
## Return to Education Analysis:
cat("For non-blacks:", round(educ_nonblack * 100, 2), "% per year of education\n")
## For non-blacks: 6.71 % per year of education
cat("For blacks:", round(educ_black * 100, 2), "% per year of education\n")
## For blacks: 4.45 % per year of education
cat("Difference:", round(educ_interaction * 100, 2), "% per year of education\n")
## Difference: -2.26 % per year of education
cat("P-value for difference:", round(summary_model3$coefficients["educ:black", "Pr(>|t|)"], 4), "\n")
## P-value for difference: 0.2626

—- Part (iv): Four Groups Analysis —-

cat("\n=== Part (iv): Four Groups Analysis ===\n")
## 
## === Part (iv): Four Groups Analysis ===
# Create interaction variable
wage2$married_black <- wage2$married * wage2$black

model4 <- lm(log(wage) ~ black + married + married_black + 
             educ + exper + tenure + south + urban, 
             data = wage2)

# Display results
summary_model4 <- summary(model4)
cat("\nFour Groups Model Results:\n")
## 
## Four Groups Model Results:
print(summary_model4)
## 
## Call:
## lm(formula = log(wage) ~ black + married + married_black + educ + 
##     exper + tenure + south + urban, data = wage2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.98013 -0.21780  0.01057  0.24219  1.22889 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    5.403793   0.114122  47.351  < 2e-16 ***
## black         -0.240820   0.096023  -2.508 0.012314 *  
## married        0.188915   0.042878   4.406 1.18e-05 ***
## married_black  0.061354   0.103275   0.594 0.552602    
## educ           0.065475   0.006253  10.471  < 2e-16 ***
## exper          0.014146   0.003191   4.433 1.04e-05 ***
## tenure         0.011663   0.002458   4.745 2.41e-06 ***
## south         -0.091989   0.026321  -3.495 0.000497 ***
## urban          0.184350   0.026978   6.833 1.50e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3656 on 926 degrees of freedom
## Multiple R-squared:  0.2528, Adjusted R-squared:  0.2464 
## F-statistic: 39.17 on 8 and 926 DF,  p-value: < 2.2e-16
# Calculate wage differential between married blacks and married nonblacks
married_black_coef <- coef(model4)["black"] + coef(model4)["married_black"]
married_black_diff <- (exp(married_black_coef) - 1) * 100

cat("\nWage Differential for Married Blacks vs Married Non-blacks:\n")
## 
## Wage Differential for Married Blacks vs Married Non-blacks:
cat("Coefficient sum:", round(married_black_coef, 4), "\n")
## Coefficient sum: -0.1795
cat("Percentage difference:", round(married_black_diff, 2), "%\n")
## Percentage difference: -16.43 %
# Save all results to a file (optional)
sink("wage_analysis_results.txt")
cat("=== Complete Wage Analysis Results ===\n")
cat("\n=== Part (i) Results ===\n")
print(summary_model1)
cat("\n=== Part (ii) Results ===\n")
print(squared_test)
cat("\n=== Part (iii) Results ===\n")
print(summary_model3)
cat("\n=== Part (iv) Results ===\n")
print(summary_model4)
sink()