Assignment 1

1. Implement the data-generating process from the previous slide in R, Stata or Python.

N = 1000 # sample size
x = rnorm (N,0,1)
mean(x); sd(x)
## [1] 0.01878364
## [1] 1.023271
e = rnorm(N,0,1)
data = data.frame(x=x,e=e)
data$y = 2 + 1*data$x + data$e

summary(lm(y~x, data = data))
## 
## Call:
## lm(formula = y ~ x, data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.2951 -0.6759 -0.0048  0.6855  2.7192 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.00594    0.03217   62.35   <2e-16 ***
## x            1.03060    0.03145   32.77   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.017 on 998 degrees of freedom
## Multiple R-squared:  0.5183, Adjusted R-squared:  0.5178 
## F-statistic:  1074 on 1 and 998 DF,  p-value: < 2.2e-16

2. Define coefficients that differ from the ones used, and recover them in a linear regression.

## I changed the beta0 to 5 and beta1 to 3
data$y1 = 5 + 3*data$x + data$e
summary(lm(y1~x, data = data))
## 
## Call:
## lm(formula = y1 ~ x, data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.2951 -0.6759 -0.0048  0.6855  2.7192 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  5.00594    0.03217  155.60   <2e-16 ***
## x            3.03060    0.03145   96.36   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.017 on 998 degrees of freedom
## Multiple R-squared:  0.9029, Adjusted R-squared:  0.9029 
## F-statistic:  9285 on 1 and 998 DF,  p-value: < 2.2e-16

Answer:
Based on 1,000 observations, the OLS estimates of the intercept (beta0) and slope (beta1) are 5.00663 and 2.96591, respectively. Both estimates are close to their true values of 5 and 3.

3. Vary the values of your coefficients; vary the value of your R^2.

# replacement coefficient
# 4 combinations of beta0 and beta1
coef_settings = data.frame(
  group = c("small positive slope",
            "original coef",
            "large positive slope",
            "negative slope"),
  beta0 = c(2,2,3,5),
  beta1 = c(0.5,1,3,-2)
)
coef_results = data.frame()

for (i in 1:nrow(coef_settings)) {
  # generate y using the same x and error term
  y = coef_settings$beta0[i] + coef_settings$beta1[i]*x + e
  model = lm(y~x)
  model_summary = summary(model)
  
  # store the results
  coef_results = rbind(
    coef_results,
    data.frame(
      group = coef_settings$group[i],
      true_beta0 = coef_settings$beta0[i],
      estimated_beta0 = coef(model)[1],
      true_beta1 = coef_settings$beta1[i],
      estimated_beta1 = coef(model)[2],
      standard_error = model_summary$coefficients["x","Std. Error"],
      R_squared = model_summary$r.squared
    )
  )
}
coef_results[, 2:7] = round(coef_results[, 2:7], 3)
coef_results

Answer:
The estimated coefficients for each group are close to the true values. Since the sample size, x, and noise are identical, the standard errors of the slopes are essentially the same. The larger the value of beta1, the more variation is explained by x, and the R^2 is typically higher, changing beta0 itself does not alter the R^2.

## keep the original coefficients fixed
beta0 = 2
beta1 = 1

target_R2 = c(0.2, 0.5, 0.9)
error_sd <- abs(beta1) * sqrt((1 - target_R2) / target_R2)

r2_results <- data.frame()
r2_samples <- list()

for (i in seq_along(target_R2)) {
  
  # Generate an error term with the required standard deviation
  e_i <- rnorm(N, mean = 0, sd = error_sd[i])
  
  # Generate y
  y_i <- beta0 + beta1 * x + e_i
  
  sample_data <- data.frame(x = x, y = y_i)
  model <- lm(y ~ x, data = sample_data)
  model_summary <- summary(model)
  
  # Save the data and model for plotting later
  r2_samples[[i]] <- list(
    data = sample_data,
    model = model
  )
  
  # Store the regression results
  r2_results <- rbind(
    r2_results,
    data.frame(
      target_R2 = target_R2[i],
      error_sd = error_sd[i],
      estimated_beta1 = coef(model)[2],
      standard_error =
        model_summary$coefficients["x", "Std. Error"],
      actual_R2 = model_summary$r.squared
    )
  )
}

round(r2_results, 3)

Answer:
The lower the error, the closer the data points are to the regression line, the higher the R^2, and the more precise the slope estimate.

4. Plot the observed and predicted values for y against the observed values for x.

par(mfrow = c(1, 3))
for (i in seq_along(target_R2)) {
  
  sample_data = r2_samples[[i]]$data
  model = r2_samples[[i]]$model
  plot(
    sample_data$x,
    sample_data$y,
    col = rgb(0, 0, 1, 0.3),
    pch = 16,
    xlab = "Observed x",
    ylab = "Observed and predicted y",
    main = paste("Target R² =", target_R2[i])
  )
  # Add the fitted regression line
  abline(model, col = "orange", lwd = 3)
}

# Restore the normal plotting layout
par(mfrow = c(1, 1))

Answer:
As R-squared increases, the observed values become closer to the fitted line. This shows that a higher R-squared indicates less noise and a better model fit.

5. Put your simulation in a loop, store the results from each iteration, run at least 100 iterations, and interpret

set.seed(123)

number_of_simulations = 100
simulation_results = data.frame()

for (i in seq_along(target_R2)) {
  for (simulation in 1:number_of_simulations) {
    x_sim = rnorm(N, mean = 0, sd = 1)
    e_sim = rnorm(N, mean = 0, sd = error_sd[i])
    y_sim = beta0 + beta1 * x_sim + e_sim
    
    model = lm(y_sim ~ x_sim)
    model_summary = summary(model)
    
    simulation_results = rbind(
      simulation_results,
      data.frame(
        target_R2 = target_R2[i],
        simulation = simulation,
        estimated_beta1 = coef(model)[2],
        standard_error =
          model_summary$coefficients["x_sim", "Std. Error"],
        actual_R2 = model_summary$r.squared
      )
    )
  }
}

# Calculate the average result for each target R-squared
average_results = aggregate(
  cbind(estimated_beta1, standard_error, actual_R2) ~ target_R2,
  data = simulation_results,
  FUN = mean
)
round(average_results, 3)

Answer:
After 100 simulations, the average beta1 estimates remain close to the true value of 1, and the average R-squared(R^2) values match their targets. Higher R-squared values are associated with smaller standard errors and more precise estimates.

Assignment 2 - seeing the endogeneity bias

1. Simulate one regressor x (standard normal) correlated with the normal error (ρ = 0.3); define y = 2 + 1· x; regress y on x.

N = 1000 # sample size
rho = 0.3
x = rnorm (N,0,1)
e_pre = rnorm(N,0,1)
e = rho * x + sqrt(1-rho^2)*e_pre # to keep the sd of e is 1 and cor(x,e)=0.3
y = 2 + 1*x + e
data = data.frame(x=x,e=e,y=y)

cor(data$x, data$e)
## [1] 0.3067492
model = lm(y~x, data = data)
summary(model)
## 
## Call:
## lm(formula = y ~ x, data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.0957 -0.6124  0.0073  0.6089  3.4686 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.99942    0.03069   65.16   <2e-16 ***
## x            1.30580    0.03004   43.48   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9695 on 998 degrees of freedom
## Multiple R-squared:  0.6544, Adjusted R-squared:  0.6541 
## F-statistic:  1890 on 1 and 998 DF,  p-value: < 2.2e-16

Answer:
The cor(x,e) = 0.3315744, and the estimated slope coefficient beta1 is 1.34054, compared with 1.00608 in assignment1. In this task, we set the population correlation between x and e to 0.3. It means as x increase, e tends to increase as well. As a result, OLS attributes part of the variation in y arising from e to x, leading to an upward-biased estimate of the beta1. Although the coefficient is statistically significant (p < 0.05), statistical significance does not imply that the estimate is unbiased.

2. Check the correlation between your regressor and the estimated error term (residuals). Is this a meaningful analysis?

data$residual = residuals(model)
true_cor = cor(data$x, data$e) # x and true error 
residual_cor = cor(data$x, data$residual) # x and residuals

correlation_results = data.frame(
  comparison = c(
    "x and true error", "x and OLS residual"
  ),
  correlation = c(
    true_cor, residual_cor
  )
)
correlation_results

Answer:
Although x is correlated with the true error term, its correlation with the OLS residuals is approximately zero. This is not a meaningful test for endogeneity because OLS residuals are uncorrelated with the included regressor by construction.

3. Repeat your simulation. Do you get the same results?

### 
set.seed(123)
N = 1000
rho = 0.3
number_of_simulations = 100 # number of simulations
simulation_results = data.frame() # to store the results

for (i in 1:number_of_simulations){
  ## generate different random x and e in every simulation
  x_sim = rnorm(N,0,1)
  e_pre_sim = rnorm(N,0,1)
  e_sim = rho * x_sim + sqrt(1 - rho^2) * e_pre_sim # error term
  y_sim = 2 + 1*x_sim + e_sim # dependent variable
  
  model_sim = lm(y_sim ~ x_sim)
  
  ### store the result
  simulation_results = rbind(
    simulation_results,
    data.frame(
      simulation = i,
      estimated_beta1 = coef(model_sim)[2],
      cor_x_true_error = cor(x_sim, e_sim),
      cor_x_residual = cor(x_sim, residuals(model_sim))
    )
  )
}
simulation_results
## calculate average
average_results = data.frame(
  results = c(
    "Avg. Estimated beta1", "Avg. x & true error", "Avg x and OLS residual"
  ),
  average = c(
    mean(simulation_results$estimated_beta1),
    mean(simulation_results$cor_x_true_error),
    mean(simulation_results$cor_x_residual)
  )
)
average_results$average = round(average_results$average,3)
average_results

Answer:
The average estimated slope is approximately 1.3 rather than the true value of 1, showing systematic endogeneity bias. The correlation between x and the OLS residuals remains approximately 0 in every simulation, so it cannot be used to detect endogeneity.

Assignment 3 - the data-rich approach

Create two correlated variables x1 and x2 (ρ = 0.3). Let the error be uncorrelated with x1 and x2. Define y = 2 + 1· x1 + 1· x2 + e. Run two regressions: excluding x2, and including x2. Compare.

library(MASS)
set.seed(123)

N = 1000
sigma = matrix(c(1,0.3,0.3,1),2,2) # correlation matrix
rhs = mvrnorm(N,mu=c(0,0),Sigma = sigma)
x1 = rhs[,1]; x2 = rhs[,2]
e = rnorm(N)
df = data.frame(x1,x2,y=2 + 1 * x1 + 1 * x2 + e)

summary(lm(y~x1, data = df)) # x2 omitted
## 
## Call:
## lm(formula = y ~ x1, data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.5006 -0.8972 -0.0682  0.9499  4.4154 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.02135    0.04438   45.55   <2e-16 ***
## x1           1.27570    0.04645   27.46   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.403 on 998 degrees of freedom
## Multiple R-squared:  0.4304, Adjusted R-squared:  0.4299 
## F-statistic: 754.2 on 1 and 998 DF,  p-value: < 2.2e-16
summary(lm(y~x1+x2, data = df))
## 
## Call:
## lm(formula = y ~ x1 + x2, data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.8360 -0.6277 -0.0370  0.6538  3.3787 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.97907    0.03098   63.88   <2e-16 ***
## x1           0.96342    0.03380   28.51   <2e-16 ***
## x2           1.00992    0.03110   32.47   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9788 on 997 degrees of freedom
## Multiple R-squared:  0.7232, Adjusted R-squared:  0.7226 
## F-statistic:  1302 on 2 and 997 DF,  p-value: < 2.2e-16

Answer:
When x2 is omitted, the estimated coefficient of x1 is 1.276, which is larger than its true value of 1 because x1 is positively correlated with x2. After including x2, the estimated coefficients of x1 and x2 are 0.963 and 1.010, respectively, both close to their true values. The R-squared also increases from 0.430 to 0.723. This shows that including the relevant variable x2 removes the omitted-variable bias and improves model fit.

Assignment 4 - three routes to the same fixed-effects estimate

Create a panel data set with, e.g., 100 cross-sectional units i, each with its own intercept, each observed T = 100 times. Create a regressor x correlated with both the intercepts and the idiosyncratic error.

set.seed(123)
n_id = 100
n_time = 100
N = n_id * n_time

id = rep(1:n_id, each = n_time)
time = rep(1:n_time, times = n_id)

alpha_i = rnorm(n_id)
alpha = alpha_i[id]

e = rnorm(N)
v = rnorm(N) # provides additional independent random variation

x = 0.5 * alpha + 0.5 * e + v  # makes x correlated with the unit-specific intercept, idiosyncratic error
y = alpha + 1 * x + e

panel = data.frame(id = id, time = time, alpha = alpha, x = x, e = e, y = y)

round(cor(panel[, c("x", "alpha", "e")]), 3)
##           x  alpha      e
## x     1.000  0.374  0.409
## alpha 0.374  1.000 -0.017
## e     0.409 -0.017  1.000

3. Estimate a fixed-effects model three ways — dummy variables, demeaning, integrated estimator — and compare the results in one table.

dummy variables
model_dummy = lm(y~x+factor(id), data = panel)
summary(model_dummy)
## 
## Call:
## lm(formula = y ~ x + factor(id), data = panel)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.2748 -0.5984  0.0006  0.6001  3.3520 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   -0.519089   0.089470  -5.802 6.76e-09 ***
## x              1.398688   0.008010 174.622  < 2e-16 ***
## factor(id)2    0.348738   0.126529   2.756 0.005858 ** 
## factor(id)3    1.729195   0.126798  13.637  < 2e-16 ***
## factor(id)4    0.656857   0.126515   5.192 2.12e-07 ***
## factor(id)5    0.532376   0.126538   4.207 2.61e-05 ***
## factor(id)6    1.803607   0.126753  14.229  < 2e-16 ***
## factor(id)7    0.976768   0.126565   7.718 1.30e-14 ***
## factor(id)8   -0.455175   0.126463  -3.599 0.000321 ***
## factor(id)9   -0.009441   0.126460  -0.075 0.940487    
## factor(id)10   0.253061   0.126472   2.001 0.045427 *  
## factor(id)11   1.446952   0.126748  11.416  < 2e-16 ***
## factor(id)12   0.763777   0.126568   6.035 1.65e-09 ***
## factor(id)13   0.766896   0.126552   6.060 1.41e-09 ***
## factor(id)14   0.683017   0.126516   5.399 6.87e-08 ***
## factor(id)15   0.084581   0.126465   0.669 0.503633    
## factor(id)16   2.054840   0.126840  16.200  < 2e-16 ***
## factor(id)17   0.971950   0.126540   7.681 1.73e-14 ***
## factor(id)18  -0.985085   0.126557  -7.784 7.75e-15 ***
## factor(id)19   1.083575   0.126664   8.555  < 2e-16 ***
## factor(id)20   0.076295   0.126459   0.603 0.546313    
## factor(id)21  -0.436289   0.126474  -3.450 0.000564 ***
## factor(id)22   0.268411   0.126468   2.122 0.033833 *  
## factor(id)23  -0.182394   0.126476  -1.442 0.149298    
## factor(id)24  -0.047393   0.126463  -0.375 0.707849    
## factor(id)25   0.112724   0.126461   0.891 0.372751    
## factor(id)26  -0.927582   0.126527  -7.331 2.46e-13 ***
## factor(id)27   1.166870   0.126626   9.215  < 2e-16 ***
## factor(id)28   0.749018   0.126482   5.922 3.29e-09 ***
## factor(id)29  -0.280104   0.126473  -2.215 0.026801 *  
## factor(id)30   1.398198   0.126810  11.026  < 2e-16 ***
## factor(id)31   0.691072   0.126558   5.461 4.86e-08 ***
## factor(id)32   0.331388   0.126469   2.620 0.008799 ** 
## factor(id)33   1.237657   0.126685   9.770  < 2e-16 ***
## factor(id)34   1.250014   0.126604   9.873  < 2e-16 ***
## factor(id)35   1.224604   0.126636   9.670  < 2e-16 ***
## factor(id)36   1.058637   0.126584   8.363  < 2e-16 ***
## factor(id)37   0.783281   0.126516   6.191 6.21e-10 ***
## factor(id)38   0.496834   0.126512   3.927 8.65e-05 ***
## factor(id)39   0.391136   0.126482   3.092 0.001991 ** 
## factor(id)40   0.092624   0.126479   0.732 0.463987    
## factor(id)41  -0.053027   0.126459  -0.419 0.674990    
## factor(id)42   0.349482   0.126474   2.763 0.005733 ** 
## factor(id)43  -0.558352   0.126469  -4.415 1.02e-05 ***
## factor(id)44   2.211159   0.126977  17.414  < 2e-16 ***
## factor(id)45   1.514616   0.126744  11.950  < 2e-16 ***
## factor(id)46  -0.290385   0.126476  -2.296 0.021699 *  
## factor(id)47   0.154174   0.126462   1.219 0.222821    
## factor(id)48   0.104031   0.126462   0.823 0.410739    
## factor(id)49   1.032240   0.126629   8.152 4.02e-16 ***
## factor(id)50   0.532131   0.126482   4.207 2.61e-05 ***
## factor(id)51   0.679451   0.126582   5.368 8.16e-08 ***
## factor(id)52   0.401858   0.126499   3.177 0.001494 ** 
## factor(id)53   0.439407   0.126495   3.474 0.000516 ***
## factor(id)54   1.668095   0.126809  13.154  < 2e-16 ***
## factor(id)55   0.399537   0.126483   3.159 0.001589 ** 
## factor(id)56   1.762394   0.126713  13.909  < 2e-16 ***
## factor(id)57  -0.661373   0.126526  -5.227 1.76e-07 ***
## factor(id)58   1.009548   0.126584   7.975 1.69e-15 ***
## factor(id)59   0.756599   0.126511   5.980 2.30e-09 ***
## factor(id)60   0.739126   0.126520   5.842 5.32e-09 ***
## factor(id)61   0.625949   0.126513   4.948 7.63e-07 ***
## factor(id)62   0.234103   0.126459   1.851 0.064169 .  
## factor(id)63   0.169531   0.126459   1.341 0.180082    
## factor(id)64  -0.211659   0.126460  -1.674 0.094217 .  
## factor(id)65  -0.234236   0.126459  -1.852 0.064017 .  
## factor(id)66   0.713424   0.126503   5.640 1.75e-08 ***
## factor(id)67   0.796454   0.126498   6.296 3.18e-10 ***
## factor(id)68   0.609302   0.126499   4.817 1.48e-06 ***
## factor(id)69   1.208065   0.126576   9.544  < 2e-16 ***
## factor(id)70   2.235855   0.126911  17.617  < 2e-16 ***
## factor(id)71   0.061422   0.126466   0.486 0.627202    
## factor(id)72  -1.185001   0.126596  -9.360  < 2e-16 ***
## factor(id)73   1.370811   0.126600  10.828  < 2e-16 ***
## factor(id)74   0.044906   0.126460   0.355 0.722524    
## factor(id)75  -0.156873   0.126459  -1.241 0.214819    
## factor(id)76   1.408373   0.126690  11.117  < 2e-16 ***
## factor(id)77   0.218361   0.126484   1.726 0.084308 .  
## factor(id)78  -0.387338   0.126478  -3.062 0.002201 ** 
## factor(id)79   0.530359   0.126491   4.193 2.78e-05 ***
## factor(id)80   0.393694   0.126472   3.113 0.001858 ** 
## factor(id)81   0.506754   0.126500   4.006 6.22e-05 ***
## factor(id)82   0.965200   0.126534   7.628 2.60e-14 ***
## factor(id)83   0.346864   0.126477   2.743 0.006108 ** 
## factor(id)84   1.073987   0.126562   8.486  < 2e-16 ***
## factor(id)85   0.414488   0.126480   3.277 0.001052 ** 
## factor(id)86   0.764736   0.126514   6.045 1.55e-09 ***
## factor(id)87   1.504853   0.126588  11.888  < 2e-16 ***
## factor(id)88   0.878129   0.126637   6.934 4.34e-12 ***
## factor(id)89   0.087210   0.126461   0.690 0.490449    
## factor(id)90   1.409419   0.126733  11.121  < 2e-16 ***
## factor(id)91   1.264387   0.126638   9.984  < 2e-16 ***
## factor(id)92   0.971690   0.126514   7.680 1.74e-14 ***
## factor(id)93   0.618836   0.126520   4.891 1.02e-06 ***
## factor(id)94   0.070489   0.126460   0.557 0.577268    
## factor(id)95   1.585326   0.126690  12.513  < 2e-16 ***
## factor(id)96   0.061977   0.126459   0.490 0.624079    
## factor(id)97   2.200579   0.126965  17.332  < 2e-16 ***
## factor(id)98   1.647033   0.126744  12.995  < 2e-16 ***
## factor(id)99   0.385264   0.126459   3.047 0.002321 ** 
## factor(id)100 -0.298829   0.126463  -2.363 0.018148 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8942 on 9899 degrees of freedom
## Multiple R-squared:  0.8436, Adjusted R-squared:  0.842 
## F-statistic:   534 on 100 and 9899 DF,  p-value: < 2.2e-16
demeaning
library(dplyr)
panel_grouped = group_by(panel, id)
panel_grouped = mutate(panel_grouped, 
                       y_dm = y - mean(y),
                       x_dm = x - mean(x))
panel = ungroup(panel_grouped)
model_demean = lm(y_dm ~ x_dm - 1, data = panel)
summary(model_demean)
## 
## Call:
## lm(formula = y_dm ~ x_dm - 1, data = panel)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.2748 -0.5984  0.0006  0.6001  3.3520 
## 
## Coefficients:
##      Estimate Std. Error t value Pr(>|t|)    
## x_dm  1.39869    0.00797   175.5   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8897 on 9999 degrees of freedom
## Multiple R-squared:  0.7549, Adjusted R-squared:  0.7549 
## F-statistic: 3.08e+04 on 1 and 9999 DF,  p-value: < 2.2e-16
integrated estimator
library(fixest)
model_fe = feols(y~x, data = panel, fixef = "id")
summary(model_fe)
## OLS estimation, Dep. Var.: y
## Observations: 10,000
## Fixed-effects: id: 100
## Standard-errors: IID 
##   Estimate Std. Error t value  Pr(>|t|)    
## x  1.39869    0.00801 174.622 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## RMSE: 0.889672     Adj. R2: 0.842043
##                  Within R2: 0.754926
comparison
comparison_results <- data.frame(
  method = c(
    "Dummy variables",
    "Manual demeaning",
    "Integrated FE estimator"
  ),
  true_beta1 = 1,

  estimated_beta1 = c(
    coef(model_dummy)["x"],
    coef(model_demean)["x_dm"],
    coef(model_fe)["x"]
  ),
  
  standard_error = c(
    coef(summary(model_dummy))["x", "Std. Error"],
    coef(summary(model_demean))["x_dm", "Std. Error"],
    sqrt(diag(vcov(model_fe)))["x"]
  )
)
round(comparison_results[, 2:4], 3)

Answer:
All three fixed-effects methods produce the same estimate of 1.399. However, it is higher than the true value of 1 because fixed effects remove the individual intercepts but do not eliminate the correlation between x and the idiosyncratic error.

Assignment 5 - the IV estimator by hand

Create a data set in which the error is correlated with the focal regressor, and an instrument is correlated with the regressor but not with the error. Estimate the model with OLS and check the magnitude of the OLS inconsistency.

set.seed(123)
N = 1000
z = rnorm(N,0,1)
e = rnorm(N)
x = 2*z+e+rnorm(N)
y = x + e + rnorm(N)

lm(y~x) # OLS biased.
## 
## Call:
## lm(formula = y ~ x)
## 
## Coefficients:
## (Intercept)            x  
##     0.02231      1.20130
x_hat = lm(x~z)$fitted.values
lm(y ~ x_hat) # manual IV
## 
## Call:
## lm(formula = y ~ x_hat)
## 
## Coefficients:
## (Intercept)        x_hat  
##     0.03106      1.04111
ols_inconsistency = coef(lm(y~x))["x"] - 1
ols_inconsistency
##         x 
## 0.2013028

Answer:
The OLS estimate is higher than the true coefficient of 1, showing an upward inconsistency caused by endogeneity. The magnitude of the OLS inconsistency is measured by beta_hat_ols - 1 (0.2013028).

Estimate the model using the manual three-step approach.

# Step 1: Regress the endogenous regressor x on the instrument z
first_stage = lm(x ~ z)
summary(first_stage)
## 
## Call:
## lm(formula = x ~ z)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.6703 -1.0018  0.0007  0.9513  4.9025 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.02124    0.04502   0.472    0.637    
## z            2.06898    0.04541  45.558   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.423 on 998 degrees of freedom
## Multiple R-squared:  0.6753, Adjusted R-squared:  0.675 
## F-statistic:  2075 on 1 and 998 DF,  p-value: < 2.2e-16
# Step 2: Obtain the predicted values of x
x_hat = first_stage$fitted.values

# Step 3: Replace x with x_hat in the main equation
second_stage = lm(y ~ x_hat)
summary(second_stage)
## 
## Call:
## lm(formula = y ~ x_hat)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -8.0731 -1.7326  0.0644  1.6113  8.2250 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.03106    0.07862   0.395    0.693    
## x_hat        1.04111    0.03833  27.165   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.485 on 998 degrees of freedom
## Multiple R-squared:  0.4251, Adjusted R-squared:  0.4245 
## F-statistic: 737.9 on 1 and 998 DF,  p-value: < 2.2e-16
comparison = c(true_beta1 = 1, OLS = unname(coef(lm(y ~ x))["x"]),
               manual_IV = unname(coef(second_stage)["x_hat"]))
round(comparison, 3)
## true_beta1        OLS  manual_IV 
##      1.000      1.201      1.041

Answer:
This shows that a valid instrument can correct the OLS inconsistency.

Assignment 6 - control function vs. 2SLS

1. Generate a data set plagued by endogeneity.

#### common data for assignment 6,7 and 8
set.seed(123)
N = 1000
z = rnorm(N)
e = rnorm(N)

x = 2 * z + e + rnorm(N) 
y = x + e + rnorm(N) 

data678 = data.frame(y,x,z,e)
cor(x, e)   # x is endogenous
## [1] 0.4834632
cor(z, x)   # z is relevant
## [1] 0.8217593
cor(z, e)   # z is exogenous
## [1] 0.08647944

2. Estimate the model with IV (2SLS) and with a control function; compare estimates and standard errors. Compute the correct control-function standard errors, either analytically or by bootstrapping.

# 2SLS
library(AER)
model_iv = ivreg(y ~ x, instruments = ~ z, data = data678)
summary(model_iv)
## 
## Call:
## ivreg(formula = y ~ x | z, data = data678)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.4124 -0.8858  0.0114  0.9464  4.0646 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.03106    0.04360   0.712    0.476    
## x            1.04111    0.02125  48.990   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.378 on 998 degrees of freedom
## Multiple R-Squared: 0.8232,  Adjusted R-squared: 0.8231 
## Wald test:  2400 on 1 and 998 DF,  p-value: < 2.2e-16
# control function
first_stage = lm(x~z,data=data678)
data678$res_first = residuals(first_stage) # residuals
model_cf = lm(y~x + res_first, data = data678)
summary(model_cf)
## 
## Call:
## lm(formula = y ~ x + res_first, data = data678)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.6354 -0.8102 -0.0098  0.8135  3.2819 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.03106    0.03753   0.828    0.408    
## x            1.04111    0.01829  56.908   <2e-16 ***
## res_first    0.49335    0.03211  15.367   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.186 on 997 degrees of freedom
## Multiple R-squared:  0.8691, Adjusted R-squared:  0.8689 
## F-statistic:  3311 on 2 and 997 DF,  p-value: < 2.2e-16
# comparison
comparison = rbind(IV_2SLS = c(estimate = unname(coef(model_iv)["x"]),
                               standard_error = unname(coef(summary(model_iv))["x", "Std. Error"])),
                   Control_function = c(estimate = unname(coef(model_cf)["x"]),
                                        standard_error = unname(coef(summary(model_cf))["x", "Std. Error"])))
round(comparison, 3)
##                  estimate standard_error
## IV_2SLS             1.041          0.021
## Control_function    1.041          0.018

Answer:
The 2SLS and control-function approaches produce the same estimated coefficient of 1.041, which is close to the true value of 1. This is expected because the two approaches are equivalent in this linear model. However, their reported standard errors differ: the 2SLS standard error is 0.021, while the naive control-function standard error is 0.018. The control-function standard error is not yet correct because it treats the estimated first-stage residual as an observed variable and ignores first-stage estimation uncertainty.

Assignment 7 - fit under IV vs. OLS

Generate a data set plagued by endogeneity. Estimate the model with IV and with OLS. Compute and compare the MSE and the R2 for both models.

library(AER)
# OSL
model_ols7 = lm(y ~ x, data = data678)
# IV
model_iv7 = ivreg(y~x, instruments = ~z, data = data678)

summary(model_ols7)
## 
## Call:
## lm(formula = y ~ x, data = data678)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.3909 -0.8790  0.0661  0.8972  4.0068 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.02231    0.04171   0.535    0.593    
## x            1.20130    0.01671  71.886   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.319 on 998 degrees of freedom
## Multiple R-squared:  0.8381, Adjusted R-squared:  0.838 
## F-statistic:  5168 on 1 and 998 DF,  p-value: < 2.2e-16
summary(model_iv7)
## 
## Call:
## ivreg(formula = y ~ x | z, data = data678)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.4124 -0.8858  0.0114  0.9464  4.0646 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.03106    0.04360   0.712    0.476    
## x            1.04111    0.02125  48.990   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.378 on 998 degrees of freedom
## Multiple R-Squared: 0.8232,  Adjusted R-squared: 0.8231 
## Wald test:  2400 on 1 and 998 DF,  p-value: < 2.2e-16
# Predicted values
y_hat_ols = predict(model_ols7)
y_hat_iv = predict(model_iv7)

# MSE
mse_ols = mean((data678$y - y_hat_ols)^2)
mse_iv = mean((data678$y - y_hat_iv)^2)

# R-squared
total_variation = sum((data678$y - mean(data678$y))^2)
r2_ols = 1 - sum((data678$y - y_hat_ols)^2) / total_variation
r2_iv = 1 - sum((data678$y - y_hat_iv)^2) / total_variation

fit_comparison = rbind(OLS = c(MSE = mse_ols, R2 = r2_ols),
                       IV = c(MSE = mse_iv, R2 = r2_iv))
round(fit_comparison, 3)
##       MSE    R2
## OLS 1.736 0.838
## IV  1.896 0.823

Answer:
The OLS model has a lower MSE and a higher r^2 than the IV model, indicating better model fit. This is expected because OLS selects the coefficients that minimize the sum of squared residuals. However, the better fit does not mean that the OLS estimate is causally correct. The OLS coefficient is inconsistent because x is endogenous, whereas the IV coefficient is closer to the true causal effect despite its slightly worse fit.

Assignment 8 - system estimator vs. 2SLS

Create a data set with endogeneity and instruments; implement a systems-of-equations estimator and compare the estimate to IV / 2SLS.

library(systemfit)
### systems-of-equations estimator
# Define the two equations
main_equation = y ~ x
first_equation = x ~ z
system_model = systemfit(list(
    main = main_equation,
    first = first_equation
  ),
  method = "SUR",
  data = data678,
  maxit = 500
)
summary(system_model)
## 
## systemfit results 
## method: iterated SUR 
## 
## convergence achieved after 11 iterations
## 
##           N   DF     SSR detRCov   OLS-R2 McElroy-R2
## system 2000 1996 3917.85 2.84941 0.768877   0.694226
## 
##          N  DF     SSR     MSE    RMSE       R2   Adj R2
## main  1000 998 1895.57 1.89937 1.37818 0.823232 0.823054
## first 1000 998 2022.28 2.02634 1.42349 0.675288 0.674963
## 
## The covariance matrix of the residuals used for estimation
##           main    first
## main  1.899355 0.999662
## first 0.999662 2.026336
## 
## The covariance matrix of the residuals
##           main    first
## main  1.899369 0.999676
## first 0.999676 2.026336
## 
## The correlations of the residuals
##           main    first
## main  1.000000 0.509564
## first 0.509564 1.000000
## 
## 
## SUR estimates for 'main' (equation 1)
## Model Formula: y ~ x
## 
##              Estimate Std. Error  t value Pr(>|t|)    
## (Intercept) 0.0310598  0.0435910  0.71253   0.4763    
## x           1.0411121  0.0165467 62.91957   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.378176 on 998 degrees of freedom
## Number of observations: 1000 Degrees of Freedom: 998 
## SSR: 1895.570133 MSE: 1.899369 Root MSE: 1.378176 
## Multiple R-Squared: 0.823232 Adjusted R-Squared: 0.823054 
## 
## 
## SUR estimates for 'first' (equation 2)
## Model Formula: x ~ z
## 
##              Estimate Std. Error  t value Pr(>|t|)    
## (Intercept) 0.0212402  0.0450202  0.47179  0.63718    
## z           2.0689826  0.0430304 48.08187  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.423494 on 998 degrees of freedom
## Number of observations: 1000 Degrees of Freedom: 998 
## SSR: 2022.282931 MSE: 2.026336 Root MSE: 1.423494 
## Multiple R-Squared: 0.675288 Adjusted R-Squared: 0.674963
### 2SLS
library(AER)
model_iv8 = ivreg(y ~ x, instruments = ~ z, data = data678)
summary(model_iv8)
## 
## Call:
## ivreg(formula = y ~ x | z, data = data678)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.4124 -0.8858  0.0114  0.9464  4.0646 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.03106    0.04360   0.712    0.476    
## x            1.04111    0.02125  48.990   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.378 on 998 degrees of freedom
## Multiple R-Squared: 0.8232,  Adjusted R-squared: 0.8231 
## Wald test:  2400 on 1 and 998 DF,  p-value: < 2.2e-16
### comparison
comparison8 = rbind( True_value = c(estimate = 1),
                     System_estimator = c(estimate = unname(coef(system_model)["main_x"])),
                     IV_2SLS = c(estimate = unname(coef(model_iv8)["x"])))
round(comparison8, 3)
##                  estimate
## True_value          1.000
## System_estimator    1.041
## IV_2SLS             1.041

Answer:
The system estimator and 2SLS produce approximately the same estimate of the coefficient on x. Both of them are close to the true value of 1. The system estimator jointly estimates the outcome equation and the first-stage equation while allowing their errors to be correlated. In this model, the system estimator and 2SLS address the same endogeneity problem and therefore produce equivalent results.

Assignment 9 - fixed effects and instruments

1. Create a panel data set with unobserved effects in both dimensions and a regressor correlated with the idiosyncratic error.

set.seed(123)
N = 100 # number of individuals
T = 100 # number of time periods
id = rep(1:N, each = T)
time = rep(1:T, times = N)

# unobserved individual effects
alpha_i = rnorm(N)
alpha = alpha_i[id]
# unobserved time effects
lambda_t = rnorm(T)
lambda = lambda_t[time]

z = rnorm(N*T) # instrument 
epsilon = rnorm(N*T) # idiosyncratic error
u = rnorm(N*T) # independent variation in x

x = z + 0.5 * alpha + 0.5 * lambda + 0.5 * epsilon + u
y = 2 + 1 * x + alpha + lambda + epsilon

panel9 = data.frame(id,time,y,x,z,alpha,lambda,epsilon)

2. Implement: cross-section FE only; time FE only; both; instruments without FE; both FE plus instruments. Compare the results.

library(fixest)
## cross-section FE only
cross_fe = feols(y ~ x, data = panel9, fixef = "id")
summary(cross_fe)
## OLS estimation, Dep. Var.: y
## Observations: 10,000
## Fixed-effects: id: 100
## Standard-errors: IID 
##   Estimate Std. Error t value  Pr(>|t|)    
## x   1.3846   0.007971 173.706 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## RMSE: 1.25342     Adj. R2: 0.808777
##                 Within R2: 0.752974
## time FE only
time_fe = feols(y ~ x,data = panel9, fixef = "time")
summary(time_fe)
## OLS estimation, Dep. Var.: y
## Observations: 10,000
## Fixed-effects: time: 100
## Standard-errors: IID 
##   Estimate Std. Error t value  Pr(>|t|)    
## x  1.37228   0.007692 178.412 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## RMSE: 1.21004     Adj. R2: 0.821785
##                 Within R2: 0.762783
## both FE
both_fe = feols(y ~ x,data = panel9,fixef = c("id", "time"))
summary(both_fe)
## OLS estimation, Dep. Var.: y
## Observations: 10,000
## Fixed-effects: id: 100,  time: 100
## Standard-errors: IID 
##   Estimate Std. Error t value  Pr(>|t|)    
## x  1.22043   0.006307 193.495 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## RMSE: 0.937139     Adj. R2: 0.892027
##                  Within R2: 0.79255
## instruments without FE
first_stage = lm(x ~ z, data = panel9)
panel9$x_hat = fitted(first_stage)
iv_no_fe = lm(y ~ x_hat, data = panel9)
summary(iv_no_fe)
## 
## Call:
## lm(formula = y ~ x_hat, data = panel9)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -11.0217  -1.8481  -0.0291   1.7734  10.5408 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.97235    0.02691   73.29   <2e-16 ***
## x_hat        1.00601    0.02629   38.27   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.691 on 9998 degrees of freedom
## Multiple R-squared:  0.1278, Adjusted R-squared:  0.1277 
## F-statistic:  1464 on 1 and 9998 DF,  p-value: < 2.2e-16
## both FE plus instruments
first_stage_fe = lm(x ~ z + factor(id) + factor(time), data = panel9)
panel9$x_hat_fe = fitted(first_stage_fe)
both_fe_iv = lm(y ~ x_hat_fe + factor(id) + factor(time),data = panel9)
summary(both_fe_iv)
## 
## Call:
## lm(formula = y ~ x_hat_fe + factor(id) + factor(time), data = panel9)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.6434 -1.2200  0.0003  1.2101  6.1581 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      0.7763148  0.2555571   3.038 0.002390 ** 
## x_hat_fe         1.0022220  0.0178785  56.057  < 2e-16 ***
## factor(id)2      0.1451262  0.2558042   0.567 0.570501    
## factor(id)3      1.9179867  0.2568082   7.469 8.80e-14 ***
## factor(id)4      0.5654199  0.2557729   2.211 0.027084 *  
## factor(id)5      0.4046992  0.2557730   1.582 0.113623    
## factor(id)6      2.0560509  0.2567759   8.007 1.31e-15 ***
## factor(id)7      0.9045216  0.2559029   3.535 0.000410 ***
## factor(id)8     -1.0024814  0.2559801  -3.916 9.05e-05 ***
## factor(id)9     -0.3203748  0.2557771  -1.253 0.210398    
## factor(id)10     0.0005020  0.2557781   0.002 0.998434    
## factor(id)11     1.6800400  0.2559941   6.563 5.55e-11 ***
## factor(id)12     0.7747093  0.2558985   3.027 0.002473 ** 
## factor(id)13     0.7508198  0.2559956   2.933 0.003365 ** 
## factor(id)14     0.5084082  0.2558272   1.987 0.046916 *  
## factor(id)15    -0.2869538  0.2557731  -1.122 0.261929    
## factor(id)16     2.0641808  0.2565886   8.045 9.65e-16 ***
## factor(id)17     0.8047391  0.2558843   3.145 0.001666 ** 
## factor(id)18    -1.4722642  0.2561949  -5.747 9.37e-09 ***
## factor(id)19     0.9775991  0.2558195   3.821 0.000133 ***
## factor(id)20    -0.1336828  0.2558115  -0.523 0.601276    
## factor(id)21    -0.7306768  0.2559594  -2.855 0.004317 ** 
## factor(id)22    -0.0191228  0.2557896  -0.075 0.940407    
## factor(id)23    -0.7906647  0.2559430  -3.089 0.002012 ** 
## factor(id)24    -0.3982651  0.2557729  -1.557 0.119478    
## factor(id)25    -0.2607373  0.2557934  -1.019 0.308073    
## factor(id)26    -1.3041238  0.2562171  -5.090 3.65e-07 ***
## factor(id)27     1.0169009  0.2559524   3.973 7.15e-05 ***
## factor(id)28     0.4171080  0.2558498   1.630 0.103073    
## factor(id)29    -0.5786693  0.2559645  -2.261 0.023797 *  
## factor(id)30     1.7207908  0.2561398   6.718 1.94e-11 ***
## factor(id)31     0.7351670  0.2558896   2.873 0.004075 ** 
## factor(id)32     0.1558837  0.2558127   0.609 0.542296    
## factor(id)33     1.1880502  0.2559719   4.641 3.51e-06 ***
## factor(id)34     1.2469041  0.2560794   4.869 1.14e-06 ***
## factor(id)35     1.1679232  0.2560211   4.562 5.13e-06 ***
## factor(id)36     0.9894457  0.2558136   3.868 0.000111 ***
## factor(id)37     0.9830715  0.2560363   3.840 0.000124 ***
## factor(id)38     0.2956637  0.2558473   1.156 0.247862    
## factor(id)39     0.1862849  0.2557775   0.728 0.466442    
## factor(id)40    -0.0716708  0.2557742  -0.280 0.779321    
## factor(id)41    -0.3571196  0.2558621  -1.396 0.162821    
## factor(id)42     0.2240041  0.2557729   0.876 0.381164    
## factor(id)43    -0.9104289  0.2559408  -3.557 0.000377 ***
## factor(id)44     2.5724670  0.2568015  10.017  < 2e-16 ***
## factor(id)45     1.4545295  0.2561750   5.678 1.40e-08 ***
## factor(id)46    -0.7971452  0.2558562  -3.116 0.001841 ** 
## factor(id)47    -0.0520510  0.2558019  -0.203 0.838763    
## factor(id)48    -0.0113990  0.2558029  -0.045 0.964458    
## factor(id)49     1.0683868  0.2560922   4.172 3.05e-05 ***
## factor(id)50     0.4472699  0.2558551   1.748 0.080472 .  
## factor(id)51     0.6859522  0.2558106   2.681 0.007342 ** 
## factor(id)52     0.3675579  0.2557812   1.437 0.150750    
## factor(id)53     0.3778611  0.2558166   1.477 0.139687    
## factor(id)54     1.7295269  0.2562309   6.750 1.56e-11 ***
## factor(id)55     0.0089951  0.2557987   0.035 0.971949    
## factor(id)56     1.7463428  0.2562524   6.815 9.99e-12 ***
## factor(id)57    -1.1766256  0.2559808  -4.597 4.35e-06 ***
## factor(id)58     0.8831853  0.2559478   3.451 0.000562 ***
## factor(id)59     0.4619726  0.2558140   1.806 0.070966 .  
## factor(id)60     0.5861632  0.2557749   2.292 0.021943 *  
## factor(id)61     0.6015598  0.2558196   2.352 0.018718 *  
## factor(id)62    -0.2796642  0.2558738  -1.093 0.274431    
## factor(id)63     0.0307602  0.2557739   0.120 0.904277    
## factor(id)64    -0.6172479  0.2558009  -2.413 0.015840 *  
## factor(id)65    -0.7804236  0.2559257  -3.049 0.002299 ** 
## factor(id)66     0.5327045  0.2557800   2.083 0.037307 *  
## factor(id)67     0.7750527  0.2558913   3.029 0.002461 ** 
## factor(id)68     0.3066210  0.2557865   1.199 0.230659    
## factor(id)69     1.1874097  0.2559367   4.639 3.54e-06 ***
## factor(id)70     2.4637198  0.2567149   9.597  < 2e-16 ***
## factor(id)71    -0.1735975  0.2557780  -0.679 0.497342    
## factor(id)72    -2.0801075  0.2566350  -8.105 5.89e-16 ***
## factor(id)73     1.3382639  0.2561065   5.225 1.77e-07 ***
## factor(id)74    -0.3353793  0.2558484  -1.311 0.189939    
## factor(id)75    -0.3137262  0.2557757  -1.227 0.220015    
## factor(id)76     1.4813368  0.2560847   5.785 7.49e-09 ***
## factor(id)77     0.0008618  0.2557772   0.003 0.997312    
## factor(id)78    -0.8847305  0.2559439  -3.457 0.000549 ***
## factor(id)79     0.4775223  0.2558110   1.867 0.061973 .  
## factor(id)80     0.2458910  0.2557926   0.961 0.336430    
## factor(id)81     0.2713321  0.2558173   1.061 0.288876    
## factor(id)82     0.7360713  0.2559777   2.876 0.004042 ** 
## factor(id)83    -0.0830349  0.2557729  -0.325 0.745458    
## factor(id)84     0.9798296  0.2558739   3.829 0.000129 ***
## factor(id)85     0.0858638  0.2557750   0.336 0.737104    
## factor(id)86     0.4635754  0.2558006   1.812 0.069978 .  
## factor(id)87     1.6416345  0.2563278   6.404 1.58e-10 ***
## factor(id)88     0.7854998  0.2558095   3.071 0.002142 ** 
## factor(id)89     0.1155591  0.2557731   0.452 0.651421    
## factor(id)90     1.5032261  0.2560018   5.872 4.45e-09 ***
## factor(id)91     1.1946299  0.2560184   4.666 3.11e-06 ***
## factor(id)92     0.9452376  0.2558812   3.694 0.000222 ***
## factor(id)93     0.5675730  0.2558113   2.219 0.026529 *  
## factor(id)94    -0.3413048  0.2558035  -1.334 0.182154    
## factor(id)95     1.6324791  0.2562516   6.371 1.97e-10 ***
## factor(id)96    -0.2542402  0.2557819  -0.994 0.320261    
## factor(id)97     2.5247062  0.2566571   9.837  < 2e-16 ***
## factor(id)98     1.6770521  0.2563253   6.543 6.34e-11 ***
## factor(id)99     0.1283873  0.2557948   0.502 0.615738    
## factor(id)100   -0.8641026  0.2558262  -3.378 0.000734 ***
## factor(time)2    1.1571046  0.2560159   4.520 6.27e-06 ***
## factor(time)3    0.7166164  0.2559214   2.800 0.005118 ** 
## factor(time)4    0.4858013  0.2558266   1.899 0.057601 .  
## factor(time)5   -0.0103191  0.2557955  -0.040 0.967822    
## factor(time)6    0.8369784  0.2561578   3.267 0.001089 ** 
## factor(time)7    0.1631250  0.2558197   0.638 0.523713    
## factor(time)8   -0.7963015  0.2557731  -3.113 0.001855 ** 
## factor(time)9    0.5009122  0.2560429   1.956 0.050451 .  
## factor(time)10   1.6853281  0.2565673   6.569 5.33e-11 ***
## factor(time)11   0.3648228  0.2558570   1.426 0.153933    
## factor(time)12   1.5146003  0.2563335   5.909 3.56e-09 ***
## factor(time)13  -0.9165377  0.2557756  -3.583 0.000341 ***
## factor(time)14   0.8548789  0.2562049   3.337 0.000851 ***
## factor(time)15   1.3957458  0.2561443   5.449 5.19e-08 ***
## factor(time)16   1.0966334  0.2559684   4.284 1.85e-05 ***
## factor(time)17   0.9988947  0.2561005   3.900 9.67e-05 ***
## factor(time)18   0.1875669  0.2557829   0.733 0.463390    
## factor(time)19   0.0096205  0.2558476   0.038 0.970005    
## factor(time)20  -0.1392360  0.2558242  -0.544 0.586272    
## factor(time)21   1.2081749  0.2561275   4.717 2.43e-06 ***
## factor(time)22  -0.1524845  0.2557958  -0.596 0.551110    
## factor(time)23   0.3810552  0.2559807   1.489 0.136623    
## factor(time)24   0.6545223  0.2558744   2.558 0.010543 *  
## factor(time)25   2.6740050  0.2570760  10.402  < 2e-16 ***
## factor(time)26   0.2752993  0.2559275   1.076 0.282091    
## factor(time)27   1.1106805  0.2562489   4.334 1.48e-05 ***
## factor(time)28   0.8324993  0.2559935   3.252 0.001150 ** 
## factor(time)29  -0.1984980  0.2557797  -0.776 0.437738    
## factor(time)30   0.6880518  0.2558527   2.689 0.007173 ** 
## factor(time)31   2.5399138  0.2575255   9.863  < 2e-16 ***
## factor(time)32   1.3341635  0.2560663   5.210 1.92e-07 ***
## factor(time)33   0.9288413  0.2558429   3.631 0.000284 ***
## factor(time)34   0.4392729  0.2558322   1.717 0.086004 .  
## factor(time)35  -1.2182415  0.2560184  -4.758 1.98e-06 ***
## factor(time)36   1.9292161  0.2563872   7.525 5.75e-14 ***
## factor(time)37  -0.5203726  0.2557734  -2.035 0.041927 *  
## factor(time)38   1.6704885  0.2564282   6.514 7.65e-11 ***
## factor(time)39   2.8479812  0.2574546  11.062  < 2e-16 ***
## factor(time)40  -0.5238722  0.2557730  -2.048 0.040568 *  
## factor(time)41   1.5706503  0.2565369   6.123 9.56e-10 ***
## factor(time)42   0.5868674  0.2560512   2.292 0.021927 *  
## factor(time)43  -0.6765061  0.2557865  -2.645 0.008187 ** 
## factor(time)44  -0.6696915  0.2557744  -2.618 0.008851 ** 
## factor(time)45  -0.6909974  0.2558077  -2.701 0.006920 ** 
## factor(time)46   0.2931427  0.2557767   1.146 0.251786    
## factor(time)47  -0.6438702  0.2557736  -2.517 0.011840 *  
## factor(time)48   1.4405983  0.2560913   5.625 1.90e-08 ***
## factor(time)49   2.9340321  0.2577291  11.384  < 2e-16 ***
## factor(time)50  -0.4400699  0.2557782  -1.721 0.085371 .  
## factor(time)51   1.7830720  0.2565159   6.951 3.86e-12 ***
## factor(time)52   1.6610133  0.2567830   6.469 1.04e-10 ***
## factor(time)53   0.9607382  0.2561429   3.751 0.000177 ***
## factor(time)54   0.0082532  0.2558860   0.032 0.974271    
## factor(time)55   0.8068946  0.2561001   3.151 0.001634 ** 
## factor(time)56   0.7036125  0.2558866   2.750 0.005976 ** 
## factor(time)57   1.3913497  0.2563832   5.427 5.87e-08 ***
## factor(time)58   0.4650256  0.2560221   1.816 0.069347 .  
## factor(time)59   1.8391241  0.2564256   7.172 7.92e-13 ***
## factor(time)60   0.3624403  0.2558515   1.417 0.156631    
## factor(time)61   2.0662263  0.2568994   8.043 9.79e-16 ***
## factor(time)62  -0.1564624  0.2558230  -0.612 0.540814    
## factor(time)63  -0.5667277  0.2557757  -2.216 0.026733 *  
## factor(time)64   4.1712199  0.2594536  16.077  < 2e-16 ***
## factor(time)65   0.3266899  0.2557875   1.277 0.201564    
## factor(time)66   1.1795032  0.2562818   4.602 4.23e-06 ***
## factor(time)67   1.2248593  0.2562163   4.781 1.77e-06 ***
## factor(time)68   0.5299221  0.2559233   2.071 0.038420 *  
## factor(time)69   1.3255614  0.2561946   5.174 2.34e-07 ***
## factor(time)70   1.1825784  0.2561670   4.616 3.95e-06 ***
## factor(time)71   0.6086783  0.2559610   2.378 0.017425 *  
## factor(time)72   1.0859152  0.2560920   4.240 2.25e-05 ***
## factor(time)73   0.7793555  0.2559082   3.045 0.002330 ** 
## factor(time)74   3.0868171  0.2575587  11.985  < 2e-16 ***
## factor(time)75   0.0999966  0.2557957   0.391 0.695862    
## factor(time)76  -0.3943831  0.2558045  -1.542 0.123170    
## factor(time)77   0.8104569  0.2562108   3.163 0.001565 ** 
## factor(time)78   1.2279323  0.2562796   4.791 1.68e-06 ***
## factor(time)79   1.1475352  0.2560885   4.481 7.51e-06 ***
## factor(time)80   0.4142930  0.2559221   1.619 0.105517    
## factor(time)81  -0.1284706  0.2557794  -0.502 0.615488    
## factor(time)82   2.1849358  0.2569386   8.504  < 2e-16 ***
## factor(time)83   0.5285209  0.2559571   2.065 0.038961 *  
## factor(time)84   0.0538158  0.2558250   0.210 0.833390    
## factor(time)85   0.6689989  0.2558862   2.614 0.008951 ** 
## factor(time)86   0.8015937  0.2560087   3.131 0.001747 ** 
## factor(time)87   2.2772019  0.2564630   8.879  < 2e-16 ***
## factor(time)88   1.0099097  0.2560613   3.944 8.07e-05 ***
## factor(time)89   1.5582719  0.2564283   6.077 1.27e-09 ***
## factor(time)90   0.4038468  0.2558316   1.579 0.114468    
## factor(time)91   1.0701874  0.2561507   4.178 2.97e-05 ***
## factor(time)92   0.3912242  0.2559206   1.529 0.126373    
## factor(time)93   1.0060319  0.2564630   3.923 8.82e-05 ***
## factor(time)94  -0.1638000  0.2559157  -0.640 0.522152    
## factor(time)95  -0.5346546  0.2557731  -2.090 0.036612 *  
## factor(time)96   2.6945978  0.2570763  10.482  < 2e-16 ***
## factor(time)97   1.6006994  0.2563218   6.245 4.42e-10 ***
## factor(time)98  -0.2863655  0.2557729  -1.120 0.262908    
## factor(time)99   0.4577692  0.2558826   1.789 0.073649 .  
## factor(time)100 -0.1251680  0.2557735  -0.489 0.624590    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.809 on 9800 degrees of freedom
## Multiple R-squared:  0.6137, Adjusted R-squared:  0.6059 
## F-statistic: 78.25 on 199 and 9800 DF,  p-value: < 2.2e-16
### comparison
comparison9 = c(Cross_section_FE = coef(cross_fe)["x"], 
                Time_FE = coef(time_fe)["x"],
                Both_FE = coef(both_fe)["x"],
                IV_without_FE = coef(iv_no_fe)["x_hat"],
                Both_FE_plus_IV = coef(both_fe_iv)["x_hat_fe"])
round(comparison9, 3)
##       Cross_section_FE.x                Time_FE.x                Both_FE.x 
##                    1.385                    1.372                    1.220 
##      IV_without_FE.x_hat Both_FE_plus_IV.x_hat_fe 
##                    1.006                    1.002

Answer:
The cross-section FE, time FE, and two-way FE estimates are 1.385, 1.372, and 1.220, respectively, all above the true coefficient of 1. Fixed effects reduce some of the bias but cannot eliminate the correlation between x and the idiosyncratic error. In contrast, IV without FE and IV with both FE produce estimates of 1.006 and 1.002, respectively. Both are close to the true value, showing that an instrument is needed to address the remaining endogeneity.

Assignment 10 - the return to schooling (Card 1995)

1. Obtain Card’s data (e.g., davidcard.berkeley.edu; in R, the wooldridge package). Estimate the OLS and IV equations from Card (1995). (You will not get identical estimates.) Report and discuss the first stage.

## The key question of this assignment is: What is the causal effect of 1 additional year of education on salary?
library(wooldridge)
data("card")
df_card = card
# dim(card); names(card)

#### pick variables for modelling
#### delete id, wage
#### Based on the official documentation examples of ivpack from CRAN
variables = c("lwage","educ", "nearc4","exper","expersq","black","south","smsa","smsa66","reg661",
              "reg662","reg663","reg664","reg665","reg666","reg667","reg668")
# summary(df_card[,variables])
# colSums(is.na(df_card[, variables])) # check NA value

df_card_model = df_card[,variables]

############################ 1. estimate the OLS model
ols_card = lm(lwage ~ educ + exper + expersq + black + south + smsa + smsa66 + reg661 + reg662 + reg663
              + reg664 + reg665 + reg666 + reg667 + reg668, data = df_card_model )
summary(ols_card)
## 
## Call:
## lm(formula = lwage ~ educ + exper + expersq + black + south + 
##     smsa + smsa66 + reg661 + reg662 + reg663 + reg664 + reg665 + 
##     reg666 + reg667 + reg668, data = df_card_model)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.62326 -0.22141  0.02001  0.23932  1.33340 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  4.7393766  0.0715282  66.259  < 2e-16 ***
## educ         0.0746933  0.0034983  21.351  < 2e-16 ***
## exper        0.0848320  0.0066242  12.806  < 2e-16 ***
## expersq     -0.0022870  0.0003166  -7.223 6.41e-13 ***
## black       -0.1990123  0.0182483 -10.906  < 2e-16 ***
## south       -0.1479550  0.0259799  -5.695 1.35e-08 ***
## smsa         0.1363845  0.0201005   6.785 1.39e-11 ***
## smsa66       0.0262417  0.0194477   1.349 0.177327    
## reg661      -0.1185698  0.0388301  -3.054 0.002281 ** 
## reg662      -0.0222026  0.0282575  -0.786 0.432092    
## reg663       0.0259703  0.0273644   0.949 0.342670    
## reg664      -0.0634942  0.0356803  -1.780 0.075254 .  
## reg665       0.0094551  0.0361174   0.262 0.793503    
## reg666       0.0219476  0.0400984   0.547 0.584182    
## reg667      -0.0005887  0.0393793  -0.015 0.988073    
## reg668      -0.1750058  0.0463394  -3.777 0.000162 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3723 on 2994 degrees of freedom
## Multiple R-squared:  0.2998, Adjusted R-squared:  0.2963 
## F-statistic: 85.48 on 15 and 2994 DF,  p-value: < 2.2e-16
coef(summary(ols_card))["educ", ]  # 0.075
##     Estimate   Std. Error      t value     Pr(>|t|) 
## 7.469326e-02 3.498346e-03 2.135102e+01 2.892621e-94

Answer:

2. Is the instrument strong? Is it plausibly exogenous?

In OLS, one additional year of education is associated with approximately a 7.5% increase in wages.

############################ 2. estimate the first stage model
first_stage = lm(educ ~ nearc4 + exper + expersq + black + south + smsa + smsa66 + reg661 + reg662 
                 + reg663 + reg664 + reg665 + reg666 + reg667 + reg668, data = df_card_model)
summary(first_stage)
## 
## Call:
## lm(formula = educ ~ nearc4 + exper + expersq + black + south + 
##     smsa + smsa66 + reg661 + reg662 + reg663 + reg664 + reg665 + 
##     reg666 + reg667 + reg668, data = df_card_model)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -7.545 -1.370 -0.091  1.278  6.239 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 16.8485239  0.2111222  79.805  < 2e-16 ***
## nearc4       0.3198989  0.0878638   3.641 0.000276 ***
## exper       -0.4125334  0.0336996 -12.241  < 2e-16 ***
## expersq      0.0008686  0.0016504   0.526 0.598728    
## black       -0.9355287  0.0937348  -9.981  < 2e-16 ***
## south       -0.0516126  0.1354284  -0.381 0.703152    
## smsa         0.4021825  0.1048112   3.837 0.000127 ***
## smsa66       0.0254805  0.1057692   0.241 0.809644    
## reg661      -0.2102710  0.2024568  -1.039 0.299076    
## reg662      -0.2889073  0.1473395  -1.961 0.049992 *  
## reg663      -0.2382099  0.1426357  -1.670 0.095012 .  
## reg664      -0.0930890  0.1859827  -0.501 0.616742    
## reg665      -0.4828875  0.1881872  -2.566 0.010336 *  
## reg666      -0.5130857  0.2096352  -2.448 0.014442 *  
## reg667      -0.4270887  0.2056208  -2.077 0.037880 *  
## reg668       0.3136204  0.2416739   1.298 0.194490    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.941 on 2994 degrees of freedom
## Multiple R-squared:  0.4771, Adjusted R-squared:  0.4745 
## F-statistic: 182.1 on 15 and 2994 DF,  p-value: < 2.2e-16
coef(summary(first_stage))["nearc4", ]  ## 0.32
##     Estimate   Std. Error      t value     Pr(>|t|) 
## 0.3198989401 0.0878638178 3.6408495342 0.0002763401
nearc4_t = coef(summary(first_stage))["nearc4", "t value"]
first_stage_F = nearc4_t^2
first_stage_F # calculate F: 13.25579
## [1] 13.25579

Answer:
First-stage result: Growing up near a four-year college is positively associated with educational attainment. The coefficient on nearc4 is positive and statistically significant. The first-stage F-statistic exceeds the conventional threshold of 10, suggesting that nearc4 is not a weak instrument in this specification.

3. Compare the OLS and IV coefficients for the focal regressor. Is the difference in the expected direction? What does it imply?

############################ 3. estimate the IV model
library(systemfit)
iv_card = systemfit(lwage ~ educ + exper + expersq + black + south + smsa + smsa66 + reg661 + reg662 
                    + reg663 + reg664 + reg665 + reg666 + reg667 + reg668, method = "2SLS",
                    inst = ~ nearc4 + exper + expersq + black + south + smsa + smsa66 + reg661 + reg662
                    + reg663 + reg664 + reg665 + reg666 + reg667 + reg668,
                    data = df_card_model)
summary(iv_card)
## 
## systemfit results 
## method: 2SLS 
## 
##           N   DF     SSR detRCov   OLS-R2 McElroy-R2
## system 3010 2994 451.495  0.1508 0.238166   0.238166
## 
##        N   DF     SSR    MSE    RMSE       R2   Adj R2
## eq1 3010 2994 451.495 0.1508 0.38833 0.238166 0.234349
## 
## The covariance matrix of the residuals
##        eq1
## eq1 0.1508
## 
## The correlations of the residuals
##     eq1
## eq1   1
## 
## 
## 2SLS estimates for 'eq1' (equation 1)
## Model Formula: lwage ~ educ + exper + expersq + black + south + smsa + smsa66 + 
##     reg661 + reg662 + reg663 + reg664 + reg665 + reg666 + reg667 + 
##     reg668
## Instruments: ~nearc4 + exper + expersq + black + south + smsa + smsa66 + reg661 + 
##     reg662 + reg663 + reg664 + reg665 + reg666 + reg667 + reg668
## 
##                 Estimate   Std. Error  t value   Pr(>|t|)    
## (Intercept)  3.773965141  0.934947017  4.03656 5.5601e-05 ***
## educ         0.131503836  0.054963673  2.39256 0.01679262 *  
## exper        0.108271106  0.023658571  4.57640 4.9229e-06 ***
## expersq     -0.002334938  0.000333497 -7.00137 3.1162e-12 ***
## black       -0.146775747  0.053899859 -2.72312 0.00650438 ** 
## south       -0.144671501  0.027284623 -5.30231 1.2266e-07 ***
## smsa         0.111808309  0.031661988  3.53131 0.00041975 ***
## smsa66       0.018531104  0.021608589  0.85758 0.39119279    
## reg661      -0.107814233  0.041813679 -2.57844 0.00997198 ** 
## reg662      -0.007046452  0.032907270 -0.21413 0.83045984    
## reg663       0.040444546  0.031780563  1.27262 0.20325211    
## reg664      -0.057917154  0.037605896 -1.54011 0.12363963    
## reg665       0.038457680  0.046938701  0.81932 0.41267073    
## reg666       0.055088709  0.052659734  1.04613 0.29558738    
## reg667       0.026757977  0.048828702  0.54800 0.58373489    
## reg668      -0.190891226  0.050711333 -3.76427 0.00017024 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.38833 on 2994 degrees of freedom
## Number of observations: 3010 Degrees of Freedom: 2994 
## SSR: 451.494832 MSE: 0.1508 Root MSE: 0.38833 
## Multiple R-Squared: 0.238166 Adjusted R-Squared: 0.234349
iv_educ = coef(iv_card$eq[[1]])["educ"]
iv_educ  # 0.1315038
##      educ 
## 0.1315038
# compare
ols_educ = coef(ols_card)["educ"]
coefficient_comparison = c(OLS = ols_educ, IV_2SLS = iv_educ)
round(coefficient_comparison, 3)
##     OLS.educ IV_2SLS.educ 
##        0.075        0.132

Answer:
The IV estimate of 0.132 is higher than the OLS estimate of 0.075, so the difference is in the expected positive direction. This suggests that OLS may underestimate the return to education, possibly because of measurement error in schooling. It may also reflect heterogeneous returns, since the IV estimate identifies the return for individuals whose educational attainment is affected by proximity to a four-year college. Therefore, the IV estimate should be interpreted as a LATE rather than necessarily the average return for the entire population.

4. Are the results robust to reasonable changes in the covariate set?

##### add family background related variables: fatheduc，motheduc，momdad14，sinmom14，step14，libcrd14 
iv_card2 = systemfit(lwage ~ educ + exper + expersq + black + south + smsa + smsa66 + reg661 + reg662 
                    + reg663 + reg664 + reg665 + reg666 + reg667 + reg668, method = "2SLS",
                    inst = ~ nearc4 + exper + expersq + black + south + smsa + smsa66 + reg661 + reg662
                    + reg663 + reg664 + reg665 + reg666 + reg667 + reg668 + fatheduc + motheduc 
                    + momdad14 + sinmom14 + step14 + libcrd14,
                    data = df_card)
summary(iv_card2)
## 
## systemfit results 
## method: 2SLS 
## 
##           N   DF    SSR  detRCov   OLS-R2 McElroy-R2
## system 2216 2200 318.61 0.144823 0.256177   0.256177
## 
##        N   DF    SSR      MSE     RMSE       R2   Adj R2
## eq1 2216 2200 318.61 0.144823 0.380556 0.256177 0.251106
## 
## The covariance matrix of the residuals
##          eq1
## eq1 0.144823
## 
## The correlations of the residuals
##     eq1
## eq1   1
## 
## 
## 2SLS estimates for 'eq1' (equation 1)
## Model Formula: lwage ~ educ + exper + expersq + black + south + smsa + smsa66 + 
##     reg661 + reg662 + reg663 + reg664 + reg665 + reg666 + reg667 + 
##     reg668
## Instruments: ~nearc4 + exper + expersq + black + south + smsa + smsa66 + reg661 + 
##     reg662 + reg663 + reg664 + reg665 + reg666 + reg667 + reg668 + 
##     fatheduc + motheduc + momdad14 + sinmom14 + step14 + libcrd14
## 
##                 Estimate   Std. Error  t value   Pr(>|t|)    
## (Intercept)  4.144188125  0.204086800 20.30601 < 2.22e-16 ***
## educ         0.106862600  0.011609130  9.20505 < 2.22e-16 ***
## exper        0.102737533  0.009248186 11.10894 < 2.22e-16 ***
## expersq     -0.002505686  0.000402047 -6.23232 5.4918e-10 ***
## black       -0.150686671  0.025948245 -5.80720 7.2740e-09 ***
## south       -0.122298589  0.031592911 -3.87108 0.00011151 ***
## smsa         0.122349834  0.024506963  4.99245 6.4300e-07 ***
## smsa66       0.028040769  0.022959575  1.22131 0.22209936    
## reg661      -0.079285192  0.046301523 -1.71237 0.08697001 .  
## reg662       0.009775613  0.032451554  0.30124 0.76326216    
## reg663       0.043143966  0.031734670  1.35952 0.17412073    
## reg664      -0.051283106  0.041381442 -1.23928 0.21537479    
## reg665       0.015396690  0.042902133  0.35888 0.71971978    
## reg666       0.034657729  0.049698865  0.69735 0.48565463    
## reg667       0.020027778  0.046737522  0.42852 0.66831741    
## reg668      -0.160611840  0.052703700 -3.04745 0.00233537 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.380556 on 2200 degrees of freedom
## Number of observations: 2216 Degrees of Freedom: 2200 
## SSR: 318.610262 MSE: 0.144823 Root MSE: 0.380556 
## Multiple R-Squared: 0.256177 Adjusted R-Squared: 0.251106
iv_educ2 = coef(iv_card2$eq[[1]])["educ"]
iv_educ2  # 0.1068626
##      educ 
## 0.1068626
# compare
coefficient_comparison2 = c(OLS = ols_educ, IV_2SLS = iv_educ, IV_2SLS_add = iv_educ2)
round(coefficient_comparison2, 3)
##         OLS.educ     IV_2SLS.educ IV_2SLS_add.educ 
##            0.075            0.132            0.107

Answer:
The IV estimate of the return to education is approximately 0.132 in the baseline specification and 0.107 after adding family-background and regional covariates. Although both estimates remain positive, the magnitude decreases by approximately 19%. Therefore, the positive relationship is qualitatively robust, but the estimated magnitude is somewhat sensitive to the choice of covariates.

############################ estimate the first stage model again
first_stage2 = lm(educ ~ nearc4 + exper + expersq + black + south + smsa + smsa66 + reg661 + reg662 
                 + reg663 + reg664 + reg665 + reg666 + reg667 + reg668 + fatheduc + motheduc + momdad14 + sinmom14 
                 + step14 + libcrd14, data = df_card)
summary(first_stage2)
## 
## Call:
## lm(formula = educ ~ nearc4 + exper + expersq + black + south + 
##     smsa + smsa66 + reg661 + reg662 + reg663 + reg664 + reg665 + 
##     reg666 + reg667 + reg668 + fatheduc + motheduc + momdad14 + 
##     sinmom14 + step14 + libcrd14, data = df_card)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.4362 -1.3586 -0.1302  1.1978  5.9700 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 14.270491   0.519396  27.475  < 2e-16 ***
## nearc4       0.224151   0.097657   2.295 0.021810 *  
## exper       -0.379486   0.037949 -10.000  < 2e-16 ***
## expersq      0.002747   0.001950   1.409 0.158968    
## black       -0.262456   0.121626  -2.158 0.031045 *  
## south       -0.029616   0.153115  -0.193 0.846648    
## smsa         0.392388   0.115744   3.390 0.000711 ***
## smsa66      -0.240418   0.115834  -2.076 0.038052 *  
## reg661      -0.469587   0.224700  -2.090 0.036747 *  
## reg662      -0.344232   0.156624  -2.198 0.028066 *  
## reg663      -0.388867   0.152616  -2.548 0.010901 *  
## reg664      -0.113098   0.200739  -0.563 0.573214    
## reg665      -0.268505   0.207637  -1.293 0.196097    
## reg666      -0.368364   0.241405  -1.526 0.127176    
## reg667      -0.271679   0.225830  -1.203 0.229097    
## reg668       0.055531   0.255864   0.217 0.828203    
## fatheduc     0.104968   0.014527   7.226 6.84e-13 ***
## motheduc     0.122630   0.017030   7.201 8.18e-13 ***
## momdad14    -0.342701   0.426843  -0.803 0.422135    
## sinmom14    -0.259160   0.618383  -0.419 0.675189    
## step14      -1.439332   0.475112  -3.029 0.002478 ** 
## libcrd14     0.471951   0.099482   4.744 2.23e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.845 on 2194 degrees of freedom
##   (794 observations deleted due to missingness)
## Multiple R-squared:  0.4972, Adjusted R-squared:  0.4923 
## F-statistic: 103.3 on 21 and 2194 DF,  p-value: < 2.2e-16
coef(summary(first_stage2))["nearc4", ]  ## 0.224151
##   Estimate Std. Error    t value   Pr(>|t|) 
## 0.22415133 0.09765664 2.29530049 0.02180986
nearc4_t2 = coef(summary(first_stage2))["nearc4", "t value"]
first_stage_F2 = nearc4_t2^2
first_stage_F2 # calculate F: 5.268404
## [1] 5.268404

Answer:
After adding family-background covariates, the first-stage coefficient on nearc4 is 0.224 and the corresponding F-statistic decreases from 13.26 to 5.27. Thus, nearc4 remains statistically relevant, but its conditional relationship with education is relatively weak.

Assignment 11 - the IV estimator by hand

1. Create a data set in which a non-normal regressor is correlated with a normal error. Implement the endogeneity correction via Gaussian copulas (original Park & Gupta 2012 approach); make sure the model has an intercept.

library(MASS)
library(nortest)

N = 1000
Sigma = matrix(c(1, 0,0, 1),2,2)

# Generate two independent normal errors
errors = mvrnorm(N,mu = c(0, 0),Sigma = Sigma)
# Transform a normal variable into a uniform variable
u = pnorm(errors[, 1])
# Generate a non-normal endogenous regressor
p = qt(u, df = 4)
# Generate a normal error correlated with p
e = errors[, 1] + errors[, 2]

# True intercept = 1 and true coefficient of p = 1
y = 1 + p + e

cor(p, e)
## [1] 0.709856
ad.test(p)
## 
##  Anderson-Darling normality test
## 
## data:  p
## A = 5.6796, p-value = 5.373e-14
summary(lm(y ~ p))
## 
## Call:
## lm(formula = y ~ p)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.8153 -0.6471 -0.0305  0.6986  3.6080 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.96928    0.03134   30.92   <2e-16 ***
## p            1.75575    0.02374   73.97   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.991 on 998 degrees of freedom
## Multiple R-squared:  0.8457, Adjusted R-squared:  0.8456 
## F-statistic:  5471 on 1 and 998 DF,  p-value: < 2.2e-16
# Calculate the empirical distribution of p
u_p = ecdf(p)(p)
# Avoid values equal to exactly 0 or 1
u_p = pmin(pmax(u_p, 1e-6),1 - 1e-6)
# Construct the Gaussian-copula correction term
p_star = qnorm(u_p)
# Copula-corrected model
copula_model = lm(y ~ p + p_star)
summary(copula_model)
## 
## Call:
## lm(formula = y ~ p + p_star)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.8863 -0.6509 -0.0271  0.6743  3.6122 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   0.9819     0.0310  31.668  < 2e-16 ***
## p             1.0815     0.1277   8.468  < 2e-16 ***
## p_star        0.9012     0.1678   5.370 9.77e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9775 on 997 degrees of freedom
## Multiple R-squared:  0.8501, Adjusted R-squared:  0.8498 
## F-statistic:  2826 on 2 and 997 DF,  p-value: < 2.2e-16
# Original biased OLS model
ols_model = lm(y ~ p)
summary(ols_model)
## 
## Call:
## lm(formula = y ~ p)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.8153 -0.6471 -0.0305  0.6986  3.6080 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.96928    0.03134   30.92   <2e-16 ***
## p            1.75575    0.02374   73.97   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.991 on 998 degrees of freedom
## Multiple R-squared:  0.8457, Adjusted R-squared:  0.8456 
## F-statistic:  5471 on 1 and 998 DF,  p-value: < 2.2e-16
cor(p, p_star, method = "spearman")
## [1] 1

Answer:
The original OLS model overestimates the true coefficient of 1 as 1.73 because p is positively correlated with the error term. After adding the Gaussian copula correction term p_star, the estimated coefficient on p decreases to 0.905, which is much closer to the true value. This suggests that the correction term removes most of the endogeneity bias. However, the high correlation between p and p_star (cor = 1) increases the standard error (0.02494 -> 0.15075) of the corrected estimate.

2. Assess robustness: make the error non-normal

# Robustness check 1: non-normal error
set.seed(123)
N = 1000
# z creates endogeneity
z = rnorm(N)
# Generate a non-normal endogenous regressor
u = pnorm(z)
p = qt(u, df = 4)
# Generate a non-normal error
# The gamma variable has mean 2, so subtract 2 to center it
v_non_normal = rgamma(N, shape = 2, rate = 1) - 2
e_non_normal = z + v_non_normal
# True coefficient of p is 1
y = 1 + 1 * p + e_non_normal
# Construct the Gaussian-copula control variable
p_star = qnorm(rank(p, ties.method = "average") / (N + 1))
# OLS and Gaussian-copula models
ols_non_normal_error = lm(y ~ 1 + p)
copula_non_normal_error = lm(y ~ 1 + p + p_star)
# Check that p is endogenous and the error is non-normal
cor(p, e_non_normal)
## [1] 0.5922586
shapiro.test(e_non_normal)
## 
##  Shapiro-Wilk normality test
## 
## data:  e_non_normal
## W = 0.97557, p-value = 6.335e-12
# Compare the estimates
non_normal_error_result = data.frame(
  cor_p_error = cor(p, e_non_normal),
  OLS_beta = unname(coef(ols_non_normal_error)["p"]),
  copula_beta = unname(coef(copula_non_normal_error)["p"]),
  copula_SE = unname(coef(summary(copula_non_normal_error))["p", "Std. Error"])
)

round(non_normal_error_result, 3)

Answer:
The correlation between p and the non-normal error is 0.592, confirming that p is endogenous. The Shapiro–Wilk test strongly rejects the normality of the error term (p<0.001). OLS estimates the coefficient on p as 1.731, substantially above the true value of 1. After adding the Gaussian copula correction term, the estimate decreases to 0.974, with a standard error of 0.166. Thus, the copula correction removes most of the endogeneity bias in this simulation, even when the error term is non-normal.

3. apply different distributions to the endogenous regressor (normal, t, gamma, uniform).

# Robustness check 2: different distributions of p
set.seed(123)
N = 1000
# Generate two independent normal variables
Sigma = matrix(c(1, 0, 0, 1), 2, 2)
errors = mvrnorm(N, mu = c(0, 0), Sigma = Sigma)
# Uniform variable used to generate different distributions
u = pnorm(errors[, 1])
# Normal error correlated with p
e = errors[, 1] + errors[, 2]
# Four distributions of the endogenous regressor
p_set = data.frame(normal = qnorm(u), t = qt(u, df = 4), gamma = qgamma(u, shape = 2, rate = 1),uniform = u)
# Empty table for the results
results = matrix(NA, nrow = 4, ncol = 7)
rownames(results) = c("Normal", "t", "Gamma", "Uniform")
colnames(results) = c("cor_p_error","cor_p_pstar","OLS_beta","OLS_bias","copula_beta","copula_bias","copula_SE")

# Repeat the estimation for each distribution
for (i in 1:4) {
  p = p_set[, i]
  # Standardize p for a fair comparison
  p = as.numeric(scale(p))
  # True coefficient of p is 1
  y = 1 + p + e
  # Gaussian-copula correction term
  p_star = qnorm(rank(p, ties.method = "average") / (N + 1))
  # Estimate OLS and copula models
  ols_model = lm(y ~ p)
  copula_model = lm(y ~ p + p_star)
  # Store the results
  results[i, ] = c(
    cor(p, e),
    cor(p, p_star),
    coef(ols_model)["p"],
    coef(ols_model)["p"] - 1,
    coef(copula_model)["p"],
    coef(copula_model)["p"] - 1,
    coef(summary(copula_model))["p", "Std. Error"]
  )
}
round(results, 3)
##         cor_p_error cor_p_pstar OLS_beta OLS_bias copula_beta copula_bias
## Normal        0.683       0.999    1.924    0.924       2.091       1.091
## t             0.669       0.970    1.905    0.905       1.155       0.155
## Gamma         0.650       0.948    1.880    0.880       1.048       0.048
## Uniform       0.662       0.977    1.895    0.895       0.862      -0.138
##         copula_SE
## Normal      0.873
## t           0.129
## Gamma       0.098
## Uniform     0.145

Answer:
In all 4 cases, the endogenous regressor is strongly positively correlated with the error term, causing OLS to substantially overestimate the true coefficient of 1. When the endogenous regressor is normally distributed, p and p_star are almost correlated, with a correlation of 0.999. Consequently, the Caussian copula estimator performs poorly and has a large standard error.

In contrast, when p follows a non-normal distribution, the copula correlation sustantially reduces the endogeneity bias. The corrected estimates are 1.155 for the t-distribution, 1.048 for the Gamma distribution, and 0.862 for the uniform distribution. The results show that the method is reasonably robust across different non-normal distributions but does not work well when the endogenous regressor is normally distributed.

4. Implement the refined 2sCOPE approach (Yang et al.).

set.seed(123)
N = 1000
# Generate three independent normal variables
w0 = rnorm(N)
z = rnorm(N)
v = rnorm(N)
# Exogenous variable: correlated with p but not with e
w = qgamma(pnorm(w0), shape = 2, rate = 1)
# Endogenous regressor correlated with both w and e
p_latent = 0.5 * w0 + sqrt(1 - 0.5^2) * z
p = qt(pnorm(p_latent), df = 4)
# Structural error
e = z + v
# True coefficient of p is 1
y = 1 + p + 0.5 * w + e

# Copula transformations
p_star = qnorm(rank(p) / (N + 1))
w_star = qnorm(rank(w) / (N + 1))
# Stage 1
first_stage = lm(p_star ~ w_star)
epsilon_hat = residuals(first_stage)
# Stage 2
model_2sCOPE = lm(y ~ p + w + epsilon_hat)
summary(model_2sCOPE)
## 
## Call:
## lm(formula = y ~ p + w + epsilon_hat)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.8359 -0.6338 -0.0379  0.6569  3.4431 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.07383    0.08826  12.167   <2e-16 ***
## p            1.13927    0.07011  16.250   <2e-16 ***
## w            0.47002    0.04288  10.962   <2e-16 ***
## epsilon_hat  1.03867    0.10762   9.652   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9806 on 996 degrees of freedom
## Multiple R-squared:  0.8876, Adjusted R-squared:  0.8872 
## F-statistic:  2621 on 3 and 996 DF,  p-value: < 2.2e-16
# compare the estimates
comparison_2sCOPE = c(
  true_beta = 1,
  OLS_beta = coef(ols_model)["p"],
  two_sCOPE_beta = coef(model_2sCOPE)["p"]
)
round(comparison_2sCOPE, 3)
##        true_beta       OLS_beta.p two_sCOPE_beta.p 
##            1.000            1.895            1.139

Answer:
OLS estimates the coefficient as 1.895, which is substantially above the true value of 1. The 2sCOPE estimate is 1.139 and is much closer to the true value. Therefore, 2sCOPE removes most of the endogeneity bias.

5. Implement a simulation with a DGP that safely breaks the identifying assumptions.

# Both p and w are normally distributed
w_bad = w0
p_bad = 0.5 * w0 + sqrt(1 - 0.5^2) * z
# Structural error and outcome
e_bad = z + v
y_bad = 1 + p_bad + 0.5 * w_bad + e_bad
# Copula transformations
p_star_bad = qnorm(rank(p_bad) / (N + 1))
w_star_bad = qnorm(rank(w_bad) / (N + 1))
# Stage 1
first_stage_bad = lm(p_star_bad ~ w_star_bad)
epsilon_hat_bad = residuals(first_stage_bad)
# Stage 2
model_2sCOPE_bad = lm(
  y_bad ~ p_bad + w_bad + epsilon_hat_bad)
# Compare with OLS
ols_bad = lm(y_bad ~ p_bad + w_bad)
bad_result = c(
  true_beta = 1,
  OLS_beta = coef(ols_bad)["p_bad"],
  two_sCOPE_beta = coef(model_2sCOPE_bad)["p_bad"],
  two_sCOPE_SE =
    coef(summary(model_2sCOPE_bad))["p_bad", "Std. Error"]
) 
round(bad_result, 3)
##            true_beta       OLS_beta.p_bad two_sCOPE_beta.p_bad 
##                1.000                2.186                3.344 
##         two_sCOPE_SE 
##                0.830

Answer:
When both p and w are normally distributed, the identifying condition of 2sCOPE is violated. The 2sCOPE estimate is far from the true value and has a large standard error. Therefore, the method fails when there is insufficient non-normality for identification.