#==== data prepare ====
require(here)
require(magrittr)
require(tidyverse)
require(broom)
require(haven)
require(systemfit)
require(AER)
#renv::install("AER")
require(ivmodel) # for LIML fit

# 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 |>
  select(all_of(vars_list))

# ==== 模型1: 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)


# 过度识别检验 - Sargan统计量
# 方法1: 使用summary函数的diagnostics参数
smry_tsls_pp_diag <- summary(fit_tsls_pp, diagnostics = TRUE)
overid_test_pp <- smry_tsls_pp_diag$diagnostics
sargan_value_pp <- overid_test_pp["Sargan","statistic"]
sargan_pvalue_pp <- overid_test_pp["Sargan","p-value"]

# 手动计算Sargan统计量（验证）
# 获取2SLS残差
residuals_pp <- residuals(fit_tsls_pp)
n <- length(residuals_pp)

# 获取工具变量矩阵
Z_pp <- model.matrix(~ public + private + exp + exp2 + black + south + urban, data = tbl_reg)
# 获取内生变量矩阵
X_pp <- model.matrix(~ edu + exp + exp2 + black + south + urban, data = tbl_reg)

# 计算Sargan统计量
sigma2_hat <- sum(residuals_pp^2) / n
ZtZ_inv <- solve(t(Z_pp) %*% Z_pp)
sargan_manual_pp <- (t(residuals_pp) %*% Z_pp %*% ZtZ_inv %*% t(Z_pp) %*% residuals_pp) / sigma2_hat
sargan_manual_pvalue_pp <- 1 - pchisq(as.numeric(sargan_manual_pp), df = ncol(Z_pp) - ncol(X_pp))

# ==== 模型2: 2SLS估计 - 内生变量edu、exp、exp2，工具变量public、private、age、age2 ====

# 2SLS模型设定
mod_tsls_ppa <- formula("lwage ~ edu + exp + exp2 + black + south + urban | public + private + age + age2 + black + south + urban")
fit_tsls_ppa <- AER::ivreg(formula = mod_tsls_ppa, data = tbl_reg)
smry_tsls_ppa <- summary(fit_tsls_ppa)


# 过度识别检验 - Sargan统计量
# 方法1: 使用summary函数的diagnostics参数
smry_tsls_ppa_diag <- summary(fit_tsls_ppa, diagnostics = TRUE)
overid_test_ppa <- smry_tsls_ppa_diag$diagnostics
sargan_value_ppa <- overid_test_ppa["Sargan","statistic"]
sargan_pvalue_ppa <- overid_test_ppa["Sargan","p-value"]


# 手动计算Sargan统计量（验证）
# 获取2SLS残差
residuals_ppa <- residuals(fit_tsls_ppa)

# 获取工具变量矩阵
Z_ppa <- model.matrix(~ public + private + age + age2 + black + south + urban, data = tbl_reg)
# 获取内生变量矩阵
X_ppa <- model.matrix(~ edu + exp + exp2 + black + south + urban, data = tbl_reg)

# 计算Sargan统计量
sigma2_hat_ppa <- sum(residuals_ppa^2) / n
ZtZ_inv_ppa <- solve(t(Z_ppa) %*% Z_ppa)
sargan_manual_ppa <- (t(residuals_ppa) %*% Z_ppa %*% ZtZ_inv_ppa %*% t(Z_ppa) %*% residuals_ppa) / sigma2_hat_ppa
sargan_manual_pvalue_ppa <- 1 - pchisq(as.numeric(sargan_manual_ppa), df = ncol(Z_ppa) - ncol(X_ppa))