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

# 获取诊断信息
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统计量（验证）
residuals_ppa <- residuals(fit_tsls_ppa)
n <- length(residuals_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)
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))

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

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

# 获取诊断信息
smry_tsls_ppdma_diag <- summary(fit_tsls_ppdma, diagnostics = TRUE)
overid_test_ppdma <- smry_tsls_ppdma_diag$diagnostics
sargan_value_ppdma <- overid_test_ppdma["Sargan","statistic"]
sargan_pvalue_ppdma <- overid_test_ppdma["Sargan","p-value"]

# 手动计算Sargan统计量（验证）
residuals_ppdma <- residuals(fit_tsls_ppdma)
Z_ppdma <- model.matrix(~ public + private + dadedu + momedu + age + age2 + black + south + urban, data = tbl_reg)
X_ppdma <- model.matrix(~ edu + exp + exp2 + black + south + urban, data = tbl_reg)
sigma2_hat_ppdma <- sum(residuals_ppdma^2) / n
ZtZ_inv_ppdma <- solve(t(Z_ppdma) %*% Z_ppdma)
sargan_manual_ppdma <- (t(residuals_ppdma) %*% Z_ppdma %*% ZtZ_inv_ppdma %*% t(Z_ppdma) %*% residuals_ppdma) / sigma2_hat_ppdma
sargan_manual_pvalue_ppdma <- 1 - pchisq(as.numeric(sargan_manual_ppdma), df = ncol(Z_ppdma) - ncol(X_ppdma))

# ==== 子集过度识别检验：N统计量和C统计量 ====

# 根据Hansen 12.32节，我们需要将工具变量分为两部分：
# Za: 被认为与误差项不相关的工具变量 (public, private, age, age2)
# Zb: 可能与误差项相关的工具变量 (dadedu, momedu)

# 模型A: 仅使用Za作为工具变量
mod_tsls_za <- formula("lwage ~ edu + exp + exp2 + black + south + urban | public + private + age + age2 + black + south + urban")
fit_tsls_za <- AER::ivreg(formula = mod_tsls_za, data = tbl_reg)
smry_tsls_za <- summary(fit_tsls_za)

# 获取模型A的诊断信息
smry_tsls_za_diag <- summary(fit_tsls_za, diagnostics = TRUE)
overid_test_za <- smry_tsls_za_diag$diagnostics
sargan_value_za <- overid_test_za["Sargan","statistic"]
sargan_pvalue_za <- overid_test_za["Sargan","p-value"]

# 手动计算模型A的Sargan统计量
residuals_za <- residuals(fit_tsls_za)
Z_za <- model.matrix(~ public + private + age + age2 + black + south + urban, data = tbl_reg)
X_za <- model.matrix(~ edu + exp + exp2 + black + south + urban, data = tbl_reg)
sigma2_hat_za <- sum(residuals_za^2) / n
ZtZ_inv_za <- solve(t(Z_za) %*% Z_za)
sargan_manual_za <- (t(residuals_za) %*% Z_za %*% ZtZ_inv_za %*% t(Z_za) %*% residuals_za) / sigma2_hat_za
sargan_manual_pvalue_za <- 1 - pchisq(as.numeric(sargan_manual_za), df = ncol(Z_za) - ncol(X_za))

# 计算C统计量：C = S - Sa
# S: 使用全部工具变量的Sargan统计量 (模型2)
# Sa: 仅使用Za的Sargan统计量 (模型A)
C_statistic <- as.numeric(sargan_manual_ppdma) - as.numeric(sargan_manual_za)
C_pvalue <- 1 - pchisq(C_statistic, df = 2)  # 自由度 = Zb的维度 = 2 (dadedu, momedu)

# 计算N统计量 (Newey统计量)
# 根据Hansen 12.32节公式(12.72)
# N = e'R(R'R - R'X(X'X)^(-1)X'R)^(-1)R'e / sigma^2
# 其中 R = Ma * Zb, Ma = I - Za(Za'Za)^(-1)Za'

# 构建Ma矩阵
Ma <- diag(n) - Z_za %*% solve(t(Z_za) %*% Z_za) %*% t(Z_za)

# 构建Zb矩阵 (dadedu, momedu)
Zb <- model.matrix(~ dadedu + momedu, data = tbl_reg)

# 计算R = Ma * Zb
R <- Ma %*% Zb

# 计算X_hat = P * X, 其中P = Z(Z'Z)^(-1)Z'
P <- Z_ppdma %*% solve(t(Z_ppdma) %*% Z_ppdma) %*% t(Z_ppdma)
X_hat <- P %*% X_ppdma

# 计算N统计量 - 使用更稳定的数值方法
RTR <- t(R) %*% R
RTX <- t(R) %*% X_hat
XTR <- t(X_hat) %*% R
XTX <- t(X_hat) %*% X_hat

# 使用广义逆来避免奇异性问题
# 计算 (R'R - R'X(X'X)^(-1)X'R)^(-1)
middle_term <- RTR - RTX %*% MASS::ginv(XTX) %*% XTR

# 检查middle_term的条件数
cond_num <- rcond(middle_term)
if(cond_num < 1e-10) {
  #cat("警告: middle_term矩阵条件数很小，使用广义逆\n")
  middle_term_inv <- MASS::ginv(middle_term)
} else {
  middle_term_inv <- solve(middle_term)
}

N_statistic <- (t(residuals_ppdma) %*% R %*% middle_term_inv %*% t(R) %*% residuals_ppdma) / sigma2_hat_ppdma
N_pvalue <- 1 - pchisq(as.numeric(N_statistic), df = 2)  # 自由度 = Zb的维度 = 2

# 验证N统计量和C统计量的等价性

diff_N_C <- abs(as.numeric(N_statistic) - C_statistic)