# Load and prepare Card (1993) NLS Young Men proximity data (nls.dat).
library(tidyverse)
library(here)

card93_region66 <- paste0("reg66", 2:9)

card93_fb_ed <- c("daded", "momed", "nodaded", "nomomed")
card93_fb_famed <- paste0("f", 1:8)
card93_fb_struct <- c("momdad14", "sinmom14")
card93_fb_full <- c(card93_fb_ed, card93_fb_famed, card93_fb_struct)

card93_geo66 <- c(card93_region66, "smsa66r")
card93_demo76 <- c("black", "reg76r", "smsa76r")
card93_exp <- c("exp76", "exp2")

# Exogenous controls (no experience – exp enters separately in wage equations).
card93_x_col1 <- card93_demo76
card93_x_base <- c(card93_demo76, card93_geo66)

# Table II RHS (OLS): education + experience + exogenous controls.
card93_t2_col1 <- c("ed76", card93_exp, card93_x_col1)
card93_t2_base <- c("ed76", card93_exp, card93_x_base)

load_card93_data <- function(
    path = here("04-topics/rep-card1993/Rcode/proximity/nls.dat")) {
  col_names <- c(
    "id", "nearc2", "nearc4", "nearc4a", "nearc4b", "ed76", "ed66", "age76",
    "daded", "nodaded", "momed", "nomomed", "weight", "momdad14", "sinmom14",
    "step14", "reg661", "reg662", "reg663", "reg664", "reg665", "reg666",
    "reg667", "reg668", "reg669", "south66", "work76", "work78", "lwage76",
    "lwage78", "famed", "black", "smsa76r", "smsa78r", "reg76r", "reg78r",
    "reg80r", "smsa66r", "wage76", "wage78", "wage80", "noint78", "noint80",
    "enroll76", "enroll78", "enroll80", "kww", "iq", "marsta76", "marsta78",
    "marsta80", "libcrd14"
  )
  starts <- c(
    1, 7, 10, 12, 15, 18, 21, 24, 27, 33, 35, 41, 43, 56, 58, 60, 62,
    64, 66, 68, 70, 72, 74, 76, 78, 80, 82, 84, 86, 99, 112, 114, 116,
    118, 120, 122, 124, 126, 128, 134, 140, 146, 148, 150, 152, 154,
    156, 159, 163, 165, 167, 169
  )
  ends <- c(
    5, 7, 10, 13, 16, 19, 22, 25, 31, 33, 39, 41, 54, 56, 58, 60, 62,
    64, 66, 68, 70, 72, 74, 76, 78, 80, 82, 84, 97, 110, 112, 114, 116,
    118, 120, 122, 124, 126, 132, 138, 144, 146, 148, 150, 152, 154,
    157, 161, 163, 165, 167, 169
  )
  col_spec <- fwf_positions(start = starts, end = ends, col_names = col_names)

  raw <- read_fwf(path, col_spec, na = c(".", ""))

  raw |>
    mutate(
      across(
        c(
          "nearc2", "nearc4", "nearc4a", "nearc4b", "nodaded", "nomomed",
          "momdad14", "sinmom14", "step14", "reg661", "reg662", "reg663",
          "reg664", "reg665", "reg666", "reg667", "reg668", "reg669",
          "south66", "work76", "work78", "black", "smsa76r", "smsa78r",
          "reg76r", "reg78r", "reg80r", "smsa66r", "noint78", "noint80",
          "enroll76", "enroll78", "enroll80", "libcrd14"
        ),
        as.integer
      ),
      exp76 = age76 - ed76 - 6,
      exp2 = exp76^2 / 100,
      age2 = age76^2 / 100,
      age66 = age76 - 10,
      f1 = as.integer(famed == 1),
      f2 = as.integer(famed == 2),
      f3 = as.integer(famed == 3),
      f4 = as.integer(famed == 4),
      f5 = as.integer(famed == 5),
      f6 = as.integer(famed == 6),
      f7 = as.integer(famed == 7),
      f8 = as.integer(famed == 8),
      # Low family background: neither parent completed high school (famed class 9).
      lowfam = as.integer(famed == 9),
      nearc4_low = nearc4 * lowfam
    )
}

card93_wage_sample <- function(data) {
  data |> filter(!is.na(lwage76))
}

card93_coef_se <- function(model, term) {
  s <- summary(model)$coefficients
  if (!term %in% rownames(s)) {
    return(list(estimate = NA_real_, std.error = NA_real_))
  }
  list(estimate = s[term, 1], std.error = s[term, 2])
}

card93_ols <- function(data, y, xvars) {
  fml <- as.formula(paste(y, "~", paste(xvars, collapse = " + ")))
  lm(fml, data = data)
}

card93_iv_panel_b <- function(data, y, endog, instruments, controls) {
  # Panel B: education and experience endogenous; age/age2 + nearc4 as IVs.
  exog <- unique(controls)
  inst <- unique(c(instruments, "age76", "age2", exog))
  endog_vars <- unique(c(endog, "exp76", "exp2"))
  fml <- as.formula(
    paste(
      y, "~", paste(c(endog_vars, exog), collapse = " + "),
      "|", paste(inst, collapse = " + ")
    )
  )
  AER::ivreg(fml, data = data)
}

# Panel B reduced form: instrument experience with age/age2; nearc4 on RHS.
card93_rf_panel_b <- function(data, y, controls) {
  exog <- unique(controls)
  fml <- as.formula(
    paste(
      y, "~ nearc4 + exp76 + exp2 +", paste(exog, collapse = " + "),
      "| nearc4 + age76 + age2 +", paste(exog, collapse = " + ")
    )
  )
  AER::ivreg(fml, data = data)
}