Untitled

############################################
# TABLE 2 recreation (Cols 1-2) from a fresh script
# Assumes Table 1B script has ALREADY been run and these objects exist:
# - t1b_panelB   (final Panel B sample)
# - annual_close (optional, not required here)
#
# This recreates the version you ended up with:
# - Table 2 sample nested inside final Table 1B sample
# - 1976-2020
# - Firm FE + Year FE
# - State-clustered SEs
# - Col (1): baseline controls
# - Col (2): + Herfindahl
# - 4 decimal display
# - Professional gt export table
############################################

# -----------------------------
# Package checks
# -----------------------------
stopifnot(requireNamespace("dplyr", quietly = TRUE))
stopifnot(requireNamespace("fixest", quietly = TRUE))
stopifnot(requireNamespace("modelsummary", quietly = TRUE))
stopifnot(requireNamespace("tibble", quietly = TRUE))
stopifnot(requireNamespace("gt", quietly = TRUE))
stopifnot(requireNamespace("broom", quietly = TRUE))

library(dplyr)
library(fixest)
library(modelsummary)
library(tibble)
library(gt)
library(broom)

# -----------------------------
# Winsorization helper (paper text)
# -----------------------------
winsor_1_99 <- function(x) {
  if (!is.numeric(x)) return(x)
  q <- quantile(x, c(0.01, 0.99), na.rm = TRUE, names = FALSE)
  pmin(pmax(x, q[1]), q[2])
}

# -----------------------------
# 1) Build Table 2 sample from final Table 1B sample
#    (ensures Table 2 N <= Table 1B N)
# -----------------------------
t2_base <- t1b_panelB %>%
  filter(year >= 1976, year <= 2020) %>%   # Table 2 period
  mutate(
    gvkey = as.factor(gvkey),
    state = as.factor(state),
    year = as.integer(year)
  )

# Column-specific samples
t2_c1 <- t2_base %>%
  filter(
    !is.na(q),
    !is.na(cash_flow),
    !is.na(leverage)
  )

t2_c2 <- t2_base %>%
  filter(
    !is.na(q),
    !is.na(cash_flow),
    !is.na(leverage),
    !is.na(herfindahl)
  )

# Winsorize regression variables (does not change N, only values)
t2_c1 <- t2_c1 %>%
  mutate(
    rd_intensity = winsor_1_99(rd_intensity),
    q = winsor_1_99(q),
    cash_flow = winsor_1_99(cash_flow),
    ln_sales = winsor_1_99(ln_sales),
    ln_age = winsor_1_99(ln_age),
    tangibility = winsor_1_99(tangibility),
    profitability = winsor_1_99(profitability),
    leverage = winsor_1_99(leverage)
  )

t2_c2 <- t2_c2 %>%
  mutate(
    rd_intensity = winsor_1_99(rd_intensity),
    q = winsor_1_99(q),
    cash_flow = winsor_1_99(cash_flow),
    ln_sales = winsor_1_99(ln_sales),
    ln_age = winsor_1_99(ln_age),
    tangibility = winsor_1_99(tangibility),
    profitability = winsor_1_99(profitability),
    leverage = winsor_1_99(leverage),
    herfindahl = winsor_1_99(herfindahl)
  )

cat("N checks:\n")
cat("  Table 1B Panel B N:", nrow(t1b_panelB), "\n")
cat("  Table 2 Col (1) pre-FE:", nrow(t2_c1), "\n")
cat("  Table 2 Col (2) pre-FE:", nrow(t2_c2), "\n")

# -----------------------------
# 2) Estimate Table 2 columns (1) and (2)
#    Equivalent to: R&D Intensity_ijt = alpha_i + gamma_t + beta*CloseElection_jt + X_it + e_ijt
# -----------------------------
m1 <- feols(
  rd_intensity ~ close_election_dummy + q + cash_flow + ln_sales + ln_age +
    tangibility + profitability + leverage | gvkey + year,
  data = t2_c1,
  cluster = ~ state
)

m2 <- feols(
  rd_intensity ~ close_election_dummy + q + cash_flow + ln_sales + ln_age +
    tangibility + profitability + leverage + herfindahl | gvkey + year,
  data = t2_c2,
  cluster = ~ state
)

# -----------------------------
# 3) Paper-like console table (4 decimals)
# -----------------------------
add_rows <- tibble::tribble(
  ~term,      ~`(1)`, ~`(2)`,
  "Firm FE",  "Yes",  "Yes",
  "Year FE",  "Yes",  "Yes"
)

tab2_cols12 <- modelsummary(
  list("(1)" = m1, "(2)" = m2),
  stars = c("*" = .1, "**" = .05, "***" = .01),
  estimate = "{estimate}{stars}",
  statistic = "({std.error})",
  fmt = 4,
  coef_map = c(
    "close_election_dummy" = "Close election dummy",
    "q" = "q",
    "cash_flow" = "Cash flow",
    "ln_sales" = "ln(Sales)",
    "ln_age" = "ln(Age)",
    "tangibility" = "Tangibility",
    "profitability" = "Profitability",
    "leverage" = "Leverage",
    "herfindahl" = "Herfindahl"
  ),
  gof_map = tibble::tribble(
    ~raw,        ~clean, ~fmt,
    "nobs",      "N",    0,
    "r.squared", "R2",   4
  ),
  add_rows = add_rows,
  output = "markdown",
  title = "Table 2. Political uncertainty and R&D intensity: Baseline results (Cols 1-2)"
)

tab2_cols12

# -----------------------------
# 4) Professional gt export (similar style to your Table 1B export)
# -----------------------------
fmt_coef <- function(est, pval) {
  stars <- ifelse(is.na(pval), "",
                  ifelse(pval < 0.01, "***",
                         ifelse(pval < 0.05, "**",
                                ifelse(pval < 0.1, "*", ""))))
  paste0(sprintf("%.4f", est), stars)
}
fmt_se <- function(se) paste0("(", sprintf("%.4f", se), ")")

coef_order <- c("close_election_dummy", "q", "cash_flow", "ln_sales", "ln_age",
                "tangibility", "profitability", "leverage", "herfindahl")

coef_labels <- c(
  close_election_dummy = "Close election dummy",
  q = "q",
  cash_flow = "Cash flow",
  ln_sales = "ln(Sales)",
  ln_age = "ln(Age)",
  tangibility = "Tangibility",
  profitability = "Profitability",
  leverage = "Leverage",
  herfindahl = "Herfindahl"
)

tidy_to_display <- function(mod, colname) {
  broom::tidy(mod) %>%
    filter(term %in% coef_order) %>%
    mutate(
      term = factor(term, levels = coef_order),
      Variable = coef_labels[as.character(term)],
      estimate_fmt = fmt_coef(estimate, p.value),
      se_fmt = fmt_se(std.error)
    ) %>%
    arrange(term) %>%
    select(Variable, estimate_fmt, se_fmt) %>%
    rename(!!paste0(colname, "_est") := estimate_fmt,
           !!paste0(colname, "_se") := se_fmt)
}

d1 <- tidy_to_display(m1, "c1")
d2 <- tidy_to_display(m2, "c2")

rows_template <- tibble::tibble(Variable = unname(coef_labels[coef_order]))

tab_core <- rows_template %>%
  left_join(d1, by = "Variable") %>%
  left_join(d2, by = "Variable") %>%
  mutate(
    c1_est = ifelse(is.na(c1_est), "", c1_est),
    c1_se  = ifelse(is.na(c1_se),  "", c1_se),
    c2_est = ifelse(is.na(c2_est), "", c2_est),
    c2_se  = ifelse(is.na(c2_se),  "", c2_se)
  )

tab_display <- bind_rows(
  lapply(seq_len(nrow(tab_core)), function(i) {
    tibble::tibble(
      Variable = c(tab_core$Variable[i], ""),
      ` (1) ` = c(tab_core$c1_est[i], tab_core$c1_se[i]),
      ` (2) ` = c(tab_core$c2_est[i], tab_core$c2_se[i]),
      row_type = c("coef", "se")
    )
  })
)

# Robust R2 extraction
r2_1 <- tryCatch(broom::glance(m1)$r.squared[1], error = function(e) NA_real_)
r2_2 <- tryCatch(broom::glance(m2)$r.squared[1], error = function(e) NA_real_)

gof_rows <- tibble::tibble(
  Variable = c("N", "R2", "Firm FE", "Year FE"),
  ` (1) ` = c(format(as.numeric(nobs(m1)), big.mark = ","), sprintf("%.4f", r2_1), "Yes", "Yes"),
  ` (2) ` = c(format(as.numeric(nobs(m2)), big.mark = ","), sprintf("%.4f", r2_2), "Yes", "Yes"),
  row_type = c("gof", "gof", "gof", "gof")
)

tab_display <- bind_rows(tab_display, gof_rows)

gt_tab2 <- tab_display %>%
  gt() %>%
  tab_header(
    title = md("**Table 2. Political uncertainty and R&D intensity: Baseline results**"),
    subtitle = md("Dependent variable: R&D intensity")
  ) %>%
  cols_label(
    Variable = "Variable",
    ` (1) ` = "(1)",
    ` (2) ` = "(2)"
  ) %>%
  cols_align(align = "left", columns = "Variable") %>%
  cols_align(align = "center", columns = c(" (1) ", " (2) ")) %>%
  tab_style(
    style = cell_text(weight = "bold"),
    locations = cells_body(rows = row_type == "coef" & Variable != "")
  ) %>%
  tab_style(
    style = cell_text(style = "italic", color = "gray30"),
    locations = cells_body(rows = row_type == "se")
  ) %>%
  tab_style(
    style = cell_borders(sides = "top", color = "black", weight = px(1)),
    locations = cells_body(rows = row_type == "gof")
  ) %>%
  tab_options(
    table.font.size = px(18),
    heading.title.font.size = px(24),
    heading.subtitle.font.size = px(18),
    column_labels.font.weight = "bold",
    data_row.padding = px(8),
    table.border.top.width = px(2),
    table.border.bottom.width = px(2)
  ) %>%
  tab_source_note(md(
    "Standard errors in parentheses. Firm and year fixed effects included. State-clustered standard errors. Data period: 1976-2020."
  )) %>%
  cols_hide(columns = "row_type")

gt_tab2

# Optional exports (uncomment if desired)
# gtsave(gt_tab2, "table2_cols12_pretty.html")
# gtsave(gt_tab2, "table2_cols12_pretty.png")
# gtsave(gt_tab2, "table2_cols12_pretty.pdf")