#==== data prepare ==== require(here) require(magrittr) require(tidyverse) require(broom) require(haven) require(systemfit) require(AER) #renv::install("AER") require(ivmodel) # for LIML fit require(MASS) # for ginv function # read data file_path <- here::here("hansen-programs/Card1995.dta") card1995 <- haven::read_dta(file_path) tbl_tidy <- card1995 |> mutate(lwage = lwage76, edu = ed76, exp = age76 - ed76 - 6, exp2 = (exp^2)/100, age = age76, age2 = (age^2)/100, south = reg76r, urban = smsa76r, college = nearc4, public = nearc4a, private = nearc4b ) |> filter(!is.na(lwage)) %>% rename_at(all_of(c("daded","momed")), ~ all_of(c("dadedu","momedu"))) vars_list <- c("lwage", "edu", "exp", "exp2", "age", "age2", "black","south", "urban", "college","public","private","momedu","dadedu") tbl_reg <- tbl_tidy |> dplyr::select(!!vars_list) # ==== 模型1: OLS估计 - 作为对比基准 ==== # OLS模型设定（将edu视为外生变量） mod_ols <- formula("lwage ~ edu + exp + exp2 + black + south + urban") fit_ols <- lm(formula = mod_ols, data = tbl_reg) smry_ols <- summary(fit_ols) # ==== 模型2: 2SLS估计 - 内生变量edu，工具变量public、private ==== # 2SLS模型设定 mod_tsls_pp <- formula("lwage ~ edu + exp + exp2 + black + south + urban | public + private + exp + exp2 + black + south + urban") fit_tsls_pp <- AER::ivreg(formula = mod_tsls_pp, data = tbl_reg) smry_tsls_pp <- summary(fit_tsls_pp) # ==== 内生性检验方法1: Hausman检验 ==== # 提取OLS和2SLS的系数估计 ols_coef <- coef(fit_ols) tsls_coef <- coef(fit_tsls_pp) # 计算系数差异（只比较edu的系数） edu_diff <- ols_coef["edu"] - tsls_coef["edu"] # 手动计算Hausman检验统计量 # 根据Hansen教材公式(12.64)，Hausman检验统计量为： # T = (beta_2sls - beta_ols)' * (Var(beta_2sls) - Var(beta_ols))^(-1) * (beta_2sls - beta_ols) / sigma^2 # 提取系数 ols_coef_edu <- ols_coef["edu"] tsls_coef_edu <- tsls_coef["edu"] # 获取协方差矩阵 ols_var_edu <- vcov(fit_ols)["edu", "edu"] tsls_var_edu <- vcov(fit_tsls_pp)["edu", "edu"] # 计算方差差异 var_diff <- tsls_var_edu - ols_var_edu # 使用2SLS回归的残差方差作为sigma^2的估计 sigma2_hat <- summary(fit_tsls_pp)$sigma^2 # 计算Hausman检验统计量 # 注意：这里需要确保方差差异为正数（这是Hausman检验的前提条件） if (var_diff > 0) { hausman_T_stat <- (tsls_coef_edu - ols_coef_edu)^2 / (var_diff / sigma2_hat) hausman_T_pvalue <- 1 - pchisq(hausman_T_stat, df = 1) } else { # 如果方差差异为负，说明可能违反了Hausman检验的前提条件 hausman_T_stat <- NA hausman_T_pvalue <- NA } # 更准确的Hausman检验：使用矩阵形式 # 根据Hansen教材，当k2=1时，可以使用更简单的形式 # 但为了完整性，我们也提供矩阵形式的计算 # 提取完整的系数向量（只针对edu变量） beta_diff <- tsls_coef_edu - ols_coef_edu # 计算方差差异矩阵（对于单个变量，这是一个标量） var_diff_matrix <- tsls_var_edu - ols_var_edu # 使用更标准的Hausman检验公式 if (var_diff_matrix > 0) { # 标准Hausman检验统计量（不除以sigma^2） hausman_stat_standard <- beta_diff^2 / var_diff_matrix # 使用卡方分布计算p值 hausman_pvalue_standard <- 1 - pchisq(hausman_stat_standard, df = 1) # 使用F分布（有限样本调整） n <- nrow(tbl_reg) k <- length(coef(fit_ols)) # 参数个数 hausman_F_stat <- hausman_stat_standard hausman_F_pvalue <- 1 - pf(hausman_F_stat, df1 = 1, df2 = n - k) } else { hausman_stat_standard <- NA hausman_pvalue_standard <- NA hausman_F_stat <- NA hausman_F_pvalue <- NA } # ==== 内生性检验方法2: 控制函数回归方法（手动进行wald F检验） ==== # 第一步：估计简化形式（edu对工具变量的回归） reduced_form <- lm(edu ~ public + private + exp + exp2 + black + south + urban, data = tbl_reg) # 提取残差 edu_residuals <- residuals(reduced_form) # 第二步：控制函数回归 control_function_data <- tbl_reg %>% mutate(edu_resid = edu_residuals) control_function <- lm(lwage ~ edu + exp + exp2 + black + south + urban + edu_resid, data = control_function_data) # 控制函数回归结果 smry_control_function <- summary(control_function) # 控制函数回归的残差方差 sigma2_control_function <- smry_control_function$sigma^2 # 检验1方法1：检验edu_resid的系数是否显著（内生性检验）：等价于t检验（仅有一个内生变量时） edu_resid_coef <- smry_control_function$coefficients["edu_resid","Estimate"] edu_resid_se <- smry_control_function$coefficients["edu_resid", "Std. Error"] edu_resid_t <- edu_resid_coef / edu_resid_se edu_resid_t_squared <- edu_resid_t^2 # 等价于Hausman T统计量 edu_resid_p <- 2 * (1 - pt(abs(edu_resid_t), df = df.residual(control_function))) # 检验方法2：使用lmtest包进行受约束检验 require(lmtest) # 对control_function模型进行受约束检验，约束edu_resid的系数为0 # 这等价于检验原假设H0: α = 0（edu是外生的） waldtest_result <- waldtest(control_function, . ~ . - edu_resid, # 约束edu_resid系数为0 test = "F") # 提取F统计量和p值 wald_F_stat <- waldtest_result$F[2] wald_F_pvalue <- waldtest_result$`Pr(>F)`[2] # 也可以使用Chisq检验 waldtest_chisq <- waldtest(control_function, . ~ . - edu_resid, # 约束edu_resid系数为0 test = "Chisq") wald_chisq_stat <- waldtest_chisq$Chisq[2] wald_chisq_pvalue <- waldtest_chisq$`Pr(>Chisq)`[2] # 方法3：使用linearHypothesis包进行受约束检验 # 也可以使用linearHypothesis包进行受约束F检验 linearHypothesis_result <- linearHypothesis(control_function, "edu_resid = 0", # 约束edu_resid系数为0 test = "F") F_stat_linearHypothesis <- linearHypothesis_result$F[2] F_pvalue_linearHypothesis <- linearHypothesis_result$`Pr(>F)`[2] # 也可以使用linearHypothesis包进行受约束卡方检验 linearHypothesis_result2 <- linearHypothesis(control_function, "edu_resid = 0", # 约束edu_resid系数为0 test = "Chisq") Chisq_stat_linearHypothesis <- linearHypothesis_result2$Chisq[2] Chisq_pvalue_linearHypothesis <- linearHypothesis_result2$`Pr(>Chisq)`[2] # ==== 内生性检验方法3: 使用summary函数进行的endogeneity检验 (获得wald F检验结果)==== endogeneity_test <- summary(fit_tsls_pp, diagnostics = TRUE) endogeneity_test_pp <- endogeneity_test$diagnostics WH_F_value_pp <- endogeneity_test_pp["Wu-Hausman","statistic"] WH_F_pvalue_pp <- endogeneity_test_pp["Wu-Hausman","p-value"] # ==== 使用矩阵方法计算Hausman T统计量并验证与Wald W0的等价性 ==== # 重新分析Hausman T统计量的计算 # 根据Hansen教材，当k2=1时，Hausman检验统计量应该等于控制函数回归中残差项系数的t统计量平方 # 方法1：直接使用方差差异（不除以sigma^2） # 根据Hansen教材，标准Hausman检验统计量是： # T = (beta_2sls - beta_ols)^2 / (Var(beta_2sls) - Var(beta_ols)) # 计算系数差异 beta_diff <- tsls_coef_edu - ols_coef_edu # 获取OLS和2SLS的方差 ols_var_edu <- vcov(fit_ols)["edu", "edu"] tsls_var_edu <- vcov(fit_tsls_pp)["edu", "edu"] # 计算方差差异 var_diff_direct <- tsls_var_edu - ols_var_edu # 考虑使用2SLS回归的残差方差作为sigma^2的估计 sigma2_hat <- summary(fit_tsls_pp)$sigma^2 # 也可以考虑使用控制函数回归的残差方差作为sigma^2的估计 sigma2_hat <- summary(control_function)$sigma^2 # 标准Hausman检验统计量（不除以sigma^2） if (var_diff_direct > 0) { hausman_T_standard <- beta_diff^2/var_diff_direct hausman_T_standard_pvalue <- 1 - pchisq(hausman_T_standard, df = 1) # cat("标准Hausman T统计量:", hausman_T_standard, "\n") # cat("标准Hausman T检验p值:", hausman_T_standard_pvalue, "\n") } else { hausman_T_standard <- NA hausman_T_standard_pvalue <- NA # cat("方差差异为负，Hausman检验前提条件不满足\n") } # 方法2：使用矩阵方法重新计算 # 构建设计矩阵 # 重要：model.matrix(fit_ols)已经包含了截距项（第一列） # 在Hansen教材的矩阵公式中，X1应该包含所有外生变量，包括截距项 # 这是因为截距项在回归模型中是一个重要的外生变量 X_full <- model.matrix(fit_ols) # 完整设计矩阵，包含截距项 X1 <- X_full[, c("(Intercept)", "exp", "exp2", "black", "south", "urban")] # 外生变量（包含截距） X2 <- X_full[, "edu"] |> as.matrix() # 内生变量 # 工具变量矩阵也需要包含截距项 # 在2SLS中，外生变量（包括截距）作为自己的工具变量 Z1 <- X1 # 外生变量作为自己的工具变量（包含截距） Z2 <- as.matrix(tbl_reg[, c("public", "private")]) # 工具变量 Z <- cbind(Z1, Z2) # 完整工具变量矩阵（包含截距） n <- nrow(X1) # 构建投影矩阵 P_Z <- Z %*% solve(t(Z) %*% Z) %*% t(Z) # P_Z = Z(Z'Z)^(-1)Z' P_1 <- X1 %*% solve(t(X1) %*% X1) %*% t(X1) # P_1 = X1(X1'X1)^(-1)X1' M_1 <- diag(n) - P_1 # M_1 = I_n - P_1 # 计算各项 X2_PZ_P1_X2 <- t(X2) %*% (P_Z - P_1) %*% X2 X2_M1_X2 <- t(X2) %*% M_1 %*% X2 # cat("X2'(P_Z - P_1)X2:", X2_PZ_P1_X2, "\n") # cat("X2'M_1X2:", X2_M1_X2, "\n") # 考虑使用2SLS回归的残差方差作为sigma^2的估计 sigma2_hat <- summary(fit_tsls_pp)$sigma^2 # 也可以考虑使用控制函数回归的残差方差作为sigma^2的估计 sigma2_hat <- summary(control_function)$sigma^2 # 计算逆矩阵的差异（根据Hansen公式） var_diff_matrix <- solve(X2_PZ_P1_X2) - solve(X2_M1_X2) # 再次取逆矩阵 var_diff_matrix_inv <- solve(var_diff_matrix) # 计算Hausman T统计量（矩阵形式，除以sigma^2） hausman_T_matrix <- beta_diff^2*var_diff_matrix_inv/ sigma2_hat hausman_T_matrix_pvalue <- 1 - pchisq(hausman_T_matrix, df = 1) # ==== 验证与Wald W0统计量的等价性 ==== # 根据Hansen教材，Hausman T统计量应该等于控制函数回归中残差项系数的Wald统计量 # 即：Hausman T = (alpha_hat / se(alpha_hat))^2 # 验证等价性（允许小的数值误差） diff_standard <- abs(edu_resid_t_squared - hausman_T_standard) diff_matrix <- abs(edu_resid_t_squared - hausman_T_matrix)