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

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 + age + female + hsdipl\)

ln_ahe_1 <- lm(log(ahe) ~ yrseduc + age + female + hsdipl, 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.18318604  0.01655335  71.477 < 2.2e-16 ***
## yrseduc      0.08142350  0.00115179  70.693 < 2.2e-16 ***
## age          0.00748968  0.00024823  30.172 < 2.2e-16 ***
## female      -0.24409049  0.00485128 -50.315 < 2.2e-16 ***
## hsdipl       0.12699253  0.00978548  12.978 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

0.95 confidence interval for \(female\) (task 5.c)

Note that we take the coefficient estimated for the model with controls - so, our confidence interval (also heteroskedastic) takes them into account.

c_int_female <- confint(ln_ahe_1, coef = "female", vcov. = vcovHC(ln_ahe_1, type = "HC1"))

cat("\n95% Confidence Interval for the coefficient on'female':\n", c_int_female["female",])
## 
## 95% Confidence Interval for the coefficient on'female':
##  -0.2536607 -0.2345203

V-marginal effect (Task 5.d)

coef_educ <- coef(ln_ahe_1)["yrseduc"]
coef_hsdipl <- coef(ln_ahe_1)["hsdipl"]

# Calculate the estimated marginal value (V)
V_estimate <- coef_educ * 2 + coef_hsdipl
print(V_estimate)
##   yrseduc 
## 0.2898395

Task (5.e)

V confidence interval.

# Extract HC1 robust standard errors
se_coef_educ <- sqrt(diag(vcovHC(ln_ahe_1, type = "HC1")))["yrseduc"]
se_coef_hsdipl <- sqrt(diag(vcovHC(ln_ahe_1, type = "HC1")))["hsdipl"]

# Calculate the estimated SE
se_V_estimate <- sqrt((2 * se_coef_educ)^2 + se_coef_hsdipl^2)

# Calculate the critical value
alpha <- 0.05
df <- df.residual(ln_ahe_1)
t_critical <- qt(1 - alpha / 2, df)

# Calculate the margin of error
margin_of_error <- t_critical * se_V_estimate

# Calculate the lower and upper bounds
lower_bound <- V_estimate - margin_of_error
upper_bound <- V_estimate + margin_of_error

# Print
cat("95% Confidence Interval for V_estimate:", "[", lower_bound, ";", upper_bound, "]", "\n")
## 95% Confidence Interval for V_estimate: [ 0.2701354 ; 0.3095436 ]

Task 6

OLS: \(ln(ahe)\) ~ \(yrseduc + age + age^2\) (Task 6.a)

ln_ahe_2 <- lm(log(ahe) ~ yrseduc + age + I(age^2), 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)  3.7598e-01  4.4161e-02   8.5139 < 2.2e-16 ***
## yrseduc      8.7370e-02  9.8281e-04  88.8985 < 2.2e-16 ***
## age          4.4478e-02  2.1149e-03  21.0305 < 2.2e-16 ***
## I(age^2)    -4.3857e-04  2.5034e-05 -17.5187 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Test for linearity vs. quadratic relationship (task 6.b)

# Null hypothesis: The relation between ln(ahe) and age is linear
# Alternative hypothesis: The relation between ln(ahe) and age is quadratic

# F-test for comparing models
l_model <- lm(log(ahe) ~ yrseduc + age + female + hsdipl, data = ern_data)
q_model <- lm(log(ahe) ~ yrseduc + age + I(age^2) + female + hsdipl, data = ern_data)

# Perform F-test
anova_test <- anova(l_model, q_model)
print(anova_test)
## Analysis of Variance Table
## 
## Model 1: log(ahe) ~ yrseduc + age + female + hsdipl
## Model 2: log(ahe) ~ yrseduc + age + I(age^2) + female + hsdipl
##   Res.Df    RSS Df Sum of Sq      F    Pr(>F)    
## 1  37805 8314.8                                  
## 2  37804 8240.2  1    74.666 342.55 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Check significance at 5% level

# Define a function
PVal_Ind <- function(p_value, sign_level = 0.05){
    if (p_value < sign_level) {
      cat("\nReject the null hypothesis\n")}
    else {
      cat("\nFail to reject the null hypothesis\n")}
}
# Use it
PVal_Ind(anova_test$"Pr(>F)"[2], sign_level = 0.05)
## 
## Reject the null hypothesis

So, depedancy isn’t linear, it’s quadratic.

Test for quadratic relationship vs. cubic relationship (task 6.d)

# Null hypothesis: The relation between ln(ahe) and age is quadratic
# Alternative hypothesis: The relation between ln(ahe) and age is cubic

c_model <- lm(log(ahe) ~ yrseduc + age + I(age^2) + I(age^3) + female + hsdipl, data = ern_data)

# Perform F-test
anova_test_cubic <- anova(q_model, c_model)
print(anova_test_cubic)
## Analysis of Variance Table
## 
## Model 1: log(ahe) ~ yrseduc + age + I(age^2) + female + hsdipl
## Model 2: log(ahe) ~ yrseduc + age + I(age^2) + I(age^3) + female + hsdipl
##   Res.Df    RSS Df Sum of Sq      F   Pr(>F)    
## 1  37804 8240.2                                 
## 2  37803 8236.0  1    4.2276 19.405 1.06e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Check significance at 5% level
PVal_Ind(anova_test_cubic$"Pr(>F)"[2], sign_level = 0.05)
## 
## Reject the null hypothesis

Thus, the dependency is not quadratic, it’s cubic.

Test for linear dependency vs. polynomial of up to cubic degree (task 6.e)

# Null hypothesis: The relation between ln(ahe) and age is linear.
# Alternative hy[othesis: The relation between ln(ahe) and age is a polynomial of up to cubic degree.

# Fit the polynomial regression model up to cubic degree
pol_3_model <- lm(log(ahe) ~ poly(age, 3) + yrseduc + female + hsdipl, data = ern_data)

# Perform F-test
anova_test_p <- anova(l_model, pol_3_model)
print(anova_test_p)
## Analysis of Variance Table
## 
## Model 1: log(ahe) ~ yrseduc + age + female + hsdipl
## Model 2: log(ahe) ~ poly(age, 3) + yrseduc + female + hsdipl
##   Res.Df    RSS Df Sum of Sq      F    Pr(>F)    
## 1  37805 8314.8                                  
## 2  37803 8236.0  2    78.894 181.06 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Check significance at 5% level
PVal_Ind(anova_test_p$"Pr(>F)"[2], sign_level = 0.05)
## 
## Reject the null hypothesis

Additional test between linear and cubic

#Perform F-test
anova_test_add <- anova(l_model, c_model)
print(anova_test_add)
## Analysis of Variance Table
## 
## Model 1: log(ahe) ~ yrseduc + age + female + hsdipl
## Model 2: log(ahe) ~ yrseduc + age + I(age^2) + I(age^3) + female + hsdipl
##   Res.Df    RSS Df Sum of Sq      F    Pr(>F)    
## 1  37805 8314.8                                  
## 2  37803 8236.0  2    78.894 181.06 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Check significance at 5% level
PVal_Ind(anova_test_add$"Pr(>F)"[2], sign_level = 0.05)
## 
## Reject the null hypothesis

Task 7

Filling the table of marginal effect for different forms of the coefficient on “age”.

# Linear specification
print(coeftest(l_model, vcov = vcovHC, type = 'HC1'))
## 
## t test of coefficients:
## 
##                Estimate  Std. Error t value  Pr(>|t|)    
## (Intercept)  1.18318604  0.01655335  71.477 < 2.2e-16 ***
## yrseduc      0.08142350  0.00115179  70.693 < 2.2e-16 ***
## age          0.00748968  0.00024823  30.172 < 2.2e-16 ***
## female      -0.24409049  0.00485128 -50.315 < 2.2e-16 ***
## hsdipl       0.12699253  0.00978548  12.978 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Quadratic specification
print(coeftest(q_model, vcov = vcovHC, type = 'HC1'))
## 
## t test of coefficients:
## 
##                Estimate  Std. Error t value  Pr(>|t|)    
## (Intercept)  4.4075e-01  4.3181e-02  10.207 < 2.2e-16 ***
## yrseduc      8.1024e-02  1.1468e-03  70.651 < 2.2e-16 ***
## age          4.4953e-02  2.0555e-03  21.870 < 2.2e-16 ***
## I(age^2)    -4.4099e-04  2.4303e-05 -18.145 < 2.2e-16 ***
## female      -2.4473e-01  4.8329e-03 -50.639 < 2.2e-16 ***
## hsdipl       1.2240e-01  9.7704e-03  12.528 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Cubic specification
print(coeftest(c_model, vcov = vcovHC, type = 'HC1'))
## 
## t test of coefficients:
## 
##                Estimate  Std. Error  t value  Pr(>|t|)    
## (Intercept) -3.0429e-01  1.7412e-01  -1.7476   0.08054 .  
## yrseduc      8.1149e-02  1.1475e-03  70.7208 < 2.2e-16 ***
## age          1.0072e-01  1.2865e-02   7.8290 5.043e-15 ***
## I(age^2)    -1.7805e-03  3.0811e-04  -5.7787 7.589e-09 ***
## I(age^3)     1.0331e-05  2.3852e-06   4.3314 1.485e-05 ***
## female      -2.4423e-01  4.8345e-03 -50.5190 < 2.2e-16 ***
## hsdipl       1.2256e-01  9.7681e-03  12.5472 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Estimating the marginal values

For the linear specification:

Predicted increase = \(\beta_1\)

For the quadratic specification:

Predicted increase = \(\beta_1\) + \(2 \beta _2 \cdot \text{Age}\)

For the cubic specification:

Predicted increase = \(\beta_1\) + \(2 \beta _2 \cdot \text{Age}\) + \(3 \beta _3 \cdot \text{Age}^2\)

But we have discrete variable age! Better to use just differences between predicted values (for age = b)

# Linear
l_age <- coeftest(l_model, vcov = vcovHC, type = 'HC1')["age", "Estimate"]

# Quadratic
q_age <- coeftest(q_model, vcov = vcovHC, type = 'HC1')["age", "Estimate"]
q_age_2 <- coeftest(q_model, vcov = vcovHC, type = 'HC1')["I(age^2)", "Estimate"]

# Cubic
c_age <- coeftest(c_model, vcov = vcovHC, type = 'HC1')["age", "Estimate"]
c_age_2 <- coeftest(c_model, vcov = vcovHC, type = 'HC1')["I(age^2)", "Estimate"]
c_age_3 <- coeftest(c_model, vcov = vcovHC, type = 'HC1')["I(age^3)", "Estimate"]

# Calculate predicted values for age 31
age31_linear <- l_age
age31_quadratic <- q_age + q_age_2 * (31^2 - 30^2)
age31_cubic <- c_age +c_age_2 * (31^2 - 30^2) + c_age_3 * (31^3 - 30^3)

# Calculate predicted values for age 46
age46_linear <- l_age
age46_quadratic <- q_age + q_age_2 * (46^2 - 45^2)
age46_cubic <- c_age +c_age_2 * (46^2 - 45^2) + c_age_3 * (46^3 - 45^3)

# Calculate predicted values for age 61
age61_linear <- l_age
age61_quadratic <- q_age + q_age_2 * (61^2 - 60^2)
age61_cubic <- c_age +c_age_2 * (61^2 - 60^2) + c_age_3 * (61^3 - 60^3)

results <- data.frame(
  Specification = c("Linear in age", "Quadratic in age", "Cubic in age"),
  `Age 30 to 31` = c(age31_linear, age31_quadratic, age31_cubic),
  `Age 45 to 46` = c(age46_linear, age46_quadratic, age46_cubic),
  `Age 60 to 61` = c(age61_linear, age61_quadratic, age61_cubic)
)

# Print the dataframe
print(results)
##      Specification Age.30.to.31 Age.45.to.46 Age.60.to.61
## 1    Linear in age  0.007489682  0.007489682  0.007489682
## 2 Quadratic in age  0.018052510  0.004822834 -0.008406843
## 3     Cubic in age  0.020950625  0.002869873 -0.001263508