#==== 子集内生性检验代码 ====

#==== 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)





# 根据理论分析，我们需要将变量分为三组：
# X1: 外生变量 (black, south, urban)
# X2: 待检验的潜在内生变量 (edu) 
# X3: 已知内生变量 (exp, exp2)
# 工具变量Z: (public, private, age, age2) + X1

# 定义变量组
X1_vars <- c("black", "south", "urban")  # 外生协变量
X2_vars <- c("edu")  # 待检验的潜在内生变量
X3_vars <- c("exp", "exp2")  # 已知内生变量
Z_2 <- c("public", "private")  # 工具变量2
Z_3 <- c("age", "age2")  # 工具变量3
Z_vars <- c(Z_2, Z_3)  # 工具变量


# 构建模型公式
# 完整模型：Y = X1'β1 + X2'β2 + X3'β3 + e
# 其中 X2 = edu, X3 = (exp, exp2)

# 模型0：基准模型，使用工具变量 (public, private, age, age2) 的2SLS
# 把edu、exp、exp2作为内生变量，测试三者的内生性
formula_baseline <- as.formula(paste("lwage ~", 
                                paste(c(X1_vars, X2_vars, X3_vars), collapse = " + "),
                                "|", 
                                paste(c(X1_vars, Z_vars ), collapse = " + ")))

# 估计模型0：使用工具变量 (Z, X1) 的2SLS，X2和X3作为内生变量
# 即使用 (public, private, age, age2, black, south, urban, edu) 作为工具变量
model0 <- ivreg(formula_baseline, data = tbl_reg)
summary_model0 <- summary(model0, diagnostics = TRUE)
# 提取Wu-Hausman F检验结果
F_stat_model0 <- summary_model0$diagnostics["Wu-Hausman","statistic"]
F_pvalue_model0 <- summary_model0$diagnostics["Wu-Hausman","p-value"]


# 模型1: 使用工具变量 Z=(X1, X2, Z2, Z3) 的2SLS估计
# 把edu作为外生变量（H0: edu是外生变量）
# 把exp、exp2作为内生变量，测试二者的内生性
# 即使用 ( black, south, urban, edu,age, age2, public, private) 作为工具变量
formula_H0 <- as.formula(paste("lwage ~", 
                                paste(c(X1_vars, X2_vars, X3_vars), collapse = " + "),
                                "|", 
                                paste(c(X1_vars, X2_vars, Z_2, Z_3), collapse = " + ")))


# 估计模型1：使用工具变量 (X1, X2, Z2, Z3 ) 的2SLS
# cat("模型1: 2SLS估计，内生变量包含edu\n")
model1 <- ivreg(formula_H0, data = tbl_reg)
summary_model1 <- summary(model1, diagnostics = TRUE)
# 提取Wu-Hausman F检验结果
F_stat_model1 <- summary_model1$diagnostics["Wu-Hausman","statistic"]
F_pvalue_model1 <- summary_model1$diagnostics["Wu-Hausman","p-value"]


# 模型2: 使用工具变量Z= (X1, Z2, Z3) 的2SLS估计
# 把exp、exp2作为内生变量
# 即使用 ( black, south, urban, age, age2, public, private) 作为工具变量
formula_H1 <- as.formula(paste("lwage ~", 
                                      paste(c(X1_vars, X2_vars, X3_vars), collapse = " + "),
                                      "|", 
                                      paste(c(X1_vars, Z_2, Z_3), collapse = " + ")))


# 估计模型2：使用工具变量 Z 的2SLS
# cat("\n模型2: 2SLS估计，内生变量不包含edu\n")
model2 <- ivreg(formula_H1, data = tbl_reg)
summary_model2 <- summary(model2, diagnostics = TRUE)
# 提取Wu-Hausman F检验结果
F_stat_model2 <- summary_model2$diagnostics["Wu-Hausman","statistic"]
F_pvalue_model2 <- summary_model2$diagnostics["Wu-Hausman","p-value"]

# =====进行子集内生性检验=====

# 提取edu的系数估计值 β̂2
beta2_hat <- coef(model1)["edu"]
var_beta2_hat <- vcov(model1)["edu", "edu"]

# cat("\n模型1中edu的系数估计值 β̂2 =", round(beta2_hat, 6), "\n")
# cat("模型1中edu的方差估计值 =", round(var_beta2_hat, 8), "\n")

# 提取edu的系数估计值 β̃2
beta2_tilde <- coef(model2)["edu"]
var_beta2_tilde <- vcov(model2)["edu", "edu"]

# cat("\n模型2中edu的系数估计值 β̃2 =", round(beta2_tilde, 6), "\n")
# cat("模型2中edu的方差估计值 =", round(var_beta2_tilde, 8), "\n")

# 计算Durbin-Wu-Hausman检验统计量
# T = (β̂2 - β̃2)' * (Var(β̃2) - Var(β̂2))^(-1) * (β̂2 - β̃2)

diff_beta <- beta2_hat - beta2_tilde
diff_var <- var_beta2_tilde - var_beta2_hat

# 检验统计量，T统计量将服从卡方分布
T_stat <- (diff_beta^2) / diff_var
# 自由度
df <- 1  # edu是单个变量
# p值
p_value <- 1 - pchisq(T_stat, df)

# 临界值
alpha_005 <- qchisq(0.95, df)
alpha_001 <- qchisq(0.99, df)


# F版本检验
## 注意：因为即使在H0假设下（edu是外生变量），模型中还存在其他内生自变量exp、exp2
## 因此，我们不能简单地构造使用F形式的统计量，而只能使用T统计量
F_stat <- T_stat / df
F_p_value <- 1 - pf(F_stat, df, model1$df.residual)