Bank Selection into Emergency Lending Facilities: Extensive Margin Analysis

BTFP, Discount Window, and FHLB Borrowing During the March 2023 Banking Crisis

1 SECTION 0: SETUP AND DATA PREPARATION

1.1 Load Packages

library(tidyverse)
library(data.table)
library(lubridate)
library(fixest)
library(modelsummary)
library(kableExtra)
library(scales)
library(ggplot2)
library(patchwork)

cat("All packages loaded successfully.\n")

## All packages loaded successfully.

1.2 Helper Functions

# ==============================================================================
# HELPER FUNCTIONS
# ==============================================================================

# Winsorization function (2.5% / 97.5% by default)
winsorize <- function(x, probs = c(0.025, 0.975)) {
  if (all(is.na(x))) return(x)
  q <- quantile(x, probs = probs, na.rm = TRUE, names = FALSE)
  pmax(pmin(x, q[2]), q[1])
}

# Z-score standardization function
standardize_z <- function(x) {
  if (all(is.na(x))) return(x)
  (x - mean(x, na.rm = TRUE)) / sd(x, na.rm = TRUE)
}

# Safe division (avoid division by zero)
safe_div <- function(num, denom, default = NA_real_) {
  ifelse(is.na(denom) | denom == 0, default, num / denom)
}

# Size category function (3 categories)
create_size_category_3 <- function(assets_thousands) {
  assets_millions <- assets_thousands / 1000
  case_when(
    assets_millions >= 100000 ~ "Large (>$100B)",
    assets_millions >= 1000   ~ "Medium ($1B-$100B)",
    TRUE                      ~ "Small (<$1B)"
  )
}

size_levels_3 <- c("Small (<$1B)", "Medium ($1B-$100B)", "Large (>$100B)")

# Format P-Value with Stars
format_pval <- function(p) {
  case_when(
    is.na(p) ~ "",
    p < 0.01 ~ "***",
    p < 0.05 ~ "**",
    p < 0.10 ~ "*",
    TRUE ~ ""
  )
}

# ==============================================================================
# OUTPUT SAVING FUNCTIONS (Added for saving tables and figures)
# ==============================================================================

# Save kable/kableExtra table to HTML and LaTeX
save_kable_table <- function(tbl, filename, caption_text = "", notes_text = "") {
  # HTML version
  html_tbl <- tbl %>%
    kable(format = "html", caption = caption_text) %>%
    kable_styling(bootstrap_options = c("striped", "hover", "condensed"), 
                  full_width = FALSE, font_size = 10)
  save_kable(html_tbl, file = file.path(TABLE_PATH, paste0(filename, ".html")))
  
  # LaTeX version
  latex_tbl <- tbl %>%
    kable(format = "latex", caption = caption_text, booktabs = TRUE) %>%
    kable_styling(latex_options = c("hold_position", "scale_down"))
  if (notes_text != "") {
    latex_tbl <- latex_tbl %>% 
      footnote(general = notes_text, general_title = "Notes: ", threeparttable = TRUE)
  }
  writeLines(latex_tbl, file.path(TABLE_PATH, paste0(filename, ".tex")))
  cat("Saved:", filename, "(HTML + LaTeX)\n")
}

# Save modelsummary regression table to HTML and LaTeX
save_reg_table <- function(models, filename, title_text = "", notes_text = "",
                           coef_map_use = NULL, gof_map_use = NULL, add_rows_use = NULL, ...) {
  # HTML
  modelsummary(models, stars = c('*' = 0.1, '**' = 0.05, '***' = 0.01),
               coef_map = coef_map_use, gof_map = gof_map_use, add_rows = add_rows_use,
               title = title_text, notes = notes_text,
               output = file.path(TABLE_PATH, paste0(filename, ".html")), ...)
  # LaTeX
  modelsummary(models, stars = c('*' = 0.1, '**' = 0.05, '***' = 0.01),
               coef_map = coef_map_use, gof_map = gof_map_use, add_rows = add_rows_use,
               title = title_text, notes = notes_text,
               output = file.path(TABLE_PATH, paste0(filename, ".tex")), ...)
  cat("Saved:", filename, "(HTML + LaTeX)\n")
}

# Save ggplot figure to PDF and PNG
save_figure <- function(plot_obj, filename, width = 10, height = 6) {
  ggsave(file.path(FIG_PATH, paste0(filename, ".pdf")), plot = plot_obj, 
         width = width, height = height, device = "pdf")
  ggsave(file.path(FIG_PATH, paste0(filename, ".png")), plot = plot_obj, 
         width = width, height = height, dpi = 300, device = "png")
  cat("Saved:", filename, "(PDF + PNG)\n")
}

# Aliases for run_4spec_models and create_n_rows (used in later sections)
run_4spec_models <- function(data, dv, family_type = c("lpm","logit")) {
  family_type <- match.arg(family_type)
  forms <- list(
    "(1) Base" = build_formula(dv, "mtm_total + uninsured_lev + mtm_x_uninsured"),
    "(2) +Risk1" = build_formula(dv, "run_risk_1"),
    "(3) +Risk1,2" = build_formula(dv, "run_risk_1 + run_risk_2"),
    "(4) Risk2,3,4" = build_formula(dv, "run_risk_2 + run_risk_3 + run_risk_4")
  )
  if (family_type == "lpm") {
    lapply(forms, function(ff) feols(ff, data = data, vcov = "hetero"))
  } else {
    lapply(forms, function(ff) feglm(ff, data = data, family = binomial("logit"), vcov = "hetero"))
  }
}

create_n_rows <- function(data, dv, n_models = 4) {
  n_ones <- sum(data[[dv]] == 1, na.rm = TRUE)
  n_sample <- nrow(data)
  out <- data.frame(term = c(paste0("N (", dv, "=1)"), "N (Sample)"))
  for (i in 1:n_models) out[[paste0("(", i, ")")]] <- c(n_ones, n_sample)
  out
}

1.3 Paths and Key Dates

# ==============================================================================
# PATHS AND KEY DATES
# ==============================================================================

BASE_PATH <- "C:/Users/mferdo2/OneDrive - Louisiana State University/Finance_PhD/DW_Stigma_paper/Liquidity_project_2025"
DATA_PROC <- file.path(BASE_PATH, "01_data/processed")
OUTPUT_PATH <- file.path(BASE_PATH, "03_documentation/crisis_borrowing_result")
TABLE_PATH <- file.path(OUTPUT_PATH, "tables")
FIG_PATH <- file.path(OUTPUT_PATH, "figures")

for (path in c(TABLE_PATH, FIG_PATH)) {
  if (!dir.exists(path)) dir.create(path, recursive = TRUE)
}

# Baseline periods for different analyses
BASELINE_MAIN <- "2022Q4"    # For Acute, Post-Acute
BASELINE_ARB <- "2023Q3"     # For Arbitrage
BASELINE_WIND <- "2023Q4"    # For Wind-down

# Key dates
DATE_MAR01 <- as.Date("2023-03-01")
DATE_MAR09 <- as.Date("2023-03-09")
DATE_MAR10 <- as.Date("2023-03-10")
DATE_MAR12 <- as.Date("2023-03-12")
DATE_MAR13 <- as.Date("2023-03-13")
DATE_MAR14 <- as.Date("2023-03-14")

# Period definitions
ACUTE_START <- as.Date("2023-03-13")
ACUTE_END <- as.Date("2023-05-01")
POST_ACUTE_END <- as.Date("2023-10-31")
ARB_START <- as.Date("2023-11-01")
ARB_END <- as.Date("2024-01-24")
WIND_START <- as.Date("2024-01-25")
WIND_END <- as.Date("2024-03-11")
DW_DATA_END <- as.Date("2023-09-30")
OVERALL_START <- as.Date("2023-03-01")
OVERALL_END <- as.Date("2024-03-11")

1.4 Load Data

# ==============================================================================
# LOAD DATA
# ==============================================================================

call_q <- read_csv(file.path(DATA_PROC, "final_call_gsib.csv"), show_col_types = FALSE) %>%
  mutate(idrssd = as.character(idrssd))

btfp_loans_raw <- read_csv(file.path(DATA_PROC, "btfp_loan_bank_only.csv"), show_col_types = FALSE) %>%
  mutate(rssd_id = as.character(rssd_id), btfp_loan_date = mdy(btfp_loan_date))

dw_loans_raw <- read_csv(file.path(DATA_PROC, "dw_loan_bank_2023.csv"), show_col_types = FALSE) %>%
  mutate(rssd_id = as.character(rssd_id), dw_loan_date = ymd(dw_loan_date))

cat("=== DATA LOADED ===\n")

## === DATA LOADED ===

cat("Call Report:", nrow(call_q), "obs |", n_distinct(call_q$idrssd), "banks\n")

## Call Report: 75989 obs | 5074 banks

cat("BTFP Loans:", nrow(btfp_loans_raw), "| DW Loans:", nrow(dw_loans_raw), "\n")

## BTFP Loans: 6734 | DW Loans: 10008

1.5 Exclude Failed Banks and G-SIBs

# ==============================================================================
# EXCLUDE FAILED BANKS AND G-SIBs
# ==============================================================================

excluded_banks <- call_q %>%
  filter(period == BASELINE_MAIN, failed_bank == 1 | gsib == 1) %>%
  pull(idrssd)

cat("Excluded banks (failed + G-SIBs):", length(excluded_banks), "\n")

## Excluded banks (failed + G-SIBs): 41

btfp_loans <- btfp_loans_raw %>% filter(!rssd_id %in% excluded_banks)
dw_loans <- dw_loans_raw %>% filter(!rssd_id %in% excluded_banks)

1.6 Create Borrower Indicators by Period

# ==============================================================================
# HELPER FUNCTION: Create Borrower Indicator
# ==============================================================================

create_borrower_indicator <- function(loans_df, date_col, id_col, amount_col,
                                       start_date, end_date, prefix) {
  loans_df %>%
    filter(!!sym(date_col) >= start_date, !!sym(date_col) <= end_date) %>%
    group_by(!!sym(id_col)) %>%
    summarise(
      "{prefix}" := 1L,
      "{prefix}_amt" := sum(!!sym(amount_col), na.rm = TRUE),
      "{prefix}_first" := min(!!sym(date_col)),
      .groups = "drop"
    ) %>%
    rename(idrssd = !!sym(id_col))
}

# ==============================================================================
# BTFP BORROWERS BY PERIOD
# ==============================================================================

btfp_mar10 <- create_borrower_indicator(
  btfp_loans, "btfp_loan_date", "rssd_id", "btfp_loan_amount", 
  DATE_MAR10, DATE_MAR10, "btfp_mar10"
)

btfp_mar10_13 <- create_borrower_indicator(
  btfp_loans, "btfp_loan_date", "rssd_id", "btfp_loan_amount", 
  DATE_MAR10, DATE_MAR13, "btfp_mar10_13"
)

btfp_acute <- create_borrower_indicator(
  btfp_loans, "btfp_loan_date", "rssd_id", "btfp_loan_amount", 
  ACUTE_START, ACUTE_END, "btfp_acute"
)

btfp_post <- create_borrower_indicator(
  btfp_loans, "btfp_loan_date", "rssd_id", "btfp_loan_amount", 
  ACUTE_END + 1, POST_ACUTE_END, "btfp_post"
)

btfp_arb <- create_borrower_indicator(
  btfp_loans, "btfp_loan_date", "rssd_id", "btfp_loan_amount", 
  ARB_START, ARB_END, "btfp_arb"
)

btfp_wind <- create_borrower_indicator(
  btfp_loans, "btfp_loan_date", "rssd_id", "btfp_loan_amount", 
  WIND_START, WIND_END, "btfp_wind"
)

btfp_overall <- create_borrower_indicator(
  btfp_loans, "btfp_loan_date", "rssd_id", "btfp_loan_amount", 
  OVERALL_START, OVERALL_END, "btfp_overall"
)

# ==============================================================================
# DW BORROWERS BY PERIOD
# ==============================================================================

dw_prebtfp <- create_borrower_indicator(
  dw_loans, "dw_loan_date", "rssd_id", "dw_loan_amount", 
  DATE_MAR01, DATE_MAR12, "dw_prebtfp"
)

dw_mar10 <- create_borrower_indicator(
  dw_loans, "dw_loan_date", "rssd_id", "dw_loan_amount", 
  DATE_MAR10, DATE_MAR10, "dw_mar10"
)

dw_mar10_13 <- create_borrower_indicator(
  dw_loans, "dw_loan_date", "rssd_id", "dw_loan_amount", 
  DATE_MAR10, DATE_MAR13, "dw_mar10_13"
)

dw_acute <- create_borrower_indicator(
  dw_loans, "dw_loan_date", "rssd_id", "dw_loan_amount", 
  ACUTE_START, min(ACUTE_END, DW_DATA_END), "dw_acute"
)

dw_post <- create_borrower_indicator(
  dw_loans, "dw_loan_date", "rssd_id", "dw_loan_amount", 
  ACUTE_END + 1, min(POST_ACUTE_END, DW_DATA_END), "dw_post"
)

dw_overall <- create_borrower_indicator(
  dw_loans, "dw_loan_date", "rssd_id", "dw_loan_amount", 
  OVERALL_START, DW_DATA_END, "dw_overall"
)

cat("\n=== BORROWER COUNTS BY PERIOD ===\n")

## 
## === BORROWER COUNTS BY PERIOD ===

cat("BTFP: Mar10=", nrow(btfp_mar10), "| Mar10-13=", nrow(btfp_mar10_13), 
    "| Acute=", nrow(btfp_acute), "| Post=", nrow(btfp_post), 
    "| Arb=", nrow(btfp_arb), "| Wind=", nrow(btfp_wind), "\n")

## BTFP: Mar10= 0 | Mar10-13= 3 | Acute= 485 | Post= 811 | Arb= 797 | Wind= 237

cat("DW: Pre=", nrow(dw_prebtfp), "| Mar10=", nrow(dw_mar10), 
    "| Mar10-13=", nrow(dw_mar10_13), "| Acute=", nrow(dw_acute), "\n")

## DW: Pre= 106 | Mar10= 50 | Mar10-13= 94 | Acute= 417

1.7 Construct Analysis Variables

# ==============================================================================
# FUNCTION: Construct Analysis Variables
# ==============================================================================

construct_analysis_vars <- function(baseline_data) {
  
  baseline_data %>%
    mutate(
      # ========================================================================
      # RAW VARIABLES (original scale)
      # ========================================================================
      mtm_total_raw = mtm_loss_to_total_asset,
      mtm_btfp_raw = mtm_loss_omo_eligible_to_total_asset,
      mtm_other_raw = mtm_loss_non_omo_eligible_to_total_asset,
      uninsured_lev_raw = uninsured_deposit_to_total_asset,
      uninsured_share_raw = uninsured_to_deposit,
      mv_asset = mm_asset,
      
      ln_assets_raw = log(total_asset),
      cash_ratio_raw = cash_to_total_asset,
      securities_ratio_raw = security_to_total_asset,
      loan_ratio_raw = total_loan_to_total_asset,
      book_equity_ratio_raw = book_equity_to_total_asset,
      tier1_ratio_raw = tier1cap_to_total_asset,
      roa_raw = roa,
      fhlb_ratio_raw = fhlb_to_total_asset,
      loan_to_deposit_raw = loan_to_deposit,
      wholesale_raw = safe_div(
        fed_fund_purchase + repo + replace_na(other_borrowed_less_than_1yr, 0),
        total_liability, 0
      ) * 100,
      
      # Jiang et al. insolvency measures
      adjusted_equity_raw = book_equity_to_total_asset - mtm_loss_to_total_asset,
      mv_adjustment_raw = if_else(mv_asset == 0 | is.na(mv_asset), NA_real_, (total_asset / mv_asset) - 1),
      idcr_1_raw = safe_div(mv_asset - 0.5 * uninsured_deposit - insured_deposit, insured_deposit),
      idcr_2_raw = safe_div(mv_asset - 1.0 * uninsured_deposit - insured_deposit, insured_deposit),
      insolvency_1_raw = safe_div((total_asset - total_liability) - 0.5 * uninsured_deposit * mv_adjustment_raw, total_asset),
      insolvency_2_raw = safe_div((total_asset - total_liability) - 1.0 * uninsured_deposit * mv_adjustment_raw, total_asset),
      
      # ========================================================================
      # WINSORIZED VARIABLES (for outlier control)
      # ========================================================================
      mtm_total_w = winsorize(mtm_total_raw),
      mtm_btfp_w = winsorize(mtm_btfp_raw),
      mtm_other_w = winsorize(mtm_other_raw),
      uninsured_lev_w = winsorize(uninsured_lev_raw),
      uninsured_share_w = winsorize(uninsured_share_raw),
      
      ln_assets_w = winsorize(ln_assets_raw),
      cash_ratio_w = winsorize(cash_ratio_raw),
      securities_ratio_w = winsorize(securities_ratio_raw),
      loan_ratio_w = winsorize(loan_ratio_raw),
      book_equity_ratio_w = winsorize(book_equity_ratio_raw),
      tier1_ratio_w = winsorize(tier1_ratio_raw),
      roa_w = winsorize(roa_raw),
      fhlb_ratio_w = winsorize(fhlb_ratio_raw),
      loan_to_deposit_w = winsorize(loan_to_deposit_raw),
      wholesale_w = winsorize(wholesale_raw),
      
      # ========================================================================
      # Z-SCORE STANDARDIZED VARIABLES (for economic interpretation)
      # ========================================================================
      mtm_total = standardize_z(mtm_total_w),
      mtm_btfp = standardize_z(mtm_btfp_w),
      mtm_other = standardize_z(mtm_other_w),
      uninsured_lev = standardize_z(uninsured_lev_w),
      uninsured_share = standardize_z(uninsured_share_w),
      
      # Interaction term (standardized)
      mtm_x_uninsured = mtm_total * uninsured_lev,
      
      # Controls (standardized)
      ln_assets = standardize_z(ln_assets_w),
      cash_ratio = standardize_z(cash_ratio_w),
      securities_ratio = standardize_z(securities_ratio_w),
      loan_ratio = standardize_z(loan_ratio_w),
      book_equity_ratio = standardize_z(book_equity_ratio_w),
      tier1_ratio = standardize_z(tier1_ratio_w),
      roa = standardize_z(roa_w),
      fhlb_ratio = standardize_z(fhlb_ratio_w),
      loan_to_deposit = standardize_z(loan_to_deposit_w),
      wholesale = standardize_z(wholesale_w),
      
      # Jiang et al. measures (winsorized)
      adjusted_equity = winsorize(adjusted_equity_raw),
      mv_adjustment = winsorize(mv_adjustment_raw),
      idcr_1 = winsorize(idcr_1_raw),
      idcr_2 = winsorize(idcr_2_raw),
      insolvency_1 = winsorize(insolvency_1_raw),
      insolvency_2 = winsorize(insolvency_2_raw),
      
      # Insolvency indicators
      mtm_insolvent = as.integer(adjusted_equity < 0),
      insolvent_idcr_s50 = as.integer(idcr_1 < 0),
      insolvent_idcr_s100 = as.integer(idcr_2 < 0),
      insolvent_cap_s50 = as.integer(insolvency_1 < 0),
      insolvent_cap_s100 = as.integer(insolvency_2 < 0),
      
      # High/Low indicators
      high_uninsured = as.integer(uninsured_lev > median(uninsured_lev, na.rm = TRUE)),
      high_mtm_loss = as.integer(mtm_total > median(mtm_total, na.rm = TRUE)),
      
      # Size category
      size_cat = factor(create_size_category_3(total_asset), levels = size_levels_3),
      
      # State for clustering
      state = if("state" %in% names(.)) state else NA_character_,
      fed_district = if("fed_district" %in% names(.)) fed_district else NA_character_
    )
}

# ==============================================================================
# FUNCTION: Add Run Risk Dummies (based on medians)
# ==============================================================================

add_run_risk_dummies <- function(data) {
  
  medians <- data %>%
    summarise(
      median_mtm = median(mtm_total_w, na.rm = TRUE),
      median_uninsured = median(uninsured_lev_w, na.rm = TRUE)
    )
  
  data %>%
    mutate(
      # Run Risk Dummies (based on WINSORIZED values before standardization)
      # Risk 1: Low MTM & Low Uninsured (Reference Category)
      # Risk 2: Low MTM & High Uninsured
      # Risk 3: High MTM & Low Uninsured
      # Risk 4: High MTM & High Uninsured
      # Risk 4: High MTM & High Uninsured (Insolvent & Illiquid - DEATH ZONE)
      run_risk_1 = as.integer(mtm_total_w < medians$median_mtm & uninsured_lev_w < medians$median_uninsured),
      run_risk_2 = as.integer(mtm_total_w < medians$median_mtm & uninsured_lev_w >= medians$median_uninsured),
      run_risk_3 = as.integer(mtm_total_w >= medians$median_mtm & uninsured_lev_w < medians$median_uninsured),
      run_risk_4 = as.integer(mtm_total_w >= medians$median_mtm & uninsured_lev_w >= medians$median_uninsured),
      
      # Store medians for reference
      median_mtm_used = medians$median_mtm,
      median_uninsured_used = medians$median_uninsured
    )
}

# ==============================================================================
# CREATE BASELINE DATASETS
# ==============================================================================

df_2022q4 <- call_q %>%
  filter(period == BASELINE_MAIN, !idrssd %in% excluded_banks, 
         !is.na(omo_eligible) & omo_eligible > 0) %>%
  construct_analysis_vars() %>%
  add_run_risk_dummies()

df_2023q3 <- call_q %>%
  filter(period == BASELINE_ARB, !idrssd %in% excluded_banks,
         !is.na(omo_eligible) & omo_eligible > 0) %>%
  construct_analysis_vars() %>%
  add_run_risk_dummies()

df_2023q4 <- call_q %>%
  filter(period == BASELINE_WIND, !idrssd %in% excluded_banks,
         !is.na(omo_eligible) & omo_eligible > 0) %>%
  construct_analysis_vars() %>%
  add_run_risk_dummies()

cat("\n=== BASELINE DATASETS ===\n")

## 
## === BASELINE DATASETS ===

cat("2022Q4:", nrow(df_2022q4), "banks\n")

## 2022Q4: 4292 banks

cat("2023Q3:", nrow(df_2023q3), "banks\n")

## 2023Q3: 4214 banks

cat("2023Q4:", nrow(df_2023q4), "banks\n")

## 2023Q4: 4197 banks

cat("\n=== RUN RISK MEDIANS (2022Q4) ===\n")

## 
## === RUN RISK MEDIANS (2022Q4) ===

cat("Median MTM Loss (winsorized):", round(df_2022q4$median_mtm_used[1], 4), "%\n")

## Median MTM Loss (winsorized): 5.32 %

cat("Median Uninsured Leverage (winsorized):", round(df_2022q4$median_uninsured_used[1], 4), "%\n")

## Median Uninsured Leverage (winsorized): 22.25 %

1.8 Join Borrower Indicators and Create Clean Samples

# ==============================================================================
# CRITICAL: CLEAN SAMPLE CONSTRUCTION
# For binary outcomes: 0 = PURE NON-BORROWERS ONLY
# Exclude other facility users from the 0 category
# ==============================================================================

# Join all borrower indicators to baseline
join_all_borrowers <- function(base_df, btfp_df, dw_df, btfp_var, dw_var) {
  
  base_df %>%
    left_join(btfp_df %>% select(idrssd, starts_with(btfp_var)), by = "idrssd") %>%
    left_join(dw_df %>% select(idrssd, starts_with(dw_var)), by = "idrssd") %>%
    mutate(
      # Fill NAs with 0
      "{btfp_var}" := replace_na(!!sym(btfp_var), 0L),
      "{dw_var}" := replace_na(!!sym(dw_var), 0L),
      
      # FHLB user (from call report)
      fhlb_user = as.integer(abnormal_fhlb_borrowing_10pct == 1),
      
      # User group classification (for descriptive tables)
      user_group = case_when(
        !!sym(btfp_var) == 1 & !!sym(dw_var) == 1 ~ "Both",
        !!sym(btfp_var) == 1 & !!sym(dw_var) == 0 ~ "BTFP_Only",
        !!sym(btfp_var) == 0 & !!sym(dw_var) == 1 ~ "DW_Only",
        TRUE ~ "Neither"
      ),
      user_group = factor(user_group, levels = c("Neither", "BTFP_Only", "DW_Only", "Both")),
      
      # Combined indicators
      any_fed = as.integer(!!sym(btfp_var) == 1 | !!sym(dw_var) == 1),
      both_fed = as.integer(!!sym(btfp_var) == 1 & !!sym(dw_var) == 1),
      all_user = as.integer(!!sym(btfp_var) == 1 | !!sym(dw_var) == 1 | fhlb_user == 1),
      
      # PURE non-user (no BTFP, no DW, no FHLB) - for clean comparison
      non_user = as.integer(!!sym(btfp_var) == 0 & !!sym(dw_var) == 0 & fhlb_user == 0)
    )
}

# ==============================================================================
# CREATE PERIOD DATASETS
# ==============================================================================

# Acute Period (Mar 13 - May 1, 2023)
df_acute <- join_all_borrowers(df_2022q4, btfp_acute, dw_acute, "btfp_acute", "dw_acute") %>%
  mutate(
    # Intensive margin variables
    btfp_pct = ifelse(btfp_acute == 1 & btfp_acute_amt > 0, 
                      100 * btfp_acute_amt / (total_asset * 1000), NA_real_),
    dw_pct = ifelse(dw_acute == 1 & dw_acute_amt > 0, 
                    100 * dw_acute_amt / (total_asset * 1000), NA_real_),
    log_btfp_amt = ifelse(btfp_acute == 1 & btfp_acute_amt > 0, log(btfp_acute_amt), NA_real_),
    log_dw_amt = ifelse(dw_acute == 1 & dw_acute_amt > 0, log(dw_acute_amt), NA_real_)
  )

# March 10 (SVB Closure Day)
df_mar10 <- join_all_borrowers(df_2022q4, btfp_mar10, dw_mar10, "btfp_mar10", "dw_mar10")

# March 10-13 (BTFP Launch Window)
df_mar10_13 <- join_all_borrowers(df_2022q4, btfp_mar10_13, dw_mar10_13, "btfp_mar10_13", "dw_mar10_13")

# Pre-BTFP (Mar 1 - Mar 12)
btfp_prebtfp <- create_borrower_indicator(
  btfp_loans, "btfp_loan_date", "rssd_id", "btfp_loan_amount", 
  DATE_MAR01, DATE_MAR12, "btfp_prebtfp"
)
df_prebtfp <- join_all_borrowers(df_2022q4, btfp_prebtfp, dw_prebtfp, "btfp_prebtfp", "dw_prebtfp")

# Post-Acute Period (May 2 - Oct 31, 2023)
df_post <- df_2022q4 %>%
  left_join(btfp_post %>% select(idrssd, btfp_post, btfp_post_amt), by = "idrssd") %>%
  left_join(dw_post %>% select(idrssd, dw_post, dw_post_amt), by = "idrssd") %>%
  mutate(
    btfp_post = replace_na(btfp_post, 0L),
    dw_post = replace_na(dw_post, 0L),
    fhlb_user = as.integer(abnormal_fhlb_borrowing_10pct == 1),
    any_fed = as.integer(btfp_post == 1 | dw_post == 1),
    non_user = as.integer(btfp_post == 0 & dw_post == 0 & fhlb_user == 0),
    user_group = case_when(
      btfp_post == 1 & dw_post == 1 ~ "Both",
      btfp_post == 1 ~ "BTFP_Only",
      dw_post == 1 ~ "DW_Only",
      TRUE ~ "Neither"
    ),
    user_group = factor(user_group, levels = c("Neither", "BTFP_Only", "DW_Only", "Both"))
  )

# Arbitrage Period (Nov 1, 2023 - Jan 24, 2024) - 2023Q3 Baseline
df_arb <- df_2023q3 %>%
  left_join(btfp_arb %>% select(idrssd, btfp_arb, btfp_arb_amt), by = "idrssd") %>%
  mutate(
    btfp_arb = replace_na(btfp_arb, 0L),
    fhlb_user = as.integer(abnormal_fhlb_borrowing_10pct == 1),
    any_fed = btfp_arb,
    non_user = as.integer(btfp_arb == 0 & fhlb_user == 0),
    user_group = factor(ifelse(btfp_arb == 1, "BTFP_Only", "Neither"),
                        levels = c("Neither", "BTFP_Only", "DW_Only", "Both"))
  )

# Wind-down Period (Jan 25 - Mar 11, 2024) - 2023Q4 Baseline
df_wind <- df_2023q4 %>%
  left_join(btfp_wind %>% select(idrssd, btfp_wind, btfp_wind_amt), by = "idrssd") %>%
  mutate(
    btfp_wind = replace_na(btfp_wind, 0L),
    fhlb_user = as.integer(abnormal_fhlb_borrowing_10pct == 1),
    any_fed = btfp_wind,
    non_user = as.integer(btfp_wind == 0 & fhlb_user == 0),
    user_group = factor(ifelse(btfp_wind == 1, "BTFP_Only", "Neither"),
                        levels = c("Neither", "BTFP_Only", "DW_Only", "Both"))
  )

# Overall Period
df_overall <- join_all_borrowers(df_2022q4, btfp_overall, dw_overall, "btfp_overall", "dw_overall")

cat("\n=== USER GROUP COUNTS (ACUTE PERIOD) ===\n")

## 
## === USER GROUP COUNTS (ACUTE PERIOD) ===

print(table(df_acute$user_group))

## 
##   Neither BTFP_Only   DW_Only      Both 
##      3531       368       299        94

cat("\nPure Non-Users:", sum(df_acute$non_user), "\n")

## 
## Pure Non-Users: 3285

1.9 Model Setup

# ==============================================================================
# MODEL SETTINGS (DEFINE ONCE, USE EVERYWHERE)
# ==============================================================================

# Standard Controls (Z-score standardized)
CONTROLS <- "ln_assets + cash_ratio + loan_to_deposit + book_equity_ratio + wholesale + roa"

# Coefficient Labels for Output
COEF_MAP <- c(
  # Key explanatory (standardized - interpret as 1 SD change)
  "mtm_total" = "MTM Loss (z)",
  "mtm_btfp" = "MTM Loss OMO (z)",
  "uninsured_lev" = "Uninsured Lev (z)",
  "mtm_x_uninsured" = "MTM × Uninsured",
  "mtm_total:uninsured_lev" = "MTM × Uninsured",
  "mtm_btfp:uninsured_lev" = "MTM_btfp × Uninsured Lev",
  
  # Run Risk dummies (Reference = Risk 1: Low MTM, Low Uninsured)
  "run_risk_1" = "Risk 1: $<$ Med MTM \\& $<$ Med Uninsured",
  "run_risk_2" = "Risk 2: $<$ Med MTM \\& $>$ Med Uninsured",
  "run_risk_3" = "Risk 3: $>$ Med MTM \\& $<$ Med Uninsured",
  "run_risk_4" = "Risk 4: $>$ Med MTM \\& $>$ Med Uninsured",
  
  # Controls (standardized)
  "ln_assets" = "Log(Assets) (z)",
  "cash_ratio" = "Cash Ratio (z)",
  "loan_to_deposit" = "Loan/Deposit (z)",
  "book_equity_ratio" = "Book Equity (z)",
  "wholesale" = "Wholesale (z)",
  "roa" = "ROA (z)"
)

# Helper Function: Build formula dynamically
build_formula <- function(dv, explanatory, controls = CONTROLS) {
  as.formula(paste(dv, "~", explanatory, "+", controls))
}

# GOF statistics
gof_lpm <- c("nobs", "r.squared", "adj.r.squared")
gof_logit <- c("nobs", "logLik", "AIC")

2 SECTION 1: DESCRIPTIVE STATISTICS

# ==============================================================================
# FUNCTION: Create Descriptive Statistics Table
# ==============================================================================

create_desc_stats <- function(data, vars, var_labels = NULL) {
  
  stats_list <- map_dfr(vars, function(v) {
    x <- data[[v]]
    x <- x[!is.na(x)]
    
    tibble(
      Variable = ifelse(!is.null(var_labels) && v %in% names(var_labels), 
                        var_labels[[v]], v),
      N = length(x),
      Mean = mean(x),
      SD = sd(x),
      Min = min(x),
      P25 = quantile(x, 0.25),
      Median = quantile(x, 0.50),
      P75 = quantile(x, 0.75),
      Max = max(x)
    )
  })
  
  return(stats_list)
}

# Variable labels for raw (non-standardized) variables
var_labels_raw <- c(
  "mtm_total_w" = "MTM Loss / Total Assets (%)",
  "mtm_btfp_w" = "MTM Loss (OMO-Eligible) / Total Assets (%)",
  "uninsured_lev_w" = "Uninsured Deposits / Total Assets (%)",
  "uninsured_share_w" = "Uninsured Deposits / Total Deposits (%)",
  "ln_assets_w" = "Log(Total Assets)",
  "cash_ratio_w" = "Cash / Total Assets (%)",
  "securities_ratio_w" = "Securities / Total Assets (%)",
  "loan_ratio_w" = "Loans / Total Assets (%)",
  "book_equity_ratio_w" = "Book Equity / Total Assets (%)",
  "roa_w" = "Return on Assets (%)",
  "fhlb_ratio_w" = "FHLB Advances / Total Assets (%)",
  "loan_to_deposit_w" = "Loans / Deposits (%)",
  "wholesale_w" = "Wholesale Funding / Total Liabilities (%)"
)

desc_vars <- c(
  "mtm_total_w", "mtm_btfp_w", "uninsured_lev_w", "uninsured_share_w",
  "ln_assets_w", "cash_ratio_w", "securities_ratio_w", "loan_ratio_w",
  "book_equity_ratio_w", "roa_w", "fhlb_ratio_w", "loan_to_deposit_w", "wholesale_w"
)

2.1 Table 1A: Summary Statistics (2022Q4 Baseline)

desc_2022q4 <- create_desc_stats(df_2022q4, desc_vars, var_labels_raw)

desc_2022q4 %>%
  mutate(across(where(is.numeric), ~round(., 3))) %>%
  kable(format = "html", caption = "Table 1A: Summary Statistics (2022Q4 Baseline)",
        col.names = c("Variable", "N", "Mean", "SD", "Min", "P25", "Median", "P75", "Max")) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE, font_size = 10) %>%
  pack_rows("Key Explanatory Variables", 1, 4) %>%
  pack_rows("Control Variables", 5, 13)

Table 1A: Summary Statistics (2022Q4 Baseline)

Variable	N	Mean	SD	Min	P25	Median	P75	Max
Key Explanatory Variables
MTM Loss / Total Assets (%)	4282	5.462	2.115	1.371	3.863	5.320	6.975	9.881
MTM Loss (OMO-Eligible) / Total Assets (%)	4282	0.655	0.712	0.000	0.136	0.405	0.932	2.941
Uninsured Deposits / Total Assets (%)	4292	23.463	11.328	4.212	15.285	22.247	30.353	51.992
Uninsured Deposits / Total Deposits (%)	4292	27.423	13.283	5.518	17.867	25.801	35.108	62.126
Control Variables
Log(Total Assets)	4292	12.869	1.370	10.445	11.903	12.721	13.650	16.433
Cash / Total Assets (%)	4292	7.789	7.503	1.098	2.667	5.030	10.051	33.276
Securities / Total Assets (%)	4292	25.436	14.934	2.551	13.704	23.323	35.090	61.559
Loans / Total Assets (%)	4292	60.259	16.657	18.300	49.604	62.303	73.421	85.811
Book Equity / Total Assets (%)	4292	9.498	3.886	3.318	7.143	8.842	10.888	22.853
Return on Assets (%)	4292	1.032	0.519	-0.112	0.704	1.014	1.319	2.381
FHLB Advances / Total Assets (%)	4292	2.515	3.644	0.000	0.000	0.346	4.073	13.950
Loans / Deposits (%)	4292	70.645	21.258	21.492	56.561	71.818	86.500	109.103
Wholesale Funding / Total Liabilities (%)	4292	0.828	1.842	0.000	0.000	0.000	0.533	8.087

# SAVE TABLE 1A
save_kable_table(desc_2022q4 %>% mutate(across(where(is.numeric), ~round(., 3))),
                 "Table_1A_Summary_Stats_2022Q4", 
                 "Table 1A: Summary Statistics (2022Q4 Baseline)",
                 "Variables winsorized at 2.5/97.5 percentiles.")

Saved: Table_1A_Summary_Stats_2022Q4 (HTML + LaTeX)

2.2 Table 1B: Summary Statistics (2023Q3 Baseline - Arbitrage)

desc_2023q3 <- create_desc_stats(df_2023q3, desc_vars, var_labels_raw)

desc_2023q3 %>%
  mutate(across(where(is.numeric), ~round(., 3))) %>%
  kable(format = "html", caption = "Table 1B: Summary Statistics (2023Q3 Baseline - Arbitrage)",
        col.names = c("Variable", "N", "Mean", "SD", "Min", "P25", "Median", "P75", "Max")) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE, font_size = 10) %>%
  pack_rows("Key Explanatory Variables", 1, 4) %>%
  pack_rows("Control Variables", 5, 13)

Table 1B: Summary Statistics (2023Q3 Baseline - Arbitrage)

Variable	N	Mean	SD	Min	P25	Median	P75	Max
Key Explanatory Variables
MTM Loss / Total Assets (%)	4197	5.381	2.113	1.375	3.759	5.224	6.905	9.760
MTM Loss (OMO-Eligible) / Total Assets (%)	4197	0.647	0.710	0.000	0.128	0.395	0.920	2.917
Uninsured Deposits / Total Assets (%)	4214	21.020	10.253	3.279	13.550	19.877	27.236	46.585
Uninsured Deposits / Total Deposits (%)	4214	24.967	12.075	4.477	16.159	23.638	32.087	55.269
Control Variables
Log(Total Assets)	4214	12.898	1.383	10.424	11.924	12.749	13.674	16.502
Cash / Total Assets (%)	4214	7.455	6.701	1.084	2.742	5.148	9.633	30.193
Securities / Total Assets (%)	4214	23.488	14.226	2.357	12.542	21.276	32.294	58.940
Loans / Total Assets (%)	4214	62.426	16.096	19.210	52.990	65.032	74.991	85.879
Book Equity / Total Assets (%)	4214	9.648	4.070	3.101	7.170	8.995	11.131	23.522
Return on Assets (%)	4214	1.090	0.719	-0.274	0.630	0.996	1.443	3.304
FHLB Advances / Total Assets (%)	4214	3.005	4.026	0.000	0.000	1.105	4.927	15.327
Loans / Deposits (%)	4214	74.319	20.864	23.348	61.302	76.806	89.898	112.184
Wholesale Funding / Total Liabilities (%)	4214	1.474	2.790	0.000	0.000	0.000	1.655	11.260

# SAVE TABLE 1B
save_kable_table(desc_2023q3 %>% mutate(across(where(is.numeric), ~round(., 3))),
                 "Table_1B_Summary_Stats_2023Q3", 
                 "Table 1B: Summary Statistics (2023Q3 Baseline - Arbitrage)",
                 "Variables winsorized at 2.5/97.5 percentiles.")

Saved: Table_1B_Summary_Stats_2023Q3 (HTML + LaTeX)

2.3 Table 1C: Summary Statistics (2023Q4 Baseline - Wind-down)

desc_2023q4 <- create_desc_stats(df_2023q4, desc_vars, var_labels_raw)

desc_2023q4 %>%
  mutate(across(where(is.numeric), ~round(., 3))) %>%
  kable(format = "html", caption = "Table 1C: Summary Statistics (2023Q4 Baseline - Wind-down)",
        col.names = c("Variable", "N", "Mean", "SD", "Min", "P25", "Median", "P75", "Max")) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE, font_size = 10) %>%
  pack_rows("Key Explanatory Variables", 1, 4) %>%
  pack_rows("Control Variables", 5, 13)

Table 1C: Summary Statistics (2023Q4 Baseline - Wind-down)

Variable	N	Mean	SD	Min	P25	Median	P75	Max
Key Explanatory Variables
MTM Loss / Total Assets (%)	4179	5.290	2.087	1.326	3.704	5.133	6.809	9.596
MTM Loss (OMO-Eligible) / Total Assets (%)	4179	0.639	0.700	0.000	0.127	0.387	0.910	2.858
Uninsured Deposits / Total Assets (%)	4197	21.108	10.417	3.518	13.435	19.943	27.172	47.013
Uninsured Deposits / Total Deposits (%)	4197	25.211	12.294	4.746	16.149	23.828	32.375	56.255
Control Variables
Log(Total Assets)	4197	12.920	1.379	10.457	11.942	12.767	13.704	16.500
Cash / Total Assets (%)	4197	7.665	6.821	1.137	2.799	5.310	10.028	30.833
Securities / Total Assets (%)	4197	23.487	14.247	2.442	12.339	21.263	32.191	59.061
Loans / Total Assets (%)	4197	62.467	15.988	19.973	52.946	65.306	75.052	85.791
Book Equity / Total Assets (%)	4197	10.169	3.840	4.416	7.851	9.424	11.493	23.761
Return on Assets (%)	4197	1.024	0.686	-0.333	0.587	0.950	1.373	3.038
FHLB Advances / Total Assets (%)	4197	2.987	4.021	0.000	0.000	1.056	4.920	15.469
Loans / Deposits (%)	4197	74.876	20.846	23.590	61.672	77.204	90.118	112.320
Wholesale Funding / Total Liabilities (%)	4197	1.604	2.932	0.000	0.000	0.000	1.989	11.749

# SAVE TABLE 1C
save_kable_table(desc_2023q4 %>% mutate(across(where(is.numeric), ~round(., 3))),
                 "Table_1C_Summary_Stats_2023Q4", 
                 "Table 1C: Summary Statistics (2023Q4 Baseline - Wind-down)",
                 "Variables winsorized at 2.5/97.5 percentiles.")

Saved: Table_1C_Summary_Stats_2023Q4 (HTML + LaTeX)

3 SECTION 2: SUMMARY TABLES BY USER GROUP (WITH P-VALUES)

# ==============================================================================
# FUNCTION: Create User Group Comparison Table with P-Values and Significance Stars
# ==============================================================================

create_user_group_table <- function(data, period_name, comparison_vars, var_labels, 
                                    return_df = FALSE) {
  
  # Get group counts
  group_counts <- data %>%
    group_by(user_group) %>%
    summarise(N = n(), .groups = "drop") %>%
    arrange(user_group)
  
  groups <- as.character(group_counts$user_group)
  
  # Calculate stats for each variable and group
  results <- map_dfr(comparison_vars, function(v) {
    
    var_label <- ifelse(v %in% names(var_labels), var_labels[[v]], v)
    row_data <- tibble(Variable = var_label)
    
    for (g in groups) {
      x <- data %>% filter(user_group == g) %>% pull(!!sym(v))
      x <- x[!is.na(x)]
      row_data[[paste0(g, "_Mean")]] <- mean(x)
      row_data[[paste0(g, "_SD")]] <- sd(x)
    }
    
    # T-test: Users vs Neither (with p-value and significance stars)
    if ("Neither" %in% groups) {
      x_neither <- data %>% filter(user_group == "Neither") %>% pull(!!sym(v))
      x_users <- data %>% filter(user_group != "Neither") %>% pull(!!sym(v))
      x_neither <- x_neither[!is.na(x_neither)]
      x_users <- x_users[!is.na(x_users)]
      
      if (length(x_neither) > 1 && length(x_users) > 1) {
        ttest <- t.test(x_users, x_neither)
        row_data$Diff <- mean(x_users) - mean(x_neither)
        row_data$T_Stat <- as.numeric(ttest$statistic)
        row_data$P_Value <- as.numeric(ttest$p.value)
        row_data$Sig <- format_pval(row_data$P_Value)
      } else {
        row_data$Diff <- NA_real_
        row_data$T_Stat <- NA_real_
        row_data$P_Value <- NA_real_
        row_data$Sig <- ""
      }
    }
    
    return(row_data)
  })
  
  # Return dataframe if requested (for saving)
  if (return_df) {
    return(list(results = results, group_counts = group_counts))
  }
  
  # Create header with N below group names
  header_spec <- c(" " = 1)
  for (i in 1:nrow(group_counts)) {
    g <- as.character(group_counts$user_group[i])
    n <- group_counts$N[i]
    header_spec[paste0(g, "\n(N=", n, ")")] <- 2
  }
  if ("Neither" %in% groups) {
    header_spec["Diff (Users-Neither)"] <- 4
  }
  
  # Format and display
  results %>%
    mutate(across(where(is.numeric), ~round(., 3))) %>%
    kable(format = "html", caption = paste0("User Group Comparison: ", period_name),
          col.names = c("Variable", 
                        rep(c("Mean", "SD"), length(groups)),
                        "Diff", "t-stat", "p-value", "Sig")) %>%
    kable_styling(bootstrap_options = c("striped", "hover", "condensed"), 
                  full_width = FALSE, font_size = 9) %>%
    add_header_above(header_spec)
}

# ==============================================================================
# FUNCTION: Save User Group Comparison Table
# ==============================================================================

save_user_group_table <- function(data, period_name, comparison_vars, var_labels, 
                                  filename, notes_text) {
  
  # Get the data
  tbl_data <- create_user_group_table(data, period_name, comparison_vars, var_labels, 
                                      return_df = TRUE)
  results <- tbl_data$results
  
  # Format for saving
  results_formatted <- results %>%
    mutate(across(where(is.numeric), ~round(., 3)))
  
  # Save HTML
  html_tbl <- results_formatted %>%
    kable(format = "html", caption = paste0("User Group Comparison: ", period_name)) %>%
    kable_styling(bootstrap_options = c("striped", "hover", "condensed"), 
                  full_width = FALSE, font_size = 9)
  save_kable(html_tbl, file = file.path(TABLE_PATH, paste0(filename, ".html")))
  
  # Save LaTeX
  latex_tbl <- results_formatted %>%
    kable(format = "latex", caption = paste0("User Group Comparison: ", period_name), 
          booktabs = TRUE) %>%
    kable_styling(latex_options = c("hold_position", "scale_down")) %>%
    footnote(general = notes_text, general_title = "Notes: ", threeparttable = TRUE)
  writeLines(latex_tbl, file.path(TABLE_PATH, paste0(filename, ".tex")))
  
  cat("Saved:", filename, "(HTML + LaTeX)\n")
}

3.1 Table 2A: Acute Period (Mar 13 - May 1, 2023)

create_user_group_table(df_acute, "Acute Period (Mar 13 - May 1, 2023)", 
                        desc_vars, var_labels_raw)

User Group Comparison: Acute Period (Mar 13 - May 1, 2023)

	Neither (N=3531)		BTFP_Only (N=368)		DW_Only (N=299)		Both (N=94)		Diff (Users-Neither)
Variable	Mean	SD	Mean	SD	Mean	SD	Mean	SD	Diff	t-stat	p-value	Sig
MTM Loss / Total Assets (%)	5.354	2.136	6.162	1.938	5.747	1.964	5.865	1.781	0.608	7.703	0.000	***
MTM Loss (OMO-Eligible) / Total Assets (%)	0.625	0.697	0.837	0.772	0.718	0.733	0.876	0.819	0.170	5.646	0.000	***
Uninsured Deposits / Total Assets (%)	22.709	11.217	25.938	10.702	26.705	11.041	31.791	12.361	4.253	9.506	0.000	***
Uninsured Deposits / Total Deposits (%)	26.506	13.109	30.401	12.688	31.438	13.129	37.436	14.488	5.172	9.776	0.000	***
Log(Total Assets)	12.686	1.305	13.477	1.261	13.801	1.338	14.396	1.451	1.032	19.272	0.000	***
Cash / Total Assets (%)	8.333	7.764	4.694	4.530	6.136	6.440	4.738	5.314	-3.067	-12.888	0.000	***
Securities / Total Assets (%)	25.241	15.187	28.037	13.085	23.897	14.001	27.478	13.851	1.101	1.974	0.049	**
Loans / Total Assets (%)	59.794	17.066	61.429	13.766	63.801	15.021	61.894	14.822	2.625	4.399	0.000	***
Book Equity / Total Assets (%)	9.699	4.020	8.232	2.988	9.060	3.150	8.289	2.457	-1.135	-8.826	0.000	***
Return on Assets (%)	1.022	0.531	1.051	0.431	1.119	0.473	1.076	0.470	0.059	3.158	0.002	***
FHLB Advances / Total Assets (%)	2.299	3.527	3.800	4.151	3.091	3.782	3.747	3.956	1.215	7.765	0.000	***
Loans / Deposits (%)	69.967	21.631	72.409	18.447	75.461	19.852	73.860	19.260	3.820	4.876	0.000	***
Wholesale Funding / Total Liabilities (%)	0.733	1.732	1.254	2.259	1.156	2.160	1.673	2.350	0.534	6.197	0.000	***

# SAVE TABLE 2A
save_user_group_table(df_acute, "Acute Period (Mar 13 - May 1, 2023)",
                      desc_vars, var_labels_raw,
                      "Table_2A_UserGroup_Acute",
                      "This table compares bank characteristics across facility user groups during the acute crisis period. User groups: Neither = banks that did not access BTFP, DW, or FHLB; BTFP_Only = BTFP users only; DW_Only = Discount Window users only; Both = used both facilities. Diff = difference between all users and non-users. t-statistics from two-sample t-test. Significance levels: *** p<0.01, ** p<0.05, * p<0.10. Variables from 2022Q4 Call Reports, winsorized at 2.5/97.5 percentiles.")

Saved: Table_2A_UserGroup_Acute (HTML + LaTeX)

3.2 Table 2B: March 10 (SVB Closure Day)

create_user_group_table(df_mar10, "March 10, 2023 (SVB Closure Day)", 
                        desc_vars, var_labels_raw)

User Group Comparison: March 10, 2023 (SVB Closure Day)

	Neither (N=4245)		DW_Only (N=47)		Diff (Users-Neither)
Variable	Mean	SD	Mean	SD	Diff	t-stat	p-value	Sig
MTM Loss / Total Assets (%)	5.460	2.116	5.686	1.983	0.226	0.776	0.442
MTM Loss (OMO-Eligible) / Total Assets (%)	0.653	0.710	0.889	0.840	0.236	1.920	0.061
Uninsured Deposits / Total Assets (%)	23.397	11.308	29.484	11.644	6.087	3.565	0.001	***
Uninsured Deposits / Total Deposits (%)	27.336	13.246	35.289	14.393	7.954	3.771	0.000	***
Log(Total Assets)	12.851	1.357	14.477	1.613	1.627	6.888	0.000	***
Cash / Total Assets (%)	7.840	7.523	3.177	2.776	-4.664	-11.074	0.000	***
Securities / Total Assets (%)	25.430	14.940	25.973	14.553	0.543	0.254	0.800
Loans / Total Assets (%)	60.209	16.671	64.755	14.862	4.545	2.082	0.043	**
Book Equity / Total Assets (%)	9.512	3.896	8.259	2.553	-1.253	-3.322	0.002	***
Return on Assets (%)	1.032	0.520	1.013	0.365	-0.019	-0.356	0.723
FHLB Advances / Total Assets (%)	2.494	3.634	4.346	4.062	1.852	3.112	0.003	***
Loans / Deposits (%)	70.568	21.272	77.529	18.934	6.960	2.503	0.016	**
Wholesale Funding / Total Liabilities (%)	0.813	1.829	2.198	2.447	1.385	3.867	0.000	***

# SAVE TABLE 2B
save_user_group_table(df_mar10, "March 10, 2023 (SVB Closure Day)",
                      desc_vars, var_labels_raw,
                      "Table_2B_UserGroup_Mar10",
                      "This table compares bank characteristics for banks that accessed emergency facilities on March 10, 2023 (SVB closure day). The BTFP was not yet operational; DW was the only Fed facility available. Diff = difference between all users and non-users. t-statistics from two-sample t-test. Significance levels: *** p<0.01, ** p<0.05, * p<0.10. Variables from 2022Q4 Call Reports.")

Saved: Table_2B_UserGroup_Mar10 (HTML + LaTeX)

3.3 Table 2C: Pre-BTFP (Mar 1 - Mar 12)

create_user_group_table(df_prebtfp, "Pre-BTFP Period (Mar 1-12, 2023)", 
                        desc_vars, var_labels_raw)

User Group Comparison: Pre-BTFP Period (Mar 1-12, 2023)

	Neither (N=4192)		DW_Only (N=100)		Diff (Users-Neither)
Variable	Mean	SD	Mean	SD	Diff	t-stat	p-value	Sig
MTM Loss / Total Assets (%)	5.454	2.113	5.829	2.169	0.375	1.710	0.090
MTM Loss (OMO-Eligible) / Total Assets (%)	0.652	0.712	0.778	0.730	0.126	1.708	0.091
Uninsured Deposits / Total Assets (%)	23.360	11.310	27.799	11.291	4.439	3.885	0.000	***
Uninsured Deposits / Total Deposits (%)	27.288	13.242	33.068	13.841	5.780	4.131	0.000	***
Log(Total Assets)	12.836	1.350	14.234	1.493	1.398	9.270	0.000	***
Cash / Total Assets (%)	7.879	7.532	4.042	4.890	-3.836	-7.632	0.000	***
Securities / Total Assets (%)	25.446	14.945	25.006	14.531	-0.440	-0.299	0.765
Loans / Total Assets (%)	60.145	16.695	65.058	14.263	4.913	3.390	0.001	***
Book Equity / Total Assets (%)	9.515	3.901	8.790	3.112	-0.725	-2.287	0.024	**
Return on Assets (%)	1.030	0.521	1.101	0.418	0.070	1.652	0.102
FHLB Advances / Total Assets (%)	2.485	3.630	3.734	4.021	1.249	3.076	0.003	***
Loans / Deposits (%)	70.479	21.290	77.591	18.659	7.112	3.754	0.000	***
Wholesale Funding / Total Liabilities (%)	0.799	1.817	2.029	2.414	1.230	5.062	0.000	***

# SAVE TABLE 2C
save_user_group_table(df_prebtfp, "Pre-BTFP Period (Mar 1-12, 2023)",
                      desc_vars, var_labels_raw,
                      "Table_2C_UserGroup_PreBTFP",
                      "This table compares bank characteristics for banks that accessed emergency facilities before BTFP announcement (March 12, 2023). Only the Discount Window was available as a Fed facility during this period. Diff = difference between all users and non-users. t-statistics from two-sample t-test. Significance levels: *** p<0.01, ** p<0.05, * p<0.10. Variables from 2022Q4 Call Reports.")

Saved: Table_2C_UserGroup_PreBTFP (HTML + LaTeX)

3.4 Table 2D: Arbitrage Period (2023Q3 Baseline)

create_user_group_table(df_arb, "Arbitrage Period (Nov 2023 - Jan 2024) [2023Q3 Baseline]", 
                        desc_vars, var_labels_raw)

User Group Comparison: Arbitrage Period (Nov 2023 - Jan 2024) [2023Q3 Baseline]

	Neither (N=3448)		BTFP_Only (N=766)		Diff (Users-Neither)
Variable	Mean	SD	Mean	SD	Diff	t-stat	p-value	Sig
MTM Loss / Total Assets (%)	5.249	2.131	5.969	1.924	0.720	9.174	0	***
MTM Loss (OMO-Eligible) / Total Assets (%)	0.612	0.690	0.802	0.773	0.190	6.256	0	***
Uninsured Deposits / Total Assets (%)	20.671	10.362	22.594	9.597	1.923	4.943	0	***
Uninsured Deposits / Total Deposits (%)	24.513	12.185	27.010	11.353	2.497	5.432	0	***
Log(Total Assets)	12.766	1.369	13.488	1.291	0.722	13.846	0	***
Cash / Total Assets (%)	8.031	7.033	4.859	4.020	-3.172	-16.850	0	***
Securities / Total Assets (%)	23.140	14.604	25.050	12.272	1.909	3.756	0	***
Loans / Total Assets (%)	62.048	16.715	64.127	12.818	2.079	3.824	0	***
Book Equity / Total Assets (%)	9.912	4.209	8.460	3.111	-1.451	-10.885	0	***
Return on Assets (%)	1.131	0.742	0.906	0.569	-0.225	-9.302	0	***
FHLB Advances / Total Assets (%)	2.876	3.976	3.584	4.196	0.708	4.265	0	***
Loans / Deposits (%)	73.680	21.529	77.198	17.279	3.518	4.859	0	***
Wholesale Funding / Total Liabilities (%)	1.087	2.472	3.213	3.408	2.126	16.334	0	***

# SAVE TABLE 2D
save_user_group_table(df_arb, "Arbitrage Period (Nov 2023 - Jan 2024)",
                      desc_vars, var_labels_raw,
                      "Table_2D_UserGroup_Arbitrage",
                      "This table compares bank characteristics during the BTFP arbitrage period (Nov 1, 2023 - Jan 24, 2024). During this period, BTFP borrowing rate fell below the Fed Funds rate, creating arbitrage opportunities. DW data not available for this period. Diff = difference between BTFP users and non-users. t-statistics from two-sample t-test. Significance levels: *** p<0.01, ** p<0.05, * p<0.10. Variables from 2023Q3 Call Reports.")

Saved: Table_2D_UserGroup_Arbitrage (HTML + LaTeX)

4 SECTION 3: EXTENSIVE MARGIN ANALYSIS - ACUTE PERIOD

4.1 3.1 Acute Period: All Specifications BTFP

# ==============================================================================
# ACUTE PERIOD ANALYSIS
# Sample Construction: DV=1 users vs. DV=0 PURE non-users only
# ==============================================================================
# -----------------------------------------------------------------------------
# 0) Helper functions: build the 4 specs + optional single spec (Risk2,3,4 only)
# -----------------------------------------------------------------------------

# Formula helpers (keeps all model definitions consistent)
f_base <- function(dv) build_formula(dv, "mtm_total + uninsured_lev + mtm_x_uninsured")
f_r1   <- function(dv) build_formula(dv, "run_risk_1")
f_r12  <- function(dv) build_formula(dv, "run_risk_1 + run_risk_2")
f_r234 <- function(dv) build_formula(dv, "run_risk_2 + run_risk_3 + run_risk_4")

# Run the 4-spec set (Base, Risk1, Risk1+2, Risk2+3+4)
run_4spec <- function(data, dv, family_type = c("lpm","logit")) {
  family_type <- match.arg(family_type)

  forms <- list(
    "(1) Base"      = f_base(dv),
    "(2) +Risk1"    = f_r1(dv),
    "(3) +Risk1,2"  = f_r12(dv),
    "(4) Risk2,3,4" = f_r234(dv)
  )

  if (family_type == "lpm") {
    lapply(forms, function(ff) feols(ff, data = data, vcov = "hetero"))
  } else {
    lapply(forms, function(ff) feglm(ff, data = data, family = binomial("logit"), vcov = "hetero"))
  }
}

# Run only Risk2+3+4 (used for AnyFed column in Table 1/2)
run_r234_only <- function(data, dv, family_type = c("lpm","logit")) {
  family_type <- match.arg(family_type)
  ff <- f_r234(dv)

  if (family_type == "lpm") {
    feols(ff, data = data, vcov = "hetero")
  } else {
    feglm(ff, data = data, family = binomial("logit"), vcov = "hetero")
  }
}

# Add N rows (DV=1 count + sample size) for a given dataset and DV
add_n_rows <- function(data, dv, n_models) {
  n_ones   <- sum(data[[dv]] == 1, na.rm = TRUE)
  n_sample <- nrow(data)

  out <- data.frame(term = c(paste0("N (", dv, "=1)"), "N (Sample)"))
  for (i in 1:n_models) out[[paste0("(", i, ")")]] <- c(n_ones, n_sample)
  out
}

4.1.1 TABLE 1 (LPM) - Acute Period

# -----------------------------------------------------------------------------
# 1) TABLE 1 (LPM): all_user (4 specs) + any_fed (Risk2,3,4 only)
#     Sample: DV=1 vs pure non-users only (non_user == 1)
# -----------------------------------------------------------------------------

df_all_acute    <- df_acute %>% filter(all_user == 1 | non_user == 1)
df_anyfed_acute <- df_acute %>% filter(any_fed  == 1 | non_user == 1)

# Run models
m_all_lpm  <- run_4spec(df_all_acute, "all_user", "lpm")
m_any_lpm  <- run_r234_only(df_anyfed_acute, "any_fed", "lpm")

# Combine into 5 columns: (1)-(4) all_user, (5) any_fed risk2,3,4
models_t1 <- c(m_all_lpm, list("(5) AnyFed: Risk2,3,4" = m_any_lpm))

# N rows for each column's own sample (so Ns are correct by column)
n_all  <- add_n_rows(df_all_acute, "all_user", n_models = 4)
n_any  <- add_n_rows(df_anyfed_acute, "any_fed",  n_models = 1)

# Build a 5-column add_rows object aligned to model columns:
# - first 4 columns use all_user sample Ns
# - 5th column uses any_fed sample Ns
nrows_t1 <- data.frame(
  term = n_all$term,
  "(1)" = n_all[["(1)"]],
  "(2)" = n_all[["(2)"]],
  "(3)" = n_all[["(3)"]],
  "(4)" = n_all[["(4)"]],
  "(5)" = n_any[["(1)"]]
)

modelsummary(
  models_t1,
  stars    = c('*' = 0.1, '**' = 0.05, '***' = 0.01),
  coef_map = COEF_MAP,
  gof_map  = gof_lpm,
  add_rows = nrows_t1,
  title    = "Table 1: Acute Period (LPM) — All Users (4 specs) + Any Fed (Risk2,3,4 only)",
  notes    = "Sample: DV=1 users vs DV=0 pure non-users. Variables z-standardized. Robust SEs.",
  output   = "kableExtra"
) %>%
  kable_styling(bootstrap_options = c("striped", "hover"), font_size = 9)

Table 1: Acute Period (LPM) — All Users (4 specs) + Any Fed (Risk2,3,4 only)

	 (1) Base	 (2) +Risk1	 (3) +Risk1,2	 (4) Risk2,3,4	 (5) AnyFed: Risk2,3,4
MTM Loss (z)	0.029***
	(0.007)
Uninsured Lev (z)	0.018**
	(0.008)
MTM × Uninsured	0.013**
	(0.006)
Risk 1: \(<\) Med MTM & \(<\) Med Uninsured		−0.028*	−0.041***
		(0.014)	(0.016)
Risk 2: \(<\) Med MTM & \(>\) Med Uninsured			−0.036**	0.011	0.001
			(0.016)	(0.017)	(0.016)
Risk 3: \(>\) Med MTM & \(<\) Med Uninsured				0.019	0.019
				(0.017)	(0.016)
Risk 4: \(>\) Med MTM & \(>\) Med Uninsured				0.076***	0.079***
				(0.020)	(0.019)
Log(Assets) (z)	0.097***	0.101***	0.103***	0.098***	0.104***
	(0.008)	(0.008)	(0.008)	(0.008)	(0.008)
Cash Ratio (z)	−0.034***	−0.042***	−0.037***	−0.038***	−0.023***
	(0.006)	(0.006)	(0.006)	(0.006)	(0.006)
Loan/Deposit (z)	0.014**	0.007	0.010	0.013*	−0.005
	(0.007)	(0.007)	(0.007)	(0.007)	(0.006)
Book Equity (z)	−0.008	−0.014**	−0.013**	−0.010*	−0.015***
	(0.006)	(0.005)	(0.006)	(0.006)	(0.005)
Wholesale (z)	0.025***	0.025***	0.024***	0.025***	0.029***
	(0.007)	(0.007)	(0.007)	(0.007)	(0.007)
ROA (z)	−0.005	−0.006	−0.004	−0.006	−0.002
	(0.006)	(0.006)	(0.006)	(0.006)	(0.005)
Num.Obs.	4282	4288	4282	4282	4036
R2	0.105	0.101	0.102	0.104	0.115
R2 Adj.	0.103	0.100	0.100	0.102	0.113
N (all_user=1)	1007.000	1007.000	1007.000	1007.000	761.000
N (Sample)	4292.000	4292.000	4292.000	4292.000	4046.000
* p < 0.1, ** p < 0.05, *** p < 0.01
Sample: DV=1 users vs DV=0 pure non-users. Variables z-standardized. Robust SEs.

# SAVE TABLE 1
save_reg_table(models_t1, "Table_1_Acute_AllUsers_LPM",
               title_text = "Table 1: Acute Period (LPM) - All Users + Any Fed",
               notes_text = "Sample: DV=1 users vs DV=0 pure non-users. Variables z-standardized. Robust SEs.",
               coef_map_use = COEF_MAP, gof_map_use = gof_lpm, add_rows_use = nrows_t1)

Saved: Table_1_Acute_AllUsers_LPM (HTML + LaTeX)

4.1.2 Table 2 (Logit) - Acute Period

# Run models
m_all_logit <- run_4spec(df_all_acute, "all_user", "logit")
m_any_logit <- run_r234_only(df_anyfed_acute, "any_fed", "logit")

models_t2 <- c(m_all_logit, list("(5) AnyFed: Risk2,3,4" = m_any_logit))

# Ns are the same construction as Table 1 (by sample/column)
nrows_t2 <- nrows_t1

modelsummary(
  models_t2,
  stars    = c('*' = 0.1, '**' = 0.05, '***' = 0.01),
  coef_map = COEF_MAP,
  gof_map  = gof_logit,
  add_rows = nrows_t2,
  title    = "Table 2: Acute Period (Logit) — All Users (4 specs) + Any Fed (Risk2,3,4 only)",
  notes    = "Sample: DV=1 users vs DV=0 pure non-users. Variables z-standardized. Robust SEs.",
  output   = "kableExtra"
) %>%
  kable_styling(bootstrap_options = c("striped", "hover"), font_size = 9)

Table 2: Acute Period (Logit) — All Users (4 specs) + Any Fed (Risk2,3,4 only)

	 (1) Base	 (2) +Risk1	 (3) +Risk1,2	 (4) Risk2,3,4	 (5) AnyFed: Risk2,3,4
MTM Loss (z)	0.184***
	(0.045)
Uninsured Lev (z)	0.127***
	(0.045)
MTM × Uninsured	0.014
	(0.039)
Risk 1: \(<\) Med MTM & \(<\) Med Uninsured		−0.278**	−0.341***
		(0.110)	(0.115)
Risk 2: \(<\) Med MTM & \(>\) Med Uninsured			−0.167*	0.198	0.226
			(0.097)	(0.126)	(0.150)
Risk 3: \(>\) Med MTM & \(<\) Med Uninsured				0.218*	0.319**
				(0.126)	(0.150)
Risk 4: \(>\) Med MTM & \(>\) Med Uninsured				0.492***	0.620***
				(0.127)	(0.149)
Log(Assets) (z)	0.522***	0.549***	0.559***	0.533***	0.663***
	(0.046)	(0.042)	(0.043)	(0.044)	(0.049)
Cash Ratio (z)	−0.349***	−0.396***	−0.366***	−0.369***	−0.285***
	(0.066)	(0.062)	(0.064)	(0.064)	(0.071)
Loan/Deposit (z)	0.143***	0.105**	0.118**	0.132***	0.036
	(0.049)	(0.048)	(0.048)	(0.048)	(0.054)
Book Equity (z)	−0.137***	−0.177***	−0.168***	−0.152***	−0.261***
	(0.050)	(0.049)	(0.049)	(0.050)	(0.062)
Wholesale (z)	0.131***	0.124***	0.122***	0.124***	0.164***
	(0.036)	(0.036)	(0.036)	(0.036)	(0.040)
ROA (z)	−0.006	−0.017	−0.010	−0.023	0.018
	(0.043)	(0.042)	(0.042)	(0.043)	(0.049)
Num.Obs.	4282	4288	4282	4282	4036
N (all_user=1)	1007.000	1007.000	1007.000	1007.000	761.000
N (Sample)	4292.000	4292.000	4292.000	4292.000	4046.000
* p < 0.1, ** p < 0.05, *** p < 0.01
Sample: DV=1 users vs DV=0 pure non-users. Variables z-standardized. Robust SEs.

# SAVE TABLE 2
save_reg_table(models_t2, "Table_2_Acute_AllUsers_Logit",
               title_text = "Table 2: Acute Period (Logit) - All Users + Any Fed",
               notes_text = "Sample: DV=1 users vs DV=0 pure non-users. Variables z-standardized. Robust SEs.",
               coef_map_use = COEF_MAP, gof_map_use = gof_logit, add_rows_use = nrows_t2)

Saved: Table_2_Acute_AllUsers_Logit (HTML + LaTeX)

4.1.3 TABLE 3 (LPM): BTFP only (4 specs) — BTFP=1 vs pure non-users- Acute Period

# -----------------------------------------------------------------------------
# 3) TABLE 3 (LPM): BTFP only (4 specs) — BTFP=1 vs pure non-users
# -----------------------------------------------------------------------------

df_btfp_acute <- df_acute %>% filter(btfp_acute == 1 | non_user == 1)

m_btfp_lpm <- run_4spec(df_btfp_acute, "btfp_acute", "lpm")
n_btfp_lpm <- add_n_rows(df_btfp_acute, "btfp_acute", n_models = 4)

modelsummary(
  m_btfp_lpm,
  stars    = c('*' = 0.1, '**' = 0.05, '***' = 0.01),
  coef_map = COEF_MAP,
  gof_map  = gof_lpm,
  add_rows = n_btfp_lpm,
  title    = "Table 3: BTFP Usage (LPM) — Acute Period (Mar 13–May 1, 2023)",
  notes    = "Sample: BTFP users vs pure non-users. Variables z-standardized. Robust SEs.",
  output   = "kableExtra"
) %>%
  kable_styling(bootstrap_options = c("striped", "hover"), font_size = 9)

Table 3: BTFP Usage (LPM) — Acute Period (Mar 13–May 1, 2023)

	 (1) Base	 (2) +Risk1	 (3) +Risk1,2	 (4) Risk2,3,4
MTM Loss (z)	0.024***
	(0.006)
Uninsured Lev (z)	0.020***
	(0.007)
MTM × Uninsured	0.018***
	(0.005)
Risk 1: \(<\) Med MTM & \(<\) Med Uninsured		−0.011	−0.025**
		(0.011)	(0.012)
Risk 2: \(<\) Med MTM & \(>\) Med Uninsured			−0.036***	−0.006
			(0.014)	(0.014)
Risk 3: \(>\) Med MTM & \(<\) Med Uninsured				0.002
				(0.013)
Risk 4: \(>\) Med MTM & \(>\) Med Uninsured				0.062***
				(0.017)
Log(Assets) (z)	0.067***	0.074***	0.076***	0.071***
	(0.007)	(0.007)	(0.007)	(0.007)
Cash Ratio (z)	−0.024***	−0.031***	−0.026***	−0.027***
	(0.005)	(0.004)	(0.005)	(0.005)
Loan/Deposit (z)	−0.008	−0.016***	−0.014**	−0.011*
	(0.005)	(0.006)	(0.006)	(0.006)
Book Equity (z)	−0.011**	−0.017***	−0.017***	−0.014***
	(0.004)	(0.004)	(0.004)	(0.004)
Wholesale (z)	0.025***	0.024***	0.024***	0.025***
	(0.006)	(0.006)	(0.006)	(0.006)
ROA (z)	−0.005	−0.005	−0.003	−0.005
	(0.005)	(0.005)	(0.005)	(0.005)
Num.Obs.	3737	3743	3737	3737
R2	0.090	0.084	0.085	0.089
R2 Adj.	0.088	0.082	0.083	0.087
N (btfp_acute=1)	462.000	462.000	462.000	462.000
N (Sample)	3747.000	3747.000	3747.000	3747.000
* p < 0.1, ** p < 0.05, *** p < 0.01
Sample: BTFP users vs pure non-users. Variables z-standardized. Robust SEs.

# SAVE TABLE 3
save_reg_table(m_btfp_lpm, "Table_3_BTFP_Acute_LPM",
               title_text = "Table 3: BTFP Usage (LPM) - Acute Period",
               notes_text = "Sample: BTFP users vs pure non-users. Variables z-standardized. Robust SEs.",
               coef_map_use = COEF_MAP, gof_map_use = gof_lpm, add_rows_use = n_btfp_lpm)

Saved: Table_3_BTFP_Acute_LPM (HTML + LaTeX)

4.1.4 TABLE 4 (LOGIT): BTFP only (4 specs) - Acute Period

# -----------------------------------------------------------------------------
# 4) TABLE 4 (LOGIT): BTFP only (4 specs)
# -----------------------------------------------------------------------------

m_btfp_logit <- run_4spec(df_btfp_acute, "btfp_acute", "logit")
n_btfp_logit <- add_n_rows(df_btfp_acute, "btfp_acute", n_models = 4)

modelsummary(
  m_btfp_logit,
  stars    = c('*' = 0.1, '**' = 0.05, '***' = 0.01),
  coef_map = COEF_MAP,
  gof_map  = gof_logit,
  add_rows = n_btfp_logit,
  title    = "Table 4: BTFP Usage (Logit) — Acute Period (Mar 13–May 1, 2023)",
  notes    = "Sample: BTFP users vs pure non-users. Variables z-standardized. Robust SEs.",
  output   = "kableExtra"
) %>%
  kable_styling(bootstrap_options = c("striped", "hover"), font_size = 9)

Table 4: BTFP Usage (Logit) — Acute Period (Mar 13–May 1, 2023)

	 (1) Base	 (2) +Risk1	 (3) +Risk1,2	 (4) Risk2,3,4
MTM Loss (z)	0.214***
	(0.062)
Uninsured Lev (z)	0.209***
	(0.063)
MTM × Uninsured	0.046
	(0.057)
Risk 1: \(<\) Med MTM & \(<\) Med Uninsured		−0.301*	−0.379**
		(0.163)	(0.168)
Risk 2: \(<\) Med MTM & \(>\) Med Uninsured			−0.223	0.191
			(0.136)	(0.188)
Risk 3: \(>\) Med MTM & \(<\) Med Uninsured				0.175
				(0.184)
Risk 4: \(>\) Med MTM & \(>\) Med Uninsured				0.591***
				(0.180)
Log(Assets) (z)	0.574***	0.636***	0.649***	0.611***
	(0.062)	(0.056)	(0.057)	(0.059)
Cash Ratio (z)	−0.517***	−0.563***	−0.521***	−0.527***
	(0.105)	(0.101)	(0.103)	(0.104)
Loan/Deposit (z)	0.001	−0.061	−0.043	−0.018
	(0.066)	(0.065)	(0.065)	(0.066)
Book Equity (z)	−0.329***	−0.398***	−0.382***	−0.359***
	(0.082)	(0.081)	(0.081)	(0.082)
Wholesale (z)	0.183***	0.168***	0.166***	0.173***
	(0.046)	(0.046)	(0.046)	(0.046)
ROA (z)	−0.016	−0.018	−0.009	−0.030
	(0.063)	(0.062)	(0.062)	(0.063)
Num.Obs.	3737	3743	3737	3737
N (btfp_acute=1)	462.000	462.000	462.000	462.000
N (Sample)	3747.000	3747.000	3747.000	3747.000
* p < 0.1, ** p < 0.05, *** p < 0.01
Sample: BTFP users vs pure non-users. Variables z-standardized. Robust SEs.

# SAVE TABLE 4
save_reg_table(m_btfp_logit, "Table_4_BTFP_Acute_Logit",
               title_text = "Table 4: BTFP Usage (Logit) - Acute Period",
               notes_text = "Sample: BTFP users vs pure non-users. Variables z-standardized. Robust SEs.",
               coef_map_use = COEF_MAP, gof_map_use = gof_logit, add_rows_use = n_btfp_logit)

Saved: Table_4_BTFP_Acute_Logit (HTML + LaTeX)

4.1.5 TABLE 5 (LPM): DW only (4 specs) — DW=1 vs pure non-users - Acute Period

# -----------------------------------------------------------------------------
# 5) TABLE 5 (LPM): DW only (4 specs) — DW=1 vs pure non-users
# -----------------------------------------------------------------------------

df_dw_acute <- df_acute %>% filter(dw_acute == 1 | non_user == 1)

m_dw_lpm <- run_4spec(df_dw_acute, "dw_acute", "lpm")
n_dw_lpm <- add_n_rows(df_dw_acute, "dw_acute", n_models = 4)

modelsummary(
  m_dw_lpm,
  stars    = c('*' = 0.1, '**' = 0.05, '***' = 0.01),
  coef_map = COEF_MAP,
  gof_map  = gof_lpm,
  add_rows = n_dw_lpm,
  title    = "Table 5: Discount Window Usage (LPM) — Acute Period (Mar 13–May 1, 2023)",
  notes    = "Sample: DW users vs pure non-users. Variables z-standardized. Robust SEs.",
  output   = "kableExtra"
) %>%
  kable_styling(bootstrap_options = c("striped", "hover"), font_size = 9)

Table 5: Discount Window Usage (LPM) — Acute Period (Mar 13–May 1, 2023)

	 (1) Base	 (2) +Risk1	 (3) +Risk1,2	 (4) Risk2,3,4
MTM Loss (z)	0.016***
	(0.006)
Uninsured Lev (z)	0.013**
	(0.006)
MTM × Uninsured	0.016***
	(0.005)
Risk 1: \(<\) Med MTM & \(<\) Med Uninsured		−0.012	−0.024**
		(0.010)	(0.012)
Risk 2: \(<\) Med MTM & \(>\) Med Uninsured			−0.030**	−0.003
			(0.013)	(0.013)
Risk 3: \(>\) Med MTM & \(<\) Med Uninsured				0.008
				(0.012)
Risk 4: \(>\) Med MTM & \(>\) Med Uninsured				0.050***
				(0.016)
Log(Assets) (z)	0.084***	0.087***	0.089***	0.086***
	(0.007)	(0.007)	(0.007)	(0.007)
Cash Ratio (z)	−0.002	−0.007	−0.003	−0.003
	(0.005)	(0.005)	(0.005)	(0.005)
Loan/Deposit (z)	−0.001	−0.007	−0.004	−0.002
	(0.005)	(0.005)	(0.005)	(0.005)
Book Equity (z)	−0.000	−0.004	−0.003	−0.002
	(0.004)	(0.004)	(0.004)	(0.004)
Wholesale (z)	0.020***	0.021***	0.020***	0.020***
	(0.006)	(0.006)	(0.006)	(0.006)
ROA (z)	−0.002	−0.001	−0.000	−0.002
	(0.005)	(0.004)	(0.005)	(0.005)
Num.Obs.	3668	3674	3668	3668
R2	0.096	0.092	0.093	0.095
R2 Adj.	0.093	0.090	0.091	0.093
N (dw_acute=1)	393.000	393.000	393.000	393.000
N (Sample)	3678.000	3678.000	3678.000	3678.000
* p < 0.1, ** p < 0.05, *** p < 0.01
Sample: DW users vs pure non-users. Variables z-standardized. Robust SEs.

# SAVE TABLE 5
save_reg_table(m_dw_lpm, "Table_5_DW_Acute_LPM",
               title_text = "Table 5: Discount Window Usage (LPM) - Acute Period",
               notes_text = "Sample: DW users vs pure non-users. Variables z-standardized. Robust SEs.",
               coef_map_use = COEF_MAP, gof_map_use = gof_lpm, add_rows_use = n_dw_lpm)

Saved: Table_5_DW_Acute_LPM (HTML + LaTeX)

4.1.6 TABLE 6 (LOGIT): DW only (4 specs) - Acute Period

# -----------------------------------------------------------------------------
# 6) TABLE 6 (LOGIT): DW only (4 specs)
# -----------------------------------------------------------------------------

m_dw_logit <- run_4spec(df_dw_acute, "dw_acute", "logit")
n_dw_logit <- add_n_rows(df_dw_acute, "dw_acute", n_models = 4)

modelsummary(
  m_dw_logit,
  stars    = c('*' = 0.1, '**' = 0.05, '***' = 0.01),
  coef_map = COEF_MAP,
  gof_map  = gof_logit,
  add_rows = n_dw_logit,
  title    = "Table 6: Discount Window Usage (Logit) — Acute Period (Mar 13–May 1, 2023)",
  notes    = "Sample: DW users vs pure non-users. Variables z-standardized. Robust SEs.",
  output   = "kableExtra"
) %>%
  kable_styling(bootstrap_options = c("striped", "hover"), font_size = 9)

Table 6: Discount Window Usage (Logit) — Acute Period (Mar 13–May 1, 2023)

	 (1) Base	 (2) +Risk1	 (3) +Risk1,2	 (4) Risk2,3,4
MTM Loss (z)	0.194***
	(0.068)
Uninsured Lev (z)	0.160***
	(0.062)
MTM × Uninsured	0.083
	(0.053)
Risk 1: \(<\) Med MTM & \(<\) Med Uninsured		−0.442**	−0.567***
		(0.181)	(0.192)
Risk 2: \(<\) Med MTM & \(>\) Med Uninsured			−0.300**	0.291
			(0.140)	(0.201)
Risk 3: \(>\) Med MTM & \(<\) Med Uninsured				0.389*
				(0.211)
Risk 4: \(>\) Med MTM & \(>\) Med Uninsured				0.730***
				(0.203)
Log(Assets) (z)	0.776***	0.800***	0.821***	0.791***
	(0.061)	(0.057)	(0.059)	(0.060)
Cash Ratio (z)	−0.088	−0.140*	−0.085	−0.088
	(0.089)	(0.083)	(0.085)	(0.086)
Loan/Deposit (z)	0.091	0.044	0.071	0.085
	(0.073)	(0.071)	(0.072)	(0.072)
Book Equity (z)	−0.167**	−0.219***	−0.203***	−0.186**
	(0.080)	(0.078)	(0.078)	(0.078)
Wholesale (z)	0.172***	0.168***	0.166***	0.169***
	(0.053)	(0.053)	(0.053)	(0.053)
ROA (z)	0.038	0.038	0.048	0.034
	(0.066)	(0.065)	(0.065)	(0.065)
Num.Obs.	3668	3674	3668	3668
N (dw_acute=1)	393.000	393.000	393.000	393.000
N (Sample)	3678.000	3678.000	3678.000	3678.000
* p < 0.1, ** p < 0.05, *** p < 0.01
Sample: DW users vs pure non-users. Variables z-standardized. Robust SEs.

# SAVE TABLE 6
save_reg_table(m_dw_logit, "Table_6_DW_Acute_Logit",
               title_text = "Table 6: Discount Window Usage (Logit) - Acute Period",
               notes_text = "Sample: DW users vs pure non-users. Variables z-standardized. Robust SEs.",
               coef_map_use = COEF_MAP, gof_map_use = gof_logit, add_rows_use = n_dw_logit)

Saved: Table_6_DW_Acute_Logit (HTML + LaTeX)

5 SECTION 4: TEMPORAL ANALYSIS

5.1 4.1 BTFP Temporal Analysis (Acute → Post → Arbitrage → Wind-down)

# ==============================================================================
# BTFP ACROSS PERIODS
# Acute, Post-Acute: 2022Q4 baseline
# Arbitrage: 2023Q3 baseline
# Wind-down: 2023Q4 baseline
# ==============================================================================

# Function to run base + risk specifications
run_temporal_pair <- function(data, dv_var, family_type = "lpm") {
  
  f_base <- build_formula(dv_var, "mtm_total + uninsured_lev + mtm_x_uninsured")
  f_risk <- build_formula(dv_var, "run_risk_2 + run_risk_3 + run_risk_4")
  
  if (family_type == "lpm") {
    m_base <- feols(f_base, data = data, vcov = "hetero")
    m_risk <- feols(f_risk, data = data, vcov = "hetero")
  } else {
    m_base <- feglm(f_base, data = data, family = binomial("logit"), vcov = "hetero")
    m_risk <- feglm(f_risk, data = data, family = binomial("logit"), vcov = "hetero")
  }
  
  list(Base = m_base, Risk = m_risk)
}

5.1.1 Table 4A: BTFP Temporal (LPM)

# Prepare samples
df_btfp_acute_s <- df_acute %>% filter(btfp_acute == 1 | non_user == 1)
df_btfp_post_s <- df_post %>% filter(btfp_post == 1 | non_user == 1)
df_btfp_arb_s <- df_arb %>% filter(btfp_arb == 1 | non_user == 1)
df_btfp_wind_s <- df_wind %>% filter(btfp_wind == 1 | non_user == 1)

# Run models
m_acute <- run_temporal_pair(df_btfp_acute_s, "btfp_acute", "lpm")
m_post <- run_temporal_pair(df_btfp_post_s, "btfp_post", "lpm")
m_arb <- run_temporal_pair(df_btfp_arb_s, "btfp_arb", "lpm")
m_wind <- run_temporal_pair(df_btfp_wind_s, "btfp_wind", "lpm")

models_btfp_temp <- list(
  "Acute (Base)" = m_acute$Base, "Acute (Risk)" = m_acute$Risk,
  "Post (Base)" = m_post$Base, "Post (Risk)" = m_post$Risk,
  "Arb (Base)" = m_arb$Base, "Arb (Risk)" = m_arb$Risk,
  "Wind (Base)" = m_wind$Base, "Wind (Risk)" = m_wind$Risk
)

# N rows
n_rows_temp <- data.frame(
  term = c("N (BTFP=1)", "N (Sample)", "Baseline"),
  `Acute (Base)` = c(sum(df_btfp_acute_s$btfp_acute), nrow(df_btfp_acute_s), "2022Q4"),
  `Acute (Risk)` = c(sum(df_btfp_acute_s$btfp_acute), nrow(df_btfp_acute_s), "2022Q4"),
  `Post (Base)` = c(sum(df_btfp_post_s$btfp_post), nrow(df_btfp_post_s), "2022Q4"),
  `Post (Risk)` = c(sum(df_btfp_post_s$btfp_post), nrow(df_btfp_post_s), "2022Q4"),
  `Arb (Base)` = c(sum(df_btfp_arb_s$btfp_arb), nrow(df_btfp_arb_s), "2023Q3"),
  `Arb (Risk)` = c(sum(df_btfp_arb_s$btfp_arb), nrow(df_btfp_arb_s), "2023Q3"),
  `Wind (Base)` = c(sum(df_btfp_wind_s$btfp_wind), nrow(df_btfp_wind_s), "2023Q4"),
  `Wind (Risk)` = c(sum(df_btfp_wind_s$btfp_wind), nrow(df_btfp_wind_s), "2023Q4"),
  check.names = FALSE
)

modelsummary(
  models_btfp_temp,
  stars = c('*' = 0.1, '**' = 0.05, '***' = 0.01),
  coef_map = COEF_MAP,
  gof_map = gof_lpm,
  add_rows = n_rows_temp,
  title = "Table 4A: BTFP Temporal Analysis (LPM)",
  notes = "Arbitrage uses 2023Q3 baseline; Wind-down uses 2023Q4 baseline.",
  output = "kableExtra"
) %>%
  kable_styling(bootstrap_options = c("striped", "hover"), font_size = 8) %>%
  add_header_above(c(" " = 1, "Acute" = 2, "Post-Acute" = 2, "Arbitrage" = 2, "Wind-down" = 2))

Table 4A: BTFP Temporal Analysis (LPM)

	Acute		Post-Acute		Arbitrage		Wind-down
	Acute (Base)	Acute (Risk)	Post (Base)	Post (Risk)	Arb (Base)	Arb (Risk)	Wind (Base)	Wind (Risk)
MTM Loss (z)	0.024***		0.025***		0.022***		0.002
	(0.006)		(0.008)		(0.007)		(0.004)
Uninsured Lev (z)	0.020***		0.018**		0.014**		0.007*
	(0.007)		(0.008)		(0.007)		(0.004)
MTM × Uninsured	0.018***		0.016***		0.008		−0.001
	(0.005)		(0.006)		(0.005)		(0.003)
Risk 2: \(<\) Med MTM & \(>\) Med Uninsured		−0.006		−0.007		0.029*		0.029***
		(0.014)		(0.019)		(0.016)		(0.010)
Risk 3: \(>\) Med MTM & \(<\) Med Uninsured		0.002		0.007		0.040**		0.019*
		(0.013)		(0.019)		(0.016)		(0.010)
Risk 4: \(>\) Med MTM & \(>\) Med Uninsured		0.062***		0.055**		0.041**		0.032***
		(0.017)		(0.022)		(0.018)		(0.011)
Log(Assets) (z)	0.067***	0.071***	0.068***	0.072***	0.054***	0.055***	−0.000	−0.001
	(0.007)	(0.007)	(0.009)	(0.009)	(0.007)	(0.007)	(0.004)	(0.004)
Cash Ratio (z)	−0.024***	−0.027***	−0.052***	−0.055***	−0.036***	−0.040***	−0.010***	−0.009***
	(0.005)	(0.005)	(0.007)	(0.006)	(0.006)	(0.006)	(0.004)	(0.004)
Loan/Deposit (z)	−0.008	−0.011*	−0.021***	−0.023***	−0.003	−0.008	−0.010**	−0.009**
	(0.005)	(0.006)	(0.007)	(0.007)	(0.006)	(0.006)	(0.004)	(0.004)
Book Equity (z)	−0.011**	−0.014***	−0.020***	−0.023***	−0.006	−0.010*	−0.006*	−0.006*
	(0.004)	(0.004)	(0.006)	(0.006)	(0.006)	(0.005)	(0.003)	(0.003)
Wholesale (z)	0.025***	0.025***	0.016**	0.015**	0.101***	0.101***	0.053***	0.053***
	(0.006)	(0.006)	(0.007)	(0.007)	(0.007)	(0.007)	(0.006)	(0.006)
ROA (z)	−0.005	−0.005	−0.013**	−0.013**	−0.024***	−0.025***	−0.007*	−0.006*
	(0.005)	(0.005)	(0.006)	(0.006)	(0.006)	(0.006)	(0.004)	(0.004)
Num.Obs.	3737	3737	3515	3515	4038	4038	4043	4043
R2	0.090	0.089	0.080	0.079	0.142	0.141	0.063	0.064
R2 Adj.	0.088	0.087	0.078	0.077	0.141	0.139	0.061	0.062
N (BTFP=1)	462	462	775	775	766	766	229	229
N (Sample)	3747	3747	3525	3525	4055	4055	4061	4061
Baseline	2022Q4	2022Q4	2022Q4	2022Q4	2023Q3	2023Q3	2023Q4	2023Q4
* p < 0.1, ** p < 0.05, *** p < 0.01
Arbitrage uses 2023Q3 baseline; Wind-down uses 2023Q4 baseline.

# SAVE TABLE 4A
save_reg_table(models_btfp_temp, "Table_4A_BTFP_Temporal_LPM",
               title_text = "Table 4A: BTFP Temporal Analysis (LPM)",
               notes_text = "Arbitrage uses 2023Q3 baseline; Wind-down uses 2023Q4 baseline.",
               coef_map_use = COEF_MAP, gof_map_use = gof_lpm, add_rows_use = n_rows_temp)

Saved: Table_4A_BTFP_Temporal_LPM (HTML + LaTeX)

5.1.2 Table 4B: BTFP Arbitrage Period - All Specifications

models_btfp_arb <- run_4spec_models(df_btfp_arb_s, "btfp_arb", "lpm")
n_rows_arb <- create_n_rows(df_btfp_arb_s, "btfp_arb")

modelsummary(
  models_btfp_arb,
  stars = c('*' = 0.1, '**' = 0.05, '***' = 0.01),
  coef_map = COEF_MAP,
  gof_map = gof_lpm,
  add_rows = n_rows_arb,
  title = "Table 4B: BTFP Arbitrage Period - All Specifications (LPM, 2023Q3 Baseline)",
  output = "kableExtra"
) %>%
  kable_styling(bootstrap_options = c("striped", "hover"), font_size = 9)

Table 4B: BTFP Arbitrage Period - All Specifications (LPM, 2023Q3 Baseline)

	 (1) Base	 (2) +Risk1	 (3) +Risk1,2	 (4) Risk2,3,4
MTM Loss (z)	0.022***
	(0.007)
Uninsured Lev (z)	0.014**
	(0.007)
MTM × Uninsured	0.008
	(0.005)
Risk 1: \(<\) Med MTM & \(<\) Med Uninsured		−0.035***	−0.040***
		(0.013)	(0.014)
Risk 2: \(<\) Med MTM & \(>\) Med Uninsured			−0.011	0.029*
			(0.016)	(0.016)
Risk 3: \(>\) Med MTM & \(<\) Med Uninsured				0.040**
				(0.016)
Risk 4: \(>\) Med MTM & \(>\) Med Uninsured				0.041**
				(0.018)
Log(Assets) (z)	0.054***	0.054***	0.055***	0.055***
	(0.007)	(0.007)	(0.007)	(0.007)
Cash Ratio (z)	−0.036***	−0.042***	−0.040***	−0.040***
	(0.006)	(0.005)	(0.006)	(0.006)
Loan/Deposit (z)	−0.003	−0.009	−0.008	−0.008
	(0.006)	(0.006)	(0.006)	(0.006)
Book Equity (z)	−0.006	−0.011**	−0.010*	−0.010*
	(0.006)	(0.005)	(0.005)	(0.005)
Wholesale (z)	0.101***	0.100***	0.101***	0.101***
	(0.007)	(0.007)	(0.007)	(0.007)
ROA (z)	−0.024***	−0.025***	−0.025***	−0.025***
	(0.006)	(0.005)	(0.006)	(0.006)
Num.Obs.	4038	4048	4038	4038
R2	0.142	0.141	0.141	0.141
R2 Adj.	0.141	0.139	0.139	0.139
N (btfp_arb=1)	766.000	766.000	766.000	766.000
N (Sample)	4055.000	4055.000	4055.000	4055.000
* p < 0.1, ** p < 0.05, *** p < 0.01

# SAVE TABLE 4B
save_reg_table(models_btfp_arb, "Table_4B_BTFP_Arbitrage_LPM",
               title_text = "Table 4B: BTFP Arbitrage Period - All Specifications (LPM)",
               notes_text = "2023Q3 Baseline. Sample: BTFP users vs pure non-users.",
               coef_map_use = COEF_MAP, gof_map_use = gof_lpm, add_rows_use = n_rows_arb)

Saved: Table_4B_BTFP_Arbitrage_LPM (HTML + LaTeX)

5.1.3 Table 4A: BTFP Temporal (LOGIT)

# ------------------------------------------------------------------------------
# Helper: run Base + Risk (Risk2,3,4) specs for LOGIT (kept generic)
# ------------------------------------------------------------------------------
run_temporal_pair <- function(data, dv_var, family_type = c("lpm","logit")) {
  family_type <- match.arg(family_type)

  f_base <- build_formula(dv_var, "mtm_total + uninsured_lev + mtm_x_uninsured")
  f_risk <- build_formula(dv_var, "run_risk_2 + run_risk_3 + run_risk_4")

  if (family_type == "lpm") {
    m_base <- feols(f_base, data = data, vcov = "hetero")
    m_risk <- feols(f_risk, data = data, vcov = "hetero")
  } else {
    m_base <- feglm(f_base, data = data, family = binomial("logit"), vcov = "hetero")
    m_risk <- feglm(f_risk, data = data, family = binomial("logit"), vcov = "hetero")
  }

  list(Base = m_base, Risk = m_risk)
}

# ==============================================================================
# Table 4A (LOGIT): BTFP Temporal Analysis
# ==============================================================================
# Prepare samples (BTFP=1 vs pure non-users)
df_btfp_acute_s <- df_acute %>% filter(btfp_acute == 1 | non_user == 1)
df_btfp_post_s  <- df_post  %>% filter(btfp_post  == 1 | non_user == 1)
df_btfp_arb_s   <- df_arb   %>% filter(btfp_arb   == 1 | non_user == 1)
df_btfp_wind_s  <- df_wind  %>% filter(btfp_wind  == 1 | non_user == 1)

# Run LOGIT models
m_acute <- run_temporal_pair(df_btfp_acute_s, "btfp_acute", "logit")
m_post  <- run_temporal_pair(df_btfp_post_s,  "btfp_post",  "logit")
m_arb   <- run_temporal_pair(df_btfp_arb_s,   "btfp_arb",   "logit")
m_wind  <- run_temporal_pair(df_btfp_wind_s,  "btfp_wind",  "logit")

# Collect models (8 columns)
models_btfp_temp_logit <- list(
  "Acute (Base)" = m_acute$Base, "Acute (Risk)" = m_acute$Risk,
  "Post (Base)"  = m_post$Base,  "Post (Risk)"  = m_post$Risk,
  "Arb (Base)"   = m_arb$Base,   "Arb (Risk)"   = m_arb$Risk,
  "Wind (Base)"  = m_wind$Base,  "Wind (Risk)"  = m_wind$Risk
)

# N rows + baseline labels (same as your LPM table, but reused here)
n_rows_temp_logit <- data.frame(
  term = c("N (BTFP=1)", "N (Sample)", "Baseline"),
  `Acute (Base)` = c(sum(df_btfp_acute_s$btfp_acute, na.rm = TRUE), nrow(df_btfp_acute_s), "2022Q4"),
  `Acute (Risk)` = c(sum(df_btfp_acute_s$btfp_acute, na.rm = TRUE), nrow(df_btfp_acute_s), "2022Q4"),
  `Post (Base)`  = c(sum(df_btfp_post_s$btfp_post,  na.rm = TRUE), nrow(df_btfp_post_s),  "2022Q4"),
  `Post (Risk)`  = c(sum(df_btfp_post_s$btfp_post,  na.rm = TRUE), nrow(df_btfp_post_s),  "2022Q4"),
  `Arb (Base)`   = c(sum(df_btfp_arb_s$btfp_arb,    na.rm = TRUE), nrow(df_btfp_arb_s),   "2023Q3"),
  `Arb (Risk)`   = c(sum(df_btfp_arb_s$btfp_arb,    na.rm = TRUE), nrow(df_btfp_arb_s),   "2023Q3"),
  `Wind (Base)`  = c(sum(df_btfp_wind_s$btfp_wind,  na.rm = TRUE), nrow(df_btfp_wind_s),  "2023Q4"),
  `Wind (Risk)`  = c(sum(df_btfp_wind_s$btfp_wind,  na.rm = TRUE), nrow(df_btfp_wind_s),  "2023Q4"),
  check.names = FALSE
)

modelsummary(
  models_btfp_temp_logit,
  stars    = c('*' = 0.1, '**' = 0.05, '***' = 0.01),
  coef_map = COEF_MAP,
  gof_map  = gof_logit,
  add_rows = n_rows_temp_logit,
  title    = "Table 4C: BTFP Temporal Analysis (Logit)",
  notes    = "Arbitrage uses 2023Q3 baseline; Wind-down uses 2023Q4 baseline.",
  output   = "kableExtra"
) %>%
  kable_styling(bootstrap_options = c("striped", "hover"), font_size = 8) %>%
  add_header_above(c(" " = 1, "Acute" = 2, "Post-Acute" = 2, "Arbitrage" = 2, "Wind-down" = 2))

Table 4C: BTFP Temporal Analysis (Logit)

	Acute		Post-Acute		Arbitrage		Wind-down
	Acute (Base)	Acute (Risk)	Post (Base)	Post (Risk)	Arb (Base)	Arb (Risk)	Wind (Base)	Wind (Risk)
MTM Loss (z)	0.214***		0.126**		0.165***		0.042
	(0.062)		(0.049)		(0.053)		(0.090)
Uninsured Lev (z)	0.209***		0.115**		0.143***		0.185**
	(0.063)		(0.051)		(0.051)		(0.087)
MTM × Uninsured	0.046		0.030		−0.014		−0.069
	(0.057)		(0.043)		(0.044)		(0.074)
Risk 2: \(<\) Med MTM & \(>\) Med Uninsured		0.191		0.053		0.401***		0.787***
		(0.188)		(0.142)		(0.148)		(0.271)
Risk 3: \(>\) Med MTM & \(<\) Med Uninsured		0.175		0.080		0.394***		0.515**
		(0.184)		(0.133)		(0.145)		(0.257)
Risk 4: \(>\) Med MTM & \(>\) Med Uninsured		0.591***		0.286**		0.383**		0.770***
		(0.180)		(0.139)		(0.149)		(0.268)
Log(Assets) (z)	0.574***	0.611***	0.378***	0.406***	0.364***	0.387***	0.013	−0.001
	(0.062)	(0.059)	(0.052)	(0.050)	(0.049)	(0.047)	(0.084)	(0.082)
Cash Ratio (z)	−0.517***	−0.527***	−0.510***	−0.522***	−0.528***	−0.547***	−0.340***	−0.305***
	(0.105)	(0.104)	(0.070)	(0.069)	(0.078)	(0.077)	(0.105)	(0.100)
Loan/Deposit (z)	0.001	−0.018	−0.104**	−0.118**	0.031	−0.006	−0.180**	−0.159*
	(0.066)	(0.066)	(0.051)	(0.051)	(0.054)	(0.053)	(0.089)	(0.084)
Book Equity (z)	−0.329***	−0.359***	−0.216***	−0.236***	−0.169***	−0.197***	−0.208**	−0.192**
	(0.082)	(0.082)	(0.056)	(0.055)	(0.059)	(0.059)	(0.098)	(0.095)
Wholesale (z)	0.183***	0.173***	0.073*	0.067*	0.546***	0.541***	0.614***	0.617***
	(0.046)	(0.046)	(0.041)	(0.040)	(0.038)	(0.038)	(0.049)	(0.049)
ROA (z)	−0.016	−0.030	−0.084*	−0.090*	−0.223***	−0.236***	−0.172*	−0.169*
	(0.063)	(0.063)	(0.049)	(0.048)	(0.058)	(0.057)	(0.098)	(0.097)
Num.Obs.	3737	3737	3515	3515	4038	4038	4043	4043
N (BTFP=1)	462	462	775	775	766	766	229	229
N (Sample)	3747	3747	3525	3525	4055	4055	4061	4061
Baseline	2022Q4	2022Q4	2022Q4	2022Q4	2023Q3	2023Q3	2023Q4	2023Q4
* p < 0.1, ** p < 0.05, *** p < 0.01
Arbitrage uses 2023Q3 baseline; Wind-down uses 2023Q4 baseline.

# SAVE TABLE 4A
save_reg_table(models_btfp_temp_logit, "Table_4C_BTFP_Temporal_LOGIT",
               title_text = "Table 4C: BTFP Temporal Analysis (LOGIT)",
               notes_text = "Arbitrage uses 2023Q3 baseline; Wind-down uses 2023Q4 baseline.",
               coef_map_use = COEF_MAP, gof_map_use = gof_lpm, add_rows_use = n_rows_temp)

Saved: Table_4C_BTFP_Temporal_LOGIT (HTML + LaTeX)

# ==============================================================================
# Table 4B (LOGIT): BTFP Arbitrage Period — All Specifications
# ==============================================================================
# NOTE: Uses your existing run_4spec_models() and create_n_rows()

models_btfp_arb_logit <- run_4spec_models(df_btfp_arb_s, "btfp_arb", "logit")
n_rows_arb_logit <- create_n_rows(df_btfp_arb_s, "btfp_arb")

modelsummary(
  models_btfp_arb_logit,
  stars    = c('*' = 0.1, '**' = 0.05, '***' = 0.01),
  coef_map = COEF_MAP,
  gof_map  = gof_logit,
  add_rows = n_rows_arb_logit,
  title    = "Table 4D: BTFP Arbitrage Period - All Specifications (Logit, 2023Q3 Baseline)",
  output   = "kableExtra"
) %>%
  kable_styling(bootstrap_options = c("striped", "hover"), font_size = 9)

Table 4D: BTFP Arbitrage Period - All Specifications (Logit, 2023Q3 Baseline)

	 (1) Base	 (2) +Risk1	 (3) +Risk1,2	 (4) Risk2,3,4
MTM Loss (z)	0.165***
	(0.053)
Uninsured Lev (z)	0.143***
	(0.051)
MTM × Uninsured	−0.014
	(0.044)
Risk 1: \(<\) Med MTM & \(<\) Med Uninsured		−0.391***	−0.389***
		(0.128)	(0.134)
Risk 2: \(<\) Med MTM & \(>\) Med Uninsured			0.013	0.401***
			(0.116)	(0.148)
Risk 3: \(>\) Med MTM & \(<\) Med Uninsured				0.394***
				(0.145)
Risk 4: \(>\) Med MTM & \(>\) Med Uninsured				0.383**
				(0.149)
Log(Assets) (z)	0.364***	0.387***	0.386***	0.387***
	(0.049)	(0.045)	(0.046)	(0.047)
Cash Ratio (z)	−0.528***	−0.551***	−0.547***	−0.547***
	(0.078)	(0.075)	(0.077)	(0.077)
Loan/Deposit (z)	0.031	−0.004	−0.005	−0.006
	(0.054)	(0.052)	(0.053)	(0.053)
Book Equity (z)	−0.169***	−0.201***	−0.196***	−0.197***
	(0.059)	(0.058)	(0.059)	(0.059)
Wholesale (z)	0.546***	0.540***	0.542***	0.541***
	(0.038)	(0.038)	(0.038)	(0.038)
ROA (z)	−0.223***	−0.231***	−0.236***	−0.236***
	(0.058)	(0.056)	(0.057)	(0.057)
Num.Obs.	4038	4048	4038	4038
N (btfp_arb=1)	766.000	766.000	766.000	766.000
N (Sample)	4055.000	4055.000	4055.000	4055.000
* p < 0.1, ** p < 0.05, *** p < 0.01

save_reg_table(models_btfp_arb_logit, "Table_4D_BTFP_Arbitrage_LOGIT",
               title_text = "Table 4D: BTFP Arbitrage Period - All Specifications (LOGIT)",
               notes_text = "2023Q3 Baseline. Sample: BTFP users vs pure non-users.",
               coef_map_use = COEF_MAP, gof_map_use = gof_lpm, add_rows_use = n_rows_arb)

Saved: Table_4D_BTFP_Arbitrage_LOGIT (HTML + LaTeX)

5.1.4 MTM loss on BTFP eligible only

# ------------------------------------------------------------------------------
# Helper: run Base + Risk (Risk2,3,4) specs for LPM (kept generic)
# ------------------------------------------------------------------------------
run_temporal_pair <- function(data, dv_var, family_type = c("lpm","logit")) {
  family_type <- match.arg(family_type)

  # This adds mtm_btfp, uninsured_lev, AND their interaction
  f_base <- build_formula(dv_var, "mtm_btfp + uninsured_lev + mtm_btfp:uninsured_lev")
  f_risk <- build_formula(dv_var, "run_risk_2 + run_risk_3 + run_risk_4")

  if (family_type == "lpm") {
    m_base <- feols(f_base, data = data, vcov = "hetero")
    m_risk <- feols(f_risk, data = data, vcov = "hetero")
  } else {
    m_base <- feglm(f_base, data = data, family = binomial("logit"), vcov = "hetero")
    m_risk <- feglm(f_risk, data = data, family = binomial("logit"), vcov = "hetero")
  }

  list(Base = m_base, Risk = m_risk)
}

# ==============================================================================
# Table 4A (LOGIT): BTFP Temporal Analysis
# ==============================================================================
# Prepare samples (BTFP=1 vs pure non-users)
df_btfp_acute_s <- df_acute %>% filter(btfp_acute == 1 | non_user == 1)
df_btfp_post_s  <- df_post  %>% filter(btfp_post  == 1 | non_user == 1)
df_btfp_arb_s   <- df_arb   %>% filter(btfp_arb   == 1 | non_user == 1)
df_btfp_wind_s  <- df_wind  %>% filter(btfp_wind  == 1 | non_user == 1)

# Run LOGIT models
m_acute <- run_temporal_pair(df_btfp_acute_s, "btfp_acute", "lpm")
m_post  <- run_temporal_pair(df_btfp_post_s,  "btfp_post",  "lpm")
m_arb   <- run_temporal_pair(df_btfp_arb_s,   "btfp_arb",   "lpm")
m_wind  <- run_temporal_pair(df_btfp_wind_s,  "btfp_wind",  "lpm")

# Collect models (8 columns)
models_btfp_temp_lpm <- list(
  "Acute (Base)" = m_acute$Base, "Acute (Risk)" = m_acute$Risk,
  "Post (Base)"  = m_post$Base,  "Post (Risk)"  = m_post$Risk,
  "Arb (Base)"   = m_arb$Base,   "Arb (Risk)"   = m_arb$Risk,
  "Wind (Base)"  = m_wind$Base,  "Wind (Risk)"  = m_wind$Risk
)

# N rows + baseline labels (same as your LPM table, but reused here)
n_rows_temp_lpm <- data.frame(
  term = c("N (BTFP=1)", "N (Sample)", "Baseline"),
  `Acute (Base)` = c(sum(df_btfp_acute_s$btfp_acute, na.rm = TRUE), nrow(df_btfp_acute_s), "2022Q4"),
  `Acute (Risk)` = c(sum(df_btfp_acute_s$btfp_acute, na.rm = TRUE), nrow(df_btfp_acute_s), "2022Q4"),
  `Post (Base)`  = c(sum(df_btfp_post_s$btfp_post,  na.rm = TRUE), nrow(df_btfp_post_s),  "2022Q4"),
  `Post (Risk)`  = c(sum(df_btfp_post_s$btfp_post,  na.rm = TRUE), nrow(df_btfp_post_s),  "2022Q4"),
  `Arb (Base)`   = c(sum(df_btfp_arb_s$btfp_arb,    na.rm = TRUE), nrow(df_btfp_arb_s),   "2023Q3"),
  `Arb (Risk)`   = c(sum(df_btfp_arb_s$btfp_arb,    na.rm = TRUE), nrow(df_btfp_arb_s),   "2023Q3"),
  `Wind (Base)`  = c(sum(df_btfp_wind_s$btfp_wind,  na.rm = TRUE), nrow(df_btfp_wind_s),  "2023Q4"),
  `Wind (Risk)`  = c(sum(df_btfp_wind_s$btfp_wind,  na.rm = TRUE), nrow(df_btfp_wind_s),  "2023Q4"),
  check.names = FALSE
)

modelsummary(
  models_btfp_temp_lpm,
  stars    = c('*' = 0.1, '**' = 0.05, '***' = 0.01),
  coef_map = COEF_MAP,
  gof_map  = gof_logit,
  add_rows = n_rows_temp_lpm,
  title    = "Table 4E: BTFP Temporal Analysis (Logit)",
  notes    = "Arbitrage uses 2023Q3 baseline; Wind-down uses 2023Q4 baseline.",
  output   = "kableExtra"
) %>%
  kable_styling(bootstrap_options = c("striped", "hover"), font_size = 8) %>%
  add_header_above(c(" " = 1, "Acute" = 2, "Post-Acute" = 2, "Arbitrage" = 2, "Wind-down" = 2))

Table 4E: BTFP Temporal Analysis (Logit)

	Acute		Post-Acute		Arbitrage		Wind-down
	Acute (Base)	Acute (Risk)	Post (Base)	Post (Risk)	Arb (Base)	Arb (Risk)	Wind (Base)	Wind (Risk)
MTM Loss OMO (z)	0.016**		0.029***		0.007		0.010**
	(0.006)		(0.008)		(0.007)		(0.005)
Uninsured Lev (z)	0.014**		0.012		0.010		0.008*
	(0.006)		(0.008)		(0.006)		(0.004)
MTM_btfp × Uninsured Lev	0.011**		−0.001		−0.005		−0.004
	(0.005)		(0.007)		(0.005)		(0.004)
Risk 2: \(<\) Med MTM & \(>\) Med Uninsured		−0.006		−0.007		0.029*		0.029***
		(0.014)		(0.019)		(0.016)		(0.010)
Risk 3: \(>\) Med MTM & \(<\) Med Uninsured		0.002		0.007		0.040**		0.019*
		(0.013)		(0.019)		(0.016)		(0.010)
Risk 4: \(>\) Med MTM & \(>\) Med Uninsured		0.062***		0.055**		0.041**		0.032***
		(0.017)		(0.022)		(0.018)		(0.011)
Log(Assets) (z)	0.066***	0.071***	0.065***	0.072***	0.054***	0.055***	−0.001	−0.001
	(0.007)	(0.007)	(0.009)	(0.009)	(0.007)	(0.007)	(0.004)	(0.004)
Cash Ratio (z)	−0.031***	−0.027***	−0.057***	−0.055***	−0.043***	−0.040***	−0.010***	−0.009***
	(0.004)	(0.005)	(0.006)	(0.006)	(0.005)	(0.006)	(0.003)	(0.004)
Loan/Deposit (z)	−0.009	−0.011*	−0.014*	−0.023***	−0.007	−0.008	−0.006	−0.009**
	(0.006)	(0.006)	(0.008)	(0.007)	(0.007)	(0.006)	(0.004)	(0.004)
Book Equity (z)	−0.014***	−0.014***	−0.023***	−0.023***	−0.010*	−0.010*	−0.006*	−0.006*
	(0.004)	(0.004)	(0.006)	(0.006)	(0.006)	(0.005)	(0.003)	(0.003)
Wholesale (z)	0.024***	0.025***	0.013*	0.015**	0.100***	0.101***	0.052***	0.053***
	(0.006)	(0.006)	(0.007)	(0.007)	(0.008)	(0.007)	(0.006)	(0.006)
ROA (z)	−0.004	−0.005	−0.010	−0.013**	−0.027***	−0.025***	−0.006	−0.006*
	(0.005)	(0.005)	(0.006)	(0.006)	(0.006)	(0.006)	(0.004)	(0.004)
Num.Obs.	3737	3737	3515	3515	4038	4038	4043	4043
N (BTFP=1)	462	462	775	775	766	766	229	229
N (Sample)	3747	3747	3525	3525	4055	4055	4061	4061
Baseline	2022Q4	2022Q4	2022Q4	2022Q4	2023Q3	2023Q3	2023Q4	2023Q4
* p < 0.1, ** p < 0.05, *** p < 0.01
Arbitrage uses 2023Q3 baseline; Wind-down uses 2023Q4 baseline.

save_reg_table(models_btfp_temp_lpm, "Table_4E_BTFP_Temporal_LPM",
               title_text = "Table 4E: BTFP Temporal Analysis (LPM)",
               notes_text = "Arbitrage uses 2023Q3 baseline; Wind-down uses 2023Q4 baseline.",
               coef_map_use = COEF_MAP, gof_map_use = gof_lpm, add_rows_use = n_rows_temp)

Saved: Table_4E_BTFP_Temporal_LPM (HTML + LaTeX)

5.2 4.2 DW Temporal Analysis

5.2.1 Table 4C: DW March 10 & March 10-13 - All Specifications

# March 10 only
df_dw_mar10_s <- df_mar10 %>% filter(dw_mar10 == 1 | non_user == 1)
models_dw_mar10 <- run_4spec_models(df_dw_mar10_s, "dw_mar10", "lpm")

# March 10-13
df_dw_mar10_13_s <- df_mar10_13 %>% filter(dw_mar10_13 == 1 | non_user == 1)
models_dw_mar10_13 <- run_4spec_models(df_dw_mar10_13_s, "dw_mar10_13", "lpm")

# Combined table
models_dw_short <- list(
  "Mar10 (1)" = models_dw_mar10$`(1) Base`,
  "Mar10 (2)" = models_dw_mar10$`(2) +Risk1`,
  "Mar10 (3)" = models_dw_mar10$`(3) +Risk1,2`,
  "Mar10 (4)" = models_dw_mar10$`(4) Risk2,3,4`,
  "Mar10-13 (1)" = models_dw_mar10_13$`(1) Base`,
  "Mar10-13 (2)" = models_dw_mar10_13$`(2) +Risk1`,
  "Mar10-13 (3)" = models_dw_mar10_13$`(3) +Risk1,2`,
  "Mar10-13 (4)" = models_dw_mar10_13$`(4) Risk2,3,4`
)

n_rows_dw_short <- data.frame(
  term = c("N (DW=1)", "N (Sample)"),
  `Mar10 (1)` = c(sum(df_dw_mar10_s$dw_mar10), nrow(df_dw_mar10_s)),
  `Mar10 (2)` = c(sum(df_dw_mar10_s$dw_mar10), nrow(df_dw_mar10_s)),
  `Mar10 (3)` = c(sum(df_dw_mar10_s$dw_mar10), nrow(df_dw_mar10_s)),
  `Mar10 (4)` = c(sum(df_dw_mar10_s$dw_mar10), nrow(df_dw_mar10_s)),
  `Mar10-13 (1)` = c(sum(df_dw_mar10_13_s$dw_mar10_13), nrow(df_dw_mar10_13_s)),
  `Mar10-13 (2)` = c(sum(df_dw_mar10_13_s$dw_mar10_13), nrow(df_dw_mar10_13_s)),
  `Mar10-13 (3)` = c(sum(df_dw_mar10_13_s$dw_mar10_13), nrow(df_dw_mar10_13_s)),
  `Mar10-13 (4)` = c(sum(df_dw_mar10_13_s$dw_mar10_13), nrow(df_dw_mar10_13_s)),
  check.names = FALSE
)

modelsummary(
  models_dw_short,
  stars = c('*' = 0.1, '**' = 0.05, '***' = 0.01),
  coef_map = COEF_MAP,
  gof_map = gof_lpm,
  add_rows = n_rows_dw_short,
  title = "Table 4C: DW Usage - March 10 and March 10-13 (LPM)",
  output = "kableExtra"
) %>%
  kable_styling(bootstrap_options = c("striped", "hover"), font_size = 8) %>%
  add_header_above(c(" " = 1, "March 10 Only" = 4, "March 10-13" = 4))

Table 4C: DW Usage - March 10 and March 10-13 (LPM)

	March 10 Only				March 10-13
	 Mar10 (1)	 Mar10 (2)	 Mar10 (3)	 Mar10 (4)	 Mar10-13 (1)	 Mar10-13 (2)	 Mar10-13 (3)	 Mar10-13 (4)
MTM Loss (z)	−0.002				−0.002
	(0.002)				(0.003)
Uninsured Lev (z)	0.001				0.008**
	(0.002)				(0.003)
MTM × Uninsured	0.001				0.005*
	(0.002)				(0.003)
Risk 1: \(<\) Med MTM & \(<\) Med Uninsured		0.003	0.003			0.006	0.006
		(0.004)	(0.004)			(0.005)	(0.005)
Risk 2: \(<\) Med MTM & \(>\) Med Uninsured			−0.000	−0.003			−0.001	−0.005
			(0.005)	(0.004)			(0.007)	(0.006)
Risk 3: \(>\) Med MTM & \(<\) Med Uninsured				−0.005				−0.012**
				(0.004)				(0.005)
Risk 4: \(>\) Med MTM & \(>\) Med Uninsured				−0.001				0.004
				(0.006)				(0.008)
Log(Assets) (z)	0.013***	0.014***	0.014***	0.014***	0.026***	0.030***	0.031***	0.029***
	(0.003)	(0.003)	(0.003)	(0.003)	(0.004)	(0.004)	(0.004)	(0.004)
Cash Ratio (z)	−0.004***	−0.003***	−0.003**	−0.003**	−0.005***	−0.004**	−0.004*	−0.004**
	(0.001)	(0.001)	(0.001)	(0.001)	(0.002)	(0.002)	(0.002)	(0.002)
Loan/Deposit (z)	−0.002	−0.002	−0.002	−0.001	−0.004	−0.005*	−0.005*	−0.004
	(0.002)	(0.002)	(0.002)	(0.002)	(0.003)	(0.003)	(0.003)	(0.003)
Book Equity (z)	−0.000	−0.000	−0.000	−0.000	0.002	0.001	0.000	0.001
	(0.001)	(0.001)	(0.001)	(0.001)	(0.002)	(0.002)	(0.002)	(0.002)
Wholesale (z)	0.007***	0.007***	0.007***	0.007***	0.008***	0.008***	0.008***	0.008***
	(0.003)	(0.003)	(0.003)	(0.003)	(0.003)	(0.003)	(0.003)	(0.003)
ROA (z)	−0.004***	−0.003***	−0.003***	−0.003***	−0.003	−0.001	−0.001	−0.002
	(0.001)	(0.001)	(0.001)	(0.001)	(0.002)	(0.002)	(0.002)	(0.002)
Num.Obs.	3986	3992	3986	3986	3988	3994	3988	3988
R2	0.024	0.024	0.024	0.024	0.047	0.044	0.044	0.045
R2 Adj.	0.022	0.022	0.022	0.022	0.045	0.042	0.042	0.043
N (DW=1)	47.000	47.000	47.000	47.000	90.000	90.000	90.000	90.000
N (Sample)	3996.000	3996.000	3996.000	3996.000	3998.000	3998.000	3998.000	3998.000
* p < 0.1, ** p < 0.05, *** p < 0.01

# SAVE TABLE 4C
save_reg_table(models_dw_short, "Table_4C_DW_Mar10_LPM",
               title_text = "Table 4C: DW Usage - March 10 and March 10-13 (LPM)",
               notes_text = "Linear Probability Model estimates for Discount Window usage during the immediate crisis period. March 10 = SVB closure day (BTFP not yet available). March 10-13 = period including BTFP announcement (March 12). DW = f(MTM Losses, Uninsured Deposits, MTM × Uninsured) + controls. Sample: DW users vs pure non-users. All continuous variables z-standardized. Heteroskedasticity-robust standard errors. Significance: *** p<0.01, ** p<0.05, * p<0.10.",
               coef_map_use = COEF_MAP, gof_map_use = gof_lpm, add_rows_use = n_rows_dw_short)

Saved: Table_4C_DW_Mar10_LPM (HTML + LaTeX)

5.2.2 Table 4D: DW Temporal (Pre-BTFP → Acute)

# Pre-BTFP
df_dw_pre_s <- df_prebtfp %>% filter(dw_prebtfp == 1 | non_user == 1)
m_pre <- run_temporal_pair(df_dw_pre_s, "dw_prebtfp", "lpm")

# Acute
m_dw_acute <- run_temporal_pair(df_dw_acute, "dw_acute", "lpm")

models_dw_temp <- list(
  "Pre-BTFP (Base)" = m_pre$Base, "Pre-BTFP (Risk)" = m_pre$Risk,
  "Acute (Base)" = m_dw_acute$Base, "Acute (Risk)" = m_dw_acute$Risk
)

n_rows_dw_temp <- data.frame(
  term = c("N (DW=1)", "N (Sample)"),
  `Pre-BTFP (Base)` = c(sum(df_dw_pre_s$dw_prebtfp), nrow(df_dw_pre_s)),
  `Pre-BTFP (Risk)` = c(sum(df_dw_pre_s$dw_prebtfp), nrow(df_dw_pre_s)),
  `Acute (Base)` = c(sum(df_dw_acute$dw_acute), nrow(df_dw_acute)),
  `Acute (Risk)` = c(sum(df_dw_acute$dw_acute), nrow(df_dw_acute)),
  check.names = FALSE
)

modelsummary(
  models_dw_temp,
  stars = c('*' = 0.1, '**' = 0.05, '***' = 0.01),
  coef_map = COEF_MAP,
  gof_map = gof_lpm,
  add_rows = n_rows_dw_temp,
  title = "Table 4D: DW Temporal Analysis (LPM)",
  output = "kableExtra"
) %>%
  kable_styling(bootstrap_options = c("striped", "hover"), font_size = 9) %>%
  add_header_above(c(" " = 1, "Pre-BTFP (Mar 1-12)" = 2, "Acute (Mar 13 - May 1)" = 2))

Table 4D: DW Temporal Analysis (LPM)

	Pre-BTFP (Mar 1-12)		Acute (Mar 13 - May 1)
	Pre-BTFP (Base)	Pre-BTFP (Risk)	Acute (Base)	Acute (Risk)
MTM Loss OMO (z)	−0.000		0.008
	(0.003)		(0.006)
Uninsured Lev (z)	0.000		0.008
	(0.003)		(0.006)
MTM_btfp × Uninsured Lev	−0.000		0.009*
	(0.002)		(0.005)
Risk 2: \(<\) Med MTM & \(>\) Med Uninsured		−0.004		−0.003
		(0.007)		(0.013)
Risk 3: \(>\) Med MTM & \(<\) Med Uninsured		0.001		0.008
		(0.006)		(0.012)
Risk 4: \(>\) Med MTM & \(>\) Med Uninsured		−0.001		0.050***
		(0.008)		(0.016)
Log(Assets) (z)	0.024***	0.024***	0.084***	0.086***
	(0.004)	(0.004)	(0.007)	(0.007)
Cash Ratio (z)	−0.005**	−0.004*	−0.007	−0.003
	(0.002)	(0.002)	(0.005)	(0.005)
Loan/Deposit (z)	−0.001	−0.001	−0.003	−0.002
	(0.003)	(0.003)	(0.006)	(0.005)
Book Equity (z)	0.001	0.001	−0.003	−0.002
	(0.002)	(0.002)	(0.004)	(0.004)
Wholesale (z)	0.013***	0.013***	0.020***	0.020***
	(0.004)	(0.004)	(0.006)	(0.006)
ROA (z)	−0.002	−0.002	−0.001	−0.002
	(0.002)	(0.002)	(0.005)	(0.005)
Num.Obs.	3988	3988	3668	3668
R2	0.035	0.036	0.093	0.095
R2 Adj.	0.033	0.033	0.091	0.093
N (DW=1)	100.000	100.000	393.000	393.000
N (Sample)	3998.000	3998.000	3678.000	3678.000
* p < 0.1, ** p < 0.05, *** p < 0.01

# SAVE TABLE 4D
save_reg_table(models_dw_temp, "Table_4D_DW_Temporal_LPM",
               title_text = "Table 4D: DW Temporal Analysis (LPM)",
               notes_text = "Linear Probability Model estimates comparing DW usage before and after BTFP introduction. Pre-BTFP (Mar 1-12, 2023): only DW available. Acute (Mar 13 - May 1, 2023): both DW and BTFP available. DW = f(MTM Losses, Uninsured Deposits, MTM × Uninsured) + controls. Sample: DW users vs pure non-users. All continuous variables z-standardized. Heteroskedasticity-robust standard errors. Significance: *** p<0.01, ** p<0.05, * p<0.10.",
               coef_map_use = COEF_MAP, gof_map_use = gof_lpm, add_rows_use = n_rows_dw_temp)

Saved: Table_4D_DW_Temporal_LPM (HTML + LaTeX)

6 SECTION 5: SIZE ANALYSIS

# ==============================================================================
# FUNCTION: Generate Size Analysis Tables
# ==============================================================================

generate_size_tables <- function(data, dv_var, dv_label) {
  
  sizes <- c("Small (<$1B)", "Medium ($1B-$100B)", "Large (>$100B)")
  short_names <- c("Small", "Medium", "Large")
  
  models_lpm_base <- list()
  models_lpm_risk <- list()
  models_log_base <- list()
  models_log_risk <- list()
  
  n_info <- data.frame(term = c("N (Users)", "N (Sample)"))
  
  for (i in seq_along(sizes)) {
    sz <- sizes[i]
    nm <- short_names[i]
    
    df_sub <- data %>% filter(size_cat == sz)
    n_users <- sum(df_sub[[dv_var]] == 1, na.rm = TRUE)
    n_total <- nrow(df_sub)
    
    f_base <- build_formula(dv_var, "mtm_total + uninsured_lev + mtm_x_uninsured")
    f_r1 <- build_formula(dv_var, "run_risk_1")
    f_r12 <- build_formula(dv_var, "run_risk_1 + run_risk_2")
    f_r234 <- build_formula(dv_var, "run_risk_2 + run_risk_3 + run_risk_4")
    
    if (n_users >= 5) {
      # LPM
      models_lpm_base[[paste0(nm, " (Base)")]] <- feols(f_base, data = df_sub, vcov = "hetero")
      models_lpm_risk[[paste0(nm, " (+R1)")]] <- feols(f_r1, data = df_sub, vcov = "hetero")
      models_lpm_risk[[paste0(nm, " (+R1,2)")]] <- feols(f_r12, data = df_sub, vcov = "hetero")
      models_lpm_risk[[paste0(nm, " (R2,3,4)")]] <- feols(f_r234, data = df_sub, vcov = "hetero")
      
      n_info[[paste0(nm, " (Base)")]] <- c(n_users, n_total)
      n_info[[paste0(nm, " (+R1)")]] <- c(n_users, n_total)
      n_info[[paste0(nm, " (+R1,2)")]] <- c(n_users, n_total)
      n_info[[paste0(nm, " (R2,3,4)")]] <- c(n_users, n_total)
    }
  }
  
  # Combine base and risk models
  all_models <- c(models_lpm_base, models_lpm_risk)
  
  list(models = all_models, n_rows = n_info)
}

6.1 5.1 Small Banks - Acute Period

# Small banks sample
df_small_btfp <- df_acute %>% filter(size_cat == "Small (<$1B)") %>% filter(btfp_acute == 1 | non_user == 1)
df_small_dw <- df_acute %>% filter(size_cat == "Small (<$1B)") %>% filter(dw_acute == 1 | non_user == 1)
df_small_anyfed <- df_acute %>% filter(size_cat == "Small (<$1B)") %>% filter(any_fed == 1 | non_user == 1)

# BTFP
if (sum(df_small_btfp$btfp_acute) >= 5) {
  models_small_btfp <- run_4spec_models(df_small_btfp, "btfp_acute", "lpm")
  n_small_btfp <- create_n_rows(df_small_btfp, "btfp_acute")
  
  modelsummary(
    models_small_btfp,
    stars = c('*' = 0.1, '**' = 0.05, '***' = 0.01),
    coef_map = COEF_MAP, gof_map = gof_lpm, add_rows = n_small_btfp,
    title = "Table 5A: BTFP Usage - Small Banks (<$1B) - Acute Period (LPM)",
    output = "kableExtra"
  ) %>% kable_styling(bootstrap_options = c("striped", "hover"), font_size = 9)
  
  # SAVE TABLE 5A
  save_reg_table(models_small_btfp, "Table_5A_BTFP_Small_LPM",
                 title_text = "Table 5A: BTFP Usage - Small Banks (<$1B) - Acute Period (LPM)",
                 notes_text = "Linear Probability Model estimates for BTFP usage among small banks (total assets < $1B). BTFP = f(MTM Losses, Uninsured Deposits, MTM × Uninsured) + controls. Sample: BTFP users vs pure non-users within size category. All continuous variables z-standardized. Heteroskedasticity-robust standard errors. Significance: *** p<0.01, ** p<0.05, * p<0.10.",
                 coef_map_use = COEF_MAP, gof_map_use = gof_lpm, add_rows_use = n_small_btfp)
}

Saved: Table_5A_BTFP_Small_LPM (HTML + LaTeX)

# DW
if (sum(df_small_dw$dw_acute) >= 5) {
  models_small_dw <- run_4spec_models(df_small_dw, "dw_acute", "lpm")
  n_small_dw <- create_n_rows(df_small_dw, "dw_acute")
  
  modelsummary(
    models_small_dw,
    stars = c('*' = 0.1, '**' = 0.05, '***' = 0.01),
    coef_map = COEF_MAP, gof_map = gof_lpm, add_rows = n_small_dw,
    title = "Table 5B: DW Usage - Small Banks (<$1B) - Acute Period (LPM)",
    output = "kableExtra"
  ) %>% kable_styling(bootstrap_options = c("striped", "hover"), font_size = 9)
  
  # SAVE TABLE 5B
  save_reg_table(models_small_dw, "Table_5B_DW_Small_LPM",
                 title_text = "Table 5B: DW Usage - Small Banks (<$1B) - Acute Period (LPM)",
                 notes_text = "Linear Probability Model estimates for DW usage among small banks. DW = f(MTM Losses, Uninsured Deposits, MTM × Uninsured) + controls. Sample: DW users vs pure non-users within size category. Heteroskedasticity-robust standard errors. Significance: *** p<0.01, ** p<0.05, * p<0.10.",
                 coef_map_use = COEF_MAP, gof_map_use = gof_lpm, add_rows_use = n_small_dw)
}

Saved: Table_5B_DW_Small_LPM (HTML + LaTeX)

# Any Fed
if (sum(df_small_anyfed$any_fed) >= 5) {
  models_small_anyfed <- run_4spec_models(df_small_anyfed, "any_fed", "lpm")
  n_small_anyfed <- create_n_rows(df_small_anyfed, "any_fed")
  
  modelsummary(
    models_small_anyfed,
    stars = c('*' = 0.1, '**' = 0.05, '***' = 0.01),
    coef_map = COEF_MAP, gof_map = gof_lpm, add_rows = n_small_anyfed,
    title = "Table 5C: Any Fed Usage - Small Banks (<$1B) - Acute Period (LPM)",
    output = "kableExtra"
  ) %>% kable_styling(bootstrap_options = c("striped", "hover"), font_size = 9)
  
  # SAVE TABLE 5C
  save_reg_table(models_small_anyfed, "Table_5C_AnyFed_Small_LPM",
                 title_text = "Table 5C: Any Fed Usage - Small Banks (<$1B) - Acute Period (LPM)",
                 notes_text = "Linear Probability Model estimates for any Fed facility usage (BTFP or DW) among small banks. Any Fed = f(MTM Losses, Uninsured Deposits, MTM × Uninsured) + controls. Sample: Any Fed users vs pure non-users within size category. Heteroskedasticity-robust standard errors. Significance: *** p<0.01, ** p<0.05, * p<0.10.",
                 coef_map_use = COEF_MAP, gof_map_use = gof_lpm, add_rows_use = n_small_anyfed)
}

Saved: Table_5C_AnyFed_Small_LPM (HTML + LaTeX)

6.2 5.2 Medium Banks - Acute Period

# Medium banks sample
df_med_btfp <- df_acute %>% filter(size_cat == "Medium ($1B-$100B)") %>% filter(btfp_acute == 1 | non_user == 1)
df_med_dw <- df_acute %>% filter(size_cat == "Medium ($1B-$100B)") %>% filter(dw_acute == 1 | non_user == 1)
df_med_anyfed <- df_acute %>% filter(size_cat == "Medium ($1B-$100B)") %>% filter(any_fed == 1 | non_user == 1)

# BTFP
if (sum(df_med_btfp$btfp_acute) >= 5) {
  models_med_btfp <- run_4spec_models(df_med_btfp, "btfp_acute", "lpm")
  n_med_btfp <- create_n_rows(df_med_btfp, "btfp_acute")
  
  modelsummary(
    models_med_btfp,
    stars = c('*' = 0.1, '**' = 0.05, '***' = 0.01),
    coef_map = COEF_MAP, gof_map = gof_lpm, add_rows = n_med_btfp,
    title = "Table 5D: BTFP Usage - Medium Banks ($1B-$100B) - Acute Period (LPM)",
    output = "kableExtra"
  ) %>% kable_styling(bootstrap_options = c("striped", "hover"), font_size = 9)
  
  # SAVE TABLE 5D
  save_reg_table(models_med_btfp, "Table_5D_BTFP_Medium_LPM",
                 title_text = "Table 5D: BTFP Usage - Medium Banks ($1B-$100B) - Acute Period (LPM)",
                 notes_text = "Linear Probability Model estimates for BTFP usage among medium banks ($1B-$100B total assets). BTFP = f(MTM Losses, Uninsured Deposits, MTM × Uninsured) + controls. Sample: BTFP users vs pure non-users within size category. All continuous variables z-standardized. Heteroskedasticity-robust standard errors. Significance: *** p<0.01, ** p<0.05, * p<0.10.",
                 coef_map_use = COEF_MAP, gof_map_use = gof_lpm, add_rows_use = n_med_btfp)
}

Saved: Table_5D_BTFP_Medium_LPM (HTML + LaTeX)

# DW
if (sum(df_med_dw$dw_acute) >= 5) {
  models_med_dw <- run_4spec_models(df_med_dw, "dw_acute", "lpm")
  n_med_dw <- create_n_rows(df_med_dw, "dw_acute")
  
  modelsummary(
    models_med_dw,
    stars = c('*' = 0.1, '**' = 0.05, '***' = 0.01),
    coef_map = COEF_MAP, gof_map = gof_lpm, add_rows = n_med_dw,
    title = "Table 5E: DW Usage - Medium Banks ($1B-$100B) - Acute Period (LPM)",
    output = "kableExtra"
  ) %>% kable_styling(bootstrap_options = c("striped", "hover"), font_size = 9)
  
  # SAVE TABLE 5E
  save_reg_table(models_med_dw, "Table_5E_DW_Medium_LPM",
                 title_text = "Table 5E: DW Usage - Medium Banks ($1B-$100B) - Acute Period (LPM)",
                 notes_text = "Linear Probability Model estimates for DW usage among medium banks. DW = f(MTM Losses, Uninsured Deposits, MTM × Uninsured) + controls. Sample: DW users vs pure non-users within size category. Heteroskedasticity-robust standard errors. Significance: *** p<0.01, ** p<0.05, * p<0.10.",
                 coef_map_use = COEF_MAP, gof_map_use = gof_lpm, add_rows_use = n_med_dw)
}

Saved: Table_5E_DW_Medium_LPM (HTML + LaTeX)

# Any Fed
if (sum(df_med_anyfed$any_fed) >= 5) {
  models_med_anyfed <- run_4spec_models(df_med_anyfed, "any_fed", "lpm")
  n_med_anyfed <- create_n_rows(df_med_anyfed, "any_fed")
  
  modelsummary(
    models_med_anyfed,
    stars = c('*' = 0.1, '**' = 0.05, '***' = 0.01),
    coef_map = COEF_MAP, gof_map = gof_lpm, add_rows = n_med_anyfed,
    title = "Table 5F: Any Fed Usage - Medium Banks ($1B-$100B) - Acute Period (LPM)",
    output = "kableExtra"
  ) %>% kable_styling(bootstrap_options = c("striped", "hover"), font_size = 9)
  
  # SAVE TABLE 5F
  save_reg_table(models_med_anyfed, "Table_5F_AnyFed_Medium_LPM",
                 title_text = "Table 5F: Any Fed Usage - Medium Banks ($1B-$100B) - Acute Period (LPM)",
                 notes_text = "Linear Probability Model estimates for any Fed facility usage (BTFP or DW) among medium banks. Any Fed = f(MTM Losses, Uninsured Deposits, MTM × Uninsured) + controls. Sample: Any Fed users vs pure non-users within size category. Heteroskedasticity-robust standard errors. Significance: *** p<0.01, ** p<0.05, * p<0.10.",
                 coef_map_use = COEF_MAP, gof_map_use = gof_lpm, add_rows_use = n_med_anyfed)
}

Saved: Table_5F_AnyFed_Medium_LPM (HTML + LaTeX)

6.3 5.3 Large Banks - Acute Period

# Large banks sample
df_large_btfp <- df_acute %>% filter(size_cat == "Large (>$100B)") %>% filter(btfp_acute == 1 | non_user == 1)
df_large_dw <- df_acute %>% filter(size_cat == "Large (>$100B)") %>% filter(dw_acute == 1 | non_user == 1)
df_large_anyfed <- df_acute %>% filter(size_cat == "Large (>$100B)") %>% filter(any_fed == 1 | non_user == 1)

# BTFP
if (sum(df_large_btfp$btfp_acute) >= 5) {
  models_large_btfp <- run_4spec_models(df_large_btfp, "btfp_acute", "lpm")
  n_large_btfp <- create_n_rows(df_large_btfp, "btfp_acute")
  
  modelsummary(
    models_large_btfp,
    stars = c('*' = 0.1, '**' = 0.05, '***' = 0.01),
    coef_map = COEF_MAP, gof_map = gof_lpm, add_rows = n_large_btfp,
    title = "Table 5G: BTFP Usage - Large Banks (>$100B) - Acute Period (LPM)",
    output = "kableExtra"
  ) %>% kable_styling(bootstrap_options = c("striped", "hover"), font_size = 9)
  
  # SAVE TABLE 5G
  save_reg_table(models_large_btfp, "Table_5G_BTFP_Large_LPM",
                 title_text = "Table 5G: BTFP Usage - Large Banks (>$100B) - Acute Period (LPM)",
                 notes_text = "Linear Probability Model estimates for BTFP usage among large banks (total assets > $100B). BTFP = f(MTM Losses, Uninsured Deposits, MTM × Uninsured) + controls. Sample: BTFP users vs pure non-users within size category. All continuous variables z-standardized. Heteroskedasticity-robust standard errors. Significance: *** p<0.01, ** p<0.05, * p<0.10.",
                 coef_map_use = COEF_MAP, gof_map_use = gof_lpm, add_rows_use = n_large_btfp)
}

Saved: Table_5G_BTFP_Large_LPM (HTML + LaTeX)

# DW
if (sum(df_large_dw$dw_acute) >= 5) {
  models_large_dw <- run_4spec_models(df_large_dw, "dw_acute", "lpm")
  n_large_dw <- create_n_rows(df_large_dw, "dw_acute")
  
  modelsummary(
    models_large_dw,
    stars = c('*' = 0.1, '**' = 0.05, '***' = 0.01),
    coef_map = COEF_MAP, gof_map = gof_lpm, add_rows = n_large_dw,
    title = "Table 5H: DW Usage - Large Banks (>$100B) - Acute Period (LPM)",
    output = "kableExtra"
  ) %>% kable_styling(bootstrap_options = c("striped", "hover"), font_size = 9)
  
  # SAVE TABLE 5H
  save_reg_table(models_large_dw, "Table_5H_DW_Large_LPM",
                 title_text = "Table 5H: DW Usage - Large Banks (>$100B) - Acute Period (LPM)",
                 notes_text = "Linear Probability Model estimates for DW usage among large banks. DW = f(MTM Losses, Uninsured Deposits, MTM × Uninsured) + controls. Sample: DW users vs pure non-users within size category. Heteroskedasticity-robust standard errors. Significance: *** p<0.01, ** p<0.05, * p<0.10.",
                 coef_map_use = COEF_MAP, gof_map_use = gof_lpm, add_rows_use = n_large_dw)
}

Saved: Table_5H_DW_Large_LPM (HTML + LaTeX)

# Any Fed
if (sum(df_large_anyfed$any_fed) >= 5) {
  models_large_anyfed <- run_4spec_models(df_large_anyfed, "any_fed", "lpm")
  n_large_anyfed <- create_n_rows(df_large_anyfed, "any_fed")
  
  modelsummary(
    models_large_anyfed,
    stars = c('*' = 0.1, '**' = 0.05, '***' = 0.01),
    coef_map = COEF_MAP, gof_map = gof_lpm, add_rows = n_large_anyfed,
    title = "Table 5I: Any Fed Usage - Large Banks (>$100B) - Acute Period (LPM)",
    output = "kableExtra"
  ) %>% kable_styling(bootstrap_options = c("striped", "hover"), font_size = 9)
  
  # SAVE TABLE 5I
  save_reg_table(models_large_anyfed, "Table_5I_AnyFed_Large_LPM",
                 title_text = "Table 5I: Any Fed Usage - Large Banks (>$100B) - Acute Period (LPM)",
                 notes_text = "Linear Probability Model estimates for any Fed facility usage (BTFP or DW) among large banks. Any Fed = f(MTM Losses, Uninsured Deposits, MTM × Uninsured) + controls. Sample: Any Fed users vs pure non-users within size category. Heteroskedasticity-robust standard errors. Significance: *** p<0.01, ** p<0.05, * p<0.10.",
                 coef_map_use = COEF_MAP, gof_map_use = gof_lpm, add_rows_use = n_large_anyfed)
}

Saved: Table_5I_AnyFed_Large_LPM (HTML + LaTeX)

6.4 5.4 Consolidated Size Comparison

# ==============================================================================
# CONSOLIDATED: Compare BTFP across size categories
# ==============================================================================

models_size_btfp <- list()
n_rows_size <- data.frame(term = c("N (BTFP=1)", "N (Sample)"))

for (sz in c("Small (<$1B)", "Medium ($1B-$100B)", "Large (>$100B)")) {
  df_sub <- df_acute %>% filter(size_cat == sz) %>% filter(btfp_acute == 1 | non_user == 1)
  n_users <- sum(df_sub$btfp_acute)
  n_total <- nrow(df_sub)
  
  short_nm <- case_when(
    sz == "Small (<$1B)" ~ "Small",
    sz == "Medium ($1B-$100B)" ~ "Medium",
    TRUE ~ "Large"
  )
  
  if (n_users >= 5) {
    f_base <- build_formula("btfp_acute", "mtm_total + uninsured_lev + mtm_x_uninsured")
    f_risk <- build_formula("btfp_acute", "run_risk_2 + run_risk_3 + run_risk_4")
    
    models_size_btfp[[paste0(short_nm, " (Base)")]] <- feols(f_base, data = df_sub, vcov = "hetero")
    models_size_btfp[[paste0(short_nm, " (Risk)")]] <- feols(f_risk, data = df_sub, vcov = "hetero")
    
    n_rows_size[[paste0(short_nm, " (Base)")]] <- c(n_users, n_total)
    n_rows_size[[paste0(short_nm, " (Risk)")]] <- c(n_users, n_total)
  }
}

if (length(models_size_btfp) > 0) {
  modelsummary(
    models_size_btfp,
    stars = c('*' = 0.1, '**' = 0.05, '***' = 0.01),
    coef_map = COEF_MAP,
    gof_map = gof_lpm,
    add_rows = n_rows_size,
    title = "Table 5J: BTFP Usage by Bank Size - Consolidated (LPM)",
    output = "kableExtra"
  ) %>%
    kable_styling(bootstrap_options = c("striped", "hover"), font_size = 8) %>%
    add_header_above(c(" " = 1, "Small" = 2, "Medium" = 2, "Large" = 2))
  
  # SAVE TABLE 5J
  save_reg_table(models_size_btfp, "Table_5J_BTFP_Size_Consolidated_LPM",
                 title_text = "Table 5J: BTFP Usage by Bank Size - Consolidated (LPM)",
                 notes_text = "Linear Probability Model estimates comparing BTFP usage across bank size categories. Small: assets < $1B; Medium: $1B - $100B; Large: > $100B. Base specification: continuous MTM + Uninsured + interaction. Risk specification: 2×2 risk category dummies. BTFP = f(Risk Categories) + controls. Sample: BTFP users vs pure non-users within each size group. All continuous variables z-standardized. Heteroskedasticity-robust standard errors. Significance: *** p<0.01, ** p<0.05, * p<0.10.",
                 coef_map_use = COEF_MAP, gof_map_use = gof_lpm, add_rows_use = n_rows_size)
}

Saved: Table_5J_BTFP_Size_Consolidated_LPM (HTML + LaTeX)

7 SECTION 6: PLOTS

7.1 6.1 Coefficient Comparison Plot

# ==============================================================================
# PLOT 1: Coefficient Comparison Across Facilities
# ==============================================================================

# Extract coefficients from risk model
extract_risk_coefs <- function(model, facility_name) {
  coefs <- coef(model)
  ses <- sqrt(diag(vcov(model)))
  
  risk_vars <- c("run_risk_2", "run_risk_3", "run_risk_4")
  risk_labels <- c("Risk 2:\nLow MTM, High Unins", 
                   "Risk 3:\nHigh MTM, Low Unins", 
                   "Risk 4:\nHigh MTM, High Unins")
  
  tibble(
    Facility = facility_name,
    Variable = risk_labels,
    Coefficient = coefs[risk_vars],
    SE = ses[risk_vars],
    CI_low = Coefficient - 1.96 * SE,
    CI_high = Coefficient + 1.96 * SE
  )
}

# Get coefficients from BTFP and DW models
coef_btfp <- extract_risk_coefs(m_btfp_lpm$`(4) Risk2,3,4`, "BTFP")
coef_dw <- extract_risk_coefs(m_dw_lpm$`(4) Risk2,3,4`, "DW")
coef_all <- bind_rows(coef_btfp, coef_dw)

# Plot
p_coef <- ggplot(coef_all, aes(x = Variable, y = Coefficient, color = Facility)) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray50") +
  geom_pointrange(aes(ymin = CI_low, ymax = CI_high), 
                  position = position_dodge(width = 0.4), size = 0.8) +
  scale_color_manual(values = c("BTFP" = "#2166AC", "DW" = "#B2182B")) +
  labs(
    title = "Effect of Bank Risk Category on Emergency Borrowing",
    subtitle = "LPM Coefficients with 95% CI (Reference: Risk 1 - Low MTM, Low Uninsured)",
    x = NULL,
    y = "Change in Probability of Borrowing",
    color = "Facility"
  ) +
  theme_minimal(base_size = 12) +
  theme(
    legend.position = "bottom",
    plot.title = element_text(face = "bold"),
    panel.grid.minor = element_blank()
  )

print(p_coef)

# SAVE FIGURE 1
save_figure(p_coef, "Figure_1_Coefficient_Comparison", width = 12, height = 8)

## Saved: Figure_1_Coefficient_Comparison (PDF + PNG)

7.2 6.2 High-Risk Quadrant Heatmap

# ==============================================================================
# PLOT 2: High-Risk Quadrant Interaction Heatmap
# ==============================================================================

# Fit continuous interaction model
fml_int <- as.formula(paste("btfp_acute ~ mtm_total * uninsured_lev +", CONTROLS))
m_cont <- feols(fml_int, data = df_btfp_acute, vcov = "hetero")

# Create prediction grid (hold controls at 0 = mean)
grid <- expand.grid(
  mtm_total = seq(min(df_btfp_acute$mtm_total, na.rm = TRUE), 
                  max(df_btfp_acute$mtm_total, na.rm = TRUE), length.out = 50),
  uninsured_lev = seq(min(df_btfp_acute$uninsured_lev, na.rm = TRUE), 
                      max(df_btfp_acute$uninsured_lev, na.rm = TRUE), length.out = 50),
  ln_assets = 0, cash_ratio = 0, loan_to_deposit = 0, 
  book_equity_ratio = 0, wholesale = 0, roa = 0
)

grid$prob <- predict(m_cont, newdata = grid)
grid$prob <- pmax(pmin(grid$prob, 1), 0)  # Bound probabilities

p_heat <- ggplot(grid, aes(x = mtm_total, y = uninsured_lev, fill = prob)) +
  geom_tile() +
  scale_fill_viridis_c(option = "magma", labels = scales::percent, limits = c(0, 0.5)) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "white", alpha = 0.5) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "white", alpha = 0.5) +
  annotate("text", x = 2, y = 2, label = "Death\nZone", color = "white", fontface = "bold", size = 5) +
  labs(
    title = "Interaction of MTM Losses and Uninsured Deposits on BTFP Usage",
    subtitle = "Predicted Probability of BTFP Borrowing (Acute Period)",
    x = "MTM Losses (Z-Score)",
    y = "Uninsured Deposits (Z-Score)",
    fill = "Pr(BTFP)"
  ) +
  theme_minimal(base_size = 12) +
  theme(
    panel.grid = element_blank(),
    plot.title = element_text(face = "bold")
  )

print(p_heat)

# SAVE FIGURE 2
save_figure(p_heat, "Figure_2_DeathZone_Heatmap", width = 10, height = 8)

## Saved: Figure_2_DeathZone_Heatmap (PDF + PNG)

7.3 6.3 Temporal Evolution Plot

# ==============================================================================
# PLOT 3: Temporal Evolution of Risk Category Effects
# ==============================================================================

# Extract Risk 4 coefficient across periods
temporal_coefs <- tibble(
  Period = c("Acute", "Post-Acute", "Arbitrage", "Wind-down"),
  Period_Order = 1:4,
  Coefficient = c(
    coef(m_acute$Risk)["run_risk_4"],
    coef(m_post$Risk)["run_risk_4"],
    coef(m_arb$Risk)["run_risk_4"],
    coef(m_wind$Risk)["run_risk_4"]
  ),
  SE = c(
    sqrt(vcov(m_acute$Risk)["run_risk_4", "run_risk_4"]),
    sqrt(vcov(m_post$Risk)["run_risk_4", "run_risk_4"]),
    sqrt(vcov(m_arb$Risk)["run_risk_4", "run_risk_4"]),
    sqrt(vcov(m_wind$Risk)["run_risk_4", "run_risk_4"])
  )
) %>%
  mutate(
    CI_low = Coefficient - 1.96 * SE,
    CI_high = Coefficient + 1.96 * SE,
    Period = factor(Period, levels = c("Acute", "Post-Acute", "Arbitrage", "Wind-down"))
  )

p_temporal <- ggplot(temporal_coefs, aes(x = Period, y = Coefficient)) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray50") +
  geom_pointrange(aes(ymin = CI_low, ymax = CI_high), 
                  color = "#2166AC", size = 1) +
  geom_line(aes(group = 1), color = "#2166AC", alpha = 0.5) +
  labs(
    title = "Evolution of High-Risk Quadrant Effect on BTFP Borrowing",
    subtitle = "Risk 4 (High MTM Loss + High Uninsured) Coefficient Over Time",
    x = "Crisis Period",
    y = "LPM Coefficient (vs. Reference: Low MTM, Low Unins)"
  ) +
  theme_minimal(base_size = 12) +
  theme(
    plot.title = element_text(face = "bold"),
    panel.grid.minor = element_blank()
  )

print(p_temporal)

# SAVE FIGURE 3
save_figure(p_temporal, "Figure_3_Temporal_Evolution", width = 12, height = 6)

## Saved: Figure_3_Temporal_Evolution (PDF + PNG)

7.4 6.4 Participation Comparison by Size

# ==============================================================================
# PLOT 4: Participation Rates by Size and Facility
# ==============================================================================

participation <- df_acute %>%
  group_by(size_cat) %>%
  summarise(
    BTFP = mean(btfp_acute, na.rm = TRUE) * 100,
    DW = mean(dw_acute, na.rm = TRUE) * 100,
    `Any Fed` = mean(any_fed, na.rm = TRUE) * 100,
    N = n(),
    .groups = "drop"
  ) %>%
  pivot_longer(cols = c(BTFP, DW, `Any Fed`), names_to = "Facility", values_to = "Rate")

p_part <- ggplot(participation, aes(x = size_cat, y = Rate, fill = Facility)) +
  geom_bar(stat = "identity", position = position_dodge(width = 0.8), width = 0.7) +
  geom_text(aes(label = sprintf("%.1f%%", Rate)), 
            position = position_dodge(width = 0.8), vjust = -0.5, size = 3) +
  scale_fill_manual(values = c("BTFP" = "#2166AC", "DW" = "#B2182B", "Any Fed" = "#7570B3")) +
  scale_y_continuous(labels = function(x) paste0(x, "%"), expand = expansion(mult = c(0, 0.15))) +
  labs(
    title = "Emergency Facility Participation by Bank Size",
    subtitle = "Acute Period (March 13 - May 1, 2023)",
    x = "Bank Size Category",
    y = "Participation Rate",
    fill = "Facility"
  ) +
  theme_minimal(base_size = 12) +
  theme(
    legend.position = "bottom",
    plot.title = element_text(face = "bold"),
    panel.grid.minor = element_blank()
  )

print(p_part)

# SAVE FIGURE 4
save_figure(p_part, "Figure_4_Participation_by_Size", width = 10, height = 6)

## Saved: Figure_4_Participation_by_Size (PDF + PNG)

7.5 6.5 Risk Category Distribution

# ==============================================================================
# PLOT 5: Risk Category Distribution by User Group
# ==============================================================================

risk_dist <- df_acute %>%
  filter(user_group != "Both") %>%  # Simplify for visualization
  group_by(user_group) %>%
  summarise(
    `Risk 1\n(Low, Low)` = mean(run_risk_1, na.rm = TRUE) * 100,
    `Risk 2\n(Low, High)` = mean(run_risk_2, na.rm = TRUE) * 100,
    `Risk 3\n(High, Low)` = mean(run_risk_3, na.rm = TRUE) * 100,
    `Risk 4\n(High, High)` = mean(run_risk_4, na.rm = TRUE) * 100,
    N = n(),
    .groups = "drop"
  ) %>%
  pivot_longer(cols = -c(user_group, N), names_to = "Risk_Category", values_to = "Pct")

p_risk <- ggplot(risk_dist, aes(x = user_group, y = Pct, fill = Risk_Category)) +
  geom_bar(stat = "identity", position = "stack") +
  scale_fill_manual(values = c(
    "Risk 1\n(Low, Low)" = "#4DAF4A",
    "Risk 2\n(Low, High)" = "#FFFF33",
    "Risk 3\n(High, Low)" = "#FF7F00",
    "Risk 4\n(High, High)" = "#E41A1C"
  )) +
  scale_y_continuous(labels = function(x) paste0(x, "%")) +
  labs(
    title = "Risk Category Composition by Facility Choice",
    subtitle = "Acute Period - Distribution of Banks Across Risk Categories",
    x = "Facility Choice",
    y = "Percentage of Banks",
    fill = "Risk Category"
  ) +
  theme_minimal(base_size = 12) +
  theme(
    legend.position = "right",
    plot.title = element_text(face = "bold")
  )

print(p_risk)

# SAVE FIGURE 5
save_figure(p_risk, "Figure_5_Risk_Distribution", width = 12, height = 6)

## Saved: Figure_5_Risk_Distribution (PDF + PNG)

8 SECTION 7: ADDITIONAL VISUALIZATIONS

9 # Run Risk = f(Fundamentals, Liquidity Mismatch, Fundamentals × Liquidity Mismatch)

9.1 7.1 Scatter Plot: MTM Losses vs Uninsured Deposits

# ==============================================================================
# PLOT 6: The Core Relationship
# X = MTM Losses (Fundamentals), Y = Uninsured Deposits (Liquidity Mismatch)
# Color = Facility Choice
# ==============================================================================

p_scatter <- df_acute %>%
  filter(user_group %in% c("Neither", "BTFP_Only", "DW_Only")) %>%
  mutate(
    user_group = factor(user_group, 
                        levels = c("Neither", "BTFP_Only", "DW_Only"),
                        labels = c("No Facility", "BTFP", "Discount Window"))
  ) %>%
  ggplot(aes(x = mtm_total, y = uninsured_lev, color = user_group)) +
  # Add quadrant lines at zero (which is median for z-scores)
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray50", alpha = 0.7) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "gray50", alpha = 0.7) +
  # Add quadrant labels
  annotate("text", x = -2, y = 2.5, label = "Risk 2\n(Low MTM, High Unins)", 
           color = "gray40", size = 3, hjust = 0) +
  annotate("text", x = 2, y = 2.5, label = "Risk 4\n(High MTM, High Unins)", 
           color = "#E41A1C", size = 3, fontface = "bold", hjust = 1) +
  annotate("text", x = -2, y = -2, label = "Risk 1\n(Low MTM, Low Unins)", 
           color = "#4DAF4A", size = 3, hjust = 0) +
  annotate("text", x = 2, y = -2, label = "Risk 3\n(High MTM, Low Unins)", 
           color = "gray40", size = 3, hjust = 1) +
  # Points
  geom_point(alpha = 0.6, size = 2) +
  scale_color_manual(values = c("No Facility" = "gray60", 
                                "BTFP" = "#2166AC", 
                                "Discount Window" = "#B2182B")) +
  labs(
    title = "The Run Risk Landscape: Where Do Facility Users Fall?",
    subtitle = "Bank Fundamentals vs Liquidity Mismatch (Z-Scores, Acute Period)",
    x = "MTM Loss (Z-Score)\n← Sound Fundamentals | Weak Fundamentals →",
    y = "Uninsured Deposits (Z-Score)\n← Low Run Risk | High Run Risk →",
    color = "Facility Used",
    caption = "Note: Dashed lines at z=0 (sample medians). Risk categories: (MTM Loss level, Uninsured Deposit level)."
  ) +
  theme_minimal(base_size = 12) +
  theme(
    legend.position = "bottom",
    plot.title = element_text(face = "bold"),
    panel.grid.minor = element_blank()
  )

print(p_scatter)

# SAVE FIGURE 6
save_figure(p_scatter, "Figure_6_Scatter_MTM_vs_Uninsured", width = 12, height = 8)

## Saved: Figure_6_Scatter_MTM_vs_Uninsured (PDF + PNG)

9.2 7.2 Marginal Effects Plot

# ==============================================================================
# PLOT 7: How Does the Effect of MTM Change with Uninsured Deposits?
# Tests the interaction: MTM × Uninsured
# ==============================================================================

# Get coefficients from base model
base_model <- feols(btfp_acute ~ mtm_total + uninsured_lev + mtm_x_uninsured + ln_assets + 
                    cash_ratio + loan_to_deposit + book_equity_ratio + wholesale + roa,
                    data = df_btfp_acute_s, vcov = "hetero")

base_coefs <- coef(base_model)
base_vcov <- vcov(base_model)

# Calculate marginal effect of MTM at different levels of uninsured
uninsured_levels <- seq(-2, 3, by = 1)
uninsured_labels <- c("-2 SD\n(Very Low)", "-1 SD\n(Low)", "0\n(Median)", 
                      "+1 SD\n(High)", "+2 SD\n(Very High)", "+3 SD\n(Extreme)")

marginal_df <- tibble(
  Uninsured_Level = factor(uninsured_labels, levels = uninsured_labels),
  Uninsured_Z = uninsured_levels,
  # Marginal effect of MTM = beta_mtm + beta_interaction * uninsured
  Marginal_Effect = base_coefs["mtm_total"] + base_coefs["mtm_x_uninsured"] * uninsured_levels
) %>%
  mutate(
    # SE via delta method
    SE = sqrt(base_vcov["mtm_total", "mtm_total"] + 
              2 * Uninsured_Z * base_vcov["mtm_total", "mtm_x_uninsured"] +
              Uninsured_Z^2 * base_vcov["mtm_x_uninsured", "mtm_x_uninsured"]),
    CI_low = Marginal_Effect - 1.96 * SE,
    CI_high = Marginal_Effect + 1.96 * SE
  )

p_marginal <- ggplot(marginal_df, aes(x = Uninsured_Level, y = Marginal_Effect)) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray50") +
  geom_pointrange(aes(ymin = CI_low, ymax = CI_high), color = "#2166AC", size = 1) +
  geom_line(aes(group = 1), color = "#2166AC", alpha = 0.5) +
  labs(
    title = "How Liquidity Risk Amplifies Fundamental Weakness",
    subtitle = "Marginal Effect of MTM Loss on BTFP Probability at Different Uninsured Deposit Levels",
    x = "Uninsured Deposit Level (Z-Score)",
    y = "Marginal Effect of 1 SD Increase in MTM Loss\non Probability of BTFP Borrowing",
    caption = "Note: Based on LPM estimates. 95% confidence intervals shown.\nPositive effect means higher MTM losses increase BTFP usage."
  ) +
  theme_minimal(base_size = 12) +
  theme(
    plot.title = element_text(face = "bold"),
    panel.grid.minor = element_blank()
  )

print(p_marginal)

# SAVE FIGURE 7
save_figure(p_marginal, "Figure_7_Marginal_Effects_Interaction", width = 12, height = 6)

## Saved: Figure_7_Marginal_Effects_Interaction (PDF + PNG)

9.3 7.3 Density Comparison: MTM Losses by User Group

# ==============================================================================
# PLOT 8: Distribution of MTM Losses by Facility Choice
# ==============================================================================

p_density_mtm <- df_acute %>%
  filter(user_group %in% c("Neither", "BTFP_Only", "DW_Only")) %>%
  mutate(user_group = factor(user_group, 
                             levels = c("Neither", "BTFP_Only", "DW_Only"),
                             labels = c("No Facility", "BTFP", "Discount Window"))) %>%
  ggplot(aes(x = mtm_total_w, fill = user_group, color = user_group)) +
  geom_density(alpha = 0.4) +
  geom_vline(xintercept = median(df_acute$mtm_total_w, na.rm = TRUE), 
             linetype = "dashed", color = "gray30") +
  annotate("text", x = median(df_acute$mtm_total_w, na.rm = TRUE) + 0.3, y = 0.15, 
           label = "Median", color = "gray30", hjust = 0) +
  scale_fill_manual(values = c("No Facility" = "gray70", "BTFP" = "#2166AC", "Discount Window" = "#B2182B")) +
  scale_color_manual(values = c("No Facility" = "gray50", "BTFP" = "#2166AC", "Discount Window" = "#B2182B")) +
  labs(
    title = "Distribution of Mark-to-Market Losses by Facility Choice",
    subtitle = "Acute Period: Comparing Fundamentals Across User Groups",
    x = "MTM Loss / Total Assets (%)",
    y = "Density",
    fill = "Facility Used",
    color = "Facility Used",
    caption = "Note: Higher values indicate larger unrealized losses (worse fundamentals)."
  ) +
  theme_minimal(base_size = 12) +
  theme(
    legend.position = "bottom",
    plot.title = element_text(face = "bold")
  )

print(p_density_mtm)

# SAVE FIGURE 8
save_figure(p_density_mtm, "Figure_8_Density_MTM_Losses", width = 12, height = 6)

## Saved: Figure_8_Density_MTM_Losses (PDF + PNG)

9.4 7.4 Density Comparison: Uninsured Deposits by User Group

# ==============================================================================
# PLOT 9: Distribution of Uninsured Deposits by Facility Choice
# ==============================================================================

p_density_unins <- df_acute %>%
  filter(user_group %in% c("Neither", "BTFP_Only", "DW_Only")) %>%
  mutate(user_group = factor(user_group, 
                             levels = c("Neither", "BTFP_Only", "DW_Only"),
                             labels = c("No Facility", "BTFP", "Discount Window"))) %>%
  ggplot(aes(x = uninsured_lev_w, fill = user_group, color = user_group)) +
  geom_density(alpha = 0.4) +
  geom_vline(xintercept = median(df_acute$uninsured_lev_w, na.rm = TRUE), 
             linetype = "dashed", color = "gray30") +
  annotate("text", x = median(df_acute$uninsured_lev_w, na.rm = TRUE) + 1, y = 0.04, 
           label = "Median", color = "gray30", hjust = 0) +
  scale_fill_manual(values = c("No Facility" = "gray70", "BTFP" = "#2166AC", "Discount Window" = "#B2182B")) +
  scale_color_manual(values = c("No Facility" = "gray50", "BTFP" = "#2166AC", "Discount Window" = "#B2182B")) +
  labs(
    title = "Distribution of Uninsured Deposits by Facility Choice",
    subtitle = "Acute Period: Comparing Liquidity Mismatch Across User Groups",
    x = "Uninsured Deposits / Total Assets (%)",
    y = "Density",
    fill = "Facility Used",
    color = "Facility Used",
    caption = "Note: Higher values indicate greater liquidity mismatch (higher run risk)."
  ) +
  theme_minimal(base_size = 12) +
  theme(
    legend.position = "bottom",
    plot.title = element_text(face = "bold")
  )

print(p_density_unins)

# SAVE FIGURE 9
save_figure(p_density_unins, "Figure_9_Density_Uninsured_Deposits", width = 12, height = 6)

## Saved: Figure_9_Density_Uninsured_Deposits (PDF + PNG)

9.5 7.5 Combined Risk Profile Panel

# ==============================================================================
# PLOT 10: Combined Panel - Key Story Visualization
# ==============================================================================

# Panel A: Average risk metrics by group
p_panel_a <- df_acute %>%
  filter(user_group %in% c("Neither", "BTFP_Only", "DW_Only")) %>%
  mutate(user_group = factor(user_group, 
                             levels = c("Neither", "BTFP_Only", "DW_Only"),
                             labels = c("No Facility", "BTFP", "DW"))) %>%
  group_by(user_group) %>%
  summarise(
    MTM_Loss = mean(mtm_total_w, na.rm = TRUE),
    MTM_SE = sd(mtm_total_w, na.rm = TRUE) / sqrt(n()),
    Uninsured = mean(uninsured_lev_w, na.rm = TRUE),
    Unins_SE = sd(uninsured_lev_w, na.rm = TRUE) / sqrt(n()),
    .groups = "drop"
  ) %>%
  pivot_longer(cols = c(MTM_Loss, Uninsured), names_to = "Metric", values_to = "Value") %>%
  mutate(
    SE = ifelse(Metric == "MTM_Loss", MTM_SE, Unins_SE),
    Metric = ifelse(Metric == "MTM_Loss", "MTM Loss (%)", "Uninsured Deposits (%)")
  ) %>%
  ggplot(aes(x = user_group, y = Value, fill = Metric)) +
  geom_col(position = position_dodge(0.8), width = 0.7) +
  geom_errorbar(aes(ymin = Value - 1.96*SE, ymax = Value + 1.96*SE),
                position = position_dodge(0.8), width = 0.2) +
  scale_fill_manual(values = c("MTM Loss (%)" = "#E41A1C", "Uninsured Deposits (%)" = "#377EB8")) +
  labs(
    title = "A. Average Risk Profile by Facility",
    x = NULL,
    y = "Mean Value (%)",
    fill = NULL
  ) +
  theme_minimal(base_size = 11) +
  theme(legend.position = "bottom")

# Panel B: High-risk quadrant share
p_panel_b <- df_acute %>%
  filter(user_group %in% c("Neither", "BTFP_Only", "DW_Only")) %>%
  group_by(user_group) %>%
  summarise(
    HighRisk = mean(run_risk_4, na.rm = TRUE) * 100,
    .groups = "drop"
  ) %>%
  mutate(user_group = factor(user_group, 
                             levels = c("Neither", "BTFP_Only", "DW_Only"),
                             labels = c("No Facility", "BTFP", "DW"))) %>%
  ggplot(aes(x = user_group, y = HighRisk, fill = user_group)) +
  geom_col() +
  geom_text(aes(label = sprintf("%.1f%%", HighRisk)), vjust = -0.5, size = 4) +
  scale_fill_manual(values = c("No Facility" = "gray70", "BTFP" = "#2166AC", "DW" = "#B2182B")) +
  scale_y_continuous(expand = expansion(mult = c(0, 0.15))) +
  labs(
    title = "B. Share in Risk 4 Quadrant",
    subtitle = "(High MTM, High Uninsured)",
    x = NULL,
    y = "Percent of Group"
  ) +
  theme_minimal(base_size = 11) +
  theme(legend.position = "none")

# Combine panels
p_combined <- (p_panel_a | p_panel_b) +
  plot_annotation(
    title = "The Run Risk Profile: How BTFP Users Differ",
    subtitle = "Acute Period (March 13 - May 1, 2023)",
    theme = theme(plot.title = element_text(face = "bold", size = 14))
  )

print(p_combined)

# SAVE FIGURE 10
save_figure(p_combined, "Figure_10_Combined_Risk_Panel", width = 14, height = 6)

## Saved: Figure_10_Combined_Risk_Panel (PDF + PNG)

# ==============================================================================
# COEFFICIENT PLOTS: Acute Period - BTFP, DW, AnyFed
# ==============================================================================

## Function to extract coefficients from a model
extract_all_coefs <- function(model, model_name) {
  coefs <- coef(model)
  ses <- sqrt(diag(vcov(model)))
  
  # Get all coefficients except intercept
  vars <- names(coefs)[names(coefs) != "(Intercept)"]
  
  tibble(
    Model = model_name,
    Variable = vars,
    Coefficient = coefs[vars],
    SE = ses[vars],
    CI_low = Coefficient - 1.96 * SE,
    CI_high = Coefficient + 1.96 * SE,
    Significant = ifelse(CI_low > 0 | CI_high < 0, "Yes", "No")
  )
}

## Variable labels for display
VAR_LABELS <- c(
  "mtm_total" = "MTM Loss (Total)",
  "uninsured_lev" = "Uninsured Deposits",
  "mtm_x_uninsured" = "MTM × Uninsured",
  "run_risk_2" = "Risk 2 (Low MTM, High Unins)",
  "run_risk_3" = "Risk 3 (High MTM, Low Unins)",
  "run_risk_4" = "Risk 4 (High MTM, High Unins)",
  "ln_assets" = "Log(Assets)",
  "cash_ratio" = "Cash Ratio",
  "loan_to_deposit" = "Loan/Deposit",
  "book_equity_ratio" = "Book Equity",
  "wholesale" = "Wholesale Funding",
  "roa" = "ROA"
)

# ------------------------------------------------------------------------------
# PLOT A1: BTFP Acute Period - All Coefficients (Base Model)
# ------------------------------------------------------------------------------

coef_btfp_base <- extract_all_coefs(m_btfp_lpm$`(1) Base`, "BTFP") %>%
  mutate(Variable = factor(Variable, levels = rev(names(VAR_LABELS)), labels = rev(VAR_LABELS)))

p_btfp_coef_base <- ggplot(coef_btfp_base, aes(x = Coefficient, y = Variable, color = Significant)) +

geom_vline(xintercept = 0, linetype = "dashed", color = "gray50") +
  geom_pointrange(aes(xmin = CI_low, xmax = CI_high), size = 0.8) +
  scale_color_manual(values = c("Yes" = "#2166AC", "No" = "gray60"), guide = "none") +
  labs(
    title = "BTFP Borrowing Determinants (Acute Period)",
    subtitle = "Base specification: MTM Loss + Uninsured + Interaction + Controls",
    x = "Coefficient Estimate (95% CI)",
    y = NULL,
    caption = "Note: LPM estimates. Blue = significant at 5% level."
  ) +
  theme_minimal(base_size = 12) +
  theme(
    plot.title = element_text(face = "bold"),
    panel.grid.major.y = element_blank()
  )

print(p_btfp_coef_base)

save_figure(p_btfp_coef_base, "Figure_A1_BTFP_Acute_Coefficients_Base", width = 10, height = 7)

## Saved: Figure_A1_BTFP_Acute_Coefficients_Base (PDF + PNG)

# ------------------------------------------------------------------------------
# PLOT A2: BTFP Acute Period - Risk Category Model
# ------------------------------------------------------------------------------

coef_btfp_risk <- extract_all_coefs(m_btfp_lpm$`(4) Risk2,3,4`, "BTFP") %>%
  mutate(Variable = factor(Variable, levels = rev(names(VAR_LABELS)), labels = rev(VAR_LABELS)))

p_btfp_coef_risk <- ggplot(coef_btfp_risk, aes(x = Coefficient, y = Variable, color = Significant)) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "gray50") +
  geom_pointrange(aes(xmin = CI_low, xmax = CI_high), size = 0.8) +
  scale_color_manual(values = c("Yes" = "#2166AC", "No" = "gray60"), guide = "none") +
  labs(
    title = "BTFP Borrowing Determinants (Acute Period)",
    subtitle = "Risk category specification (Reference: Risk 1 - Low MTM, Low Uninsured)",
    x = "Coefficient Estimate (95% CI)",
    y = NULL,
    caption = "Note: LPM estimates. Blue = significant at 5% level."
  ) +
  theme_minimal(base_size = 12) +
  theme(
    plot.title = element_text(face = "bold"),
    panel.grid.major.y = element_blank()
  )

print(p_btfp_coef_risk)

save_figure(p_btfp_coef_risk, "Figure_A2_BTFP_Acute_Coefficients_Risk", width = 10, height = 7)

## Saved: Figure_A2_BTFP_Acute_Coefficients_Risk (PDF + PNG)

# ------------------------------------------------------------------------------
# PLOT A3: DW Acute Period - All Coefficients (Base Model)
# ------------------------------------------------------------------------------

coef_dw_base <- extract_all_coefs(m_dw_lpm$`(1) Base`, "DW") %>%
  mutate(Variable = factor(Variable, levels = rev(names(VAR_LABELS)), labels = rev(VAR_LABELS)))

p_dw_coef_base <- ggplot(coef_dw_base, aes(x = Coefficient, y = Variable, color = Significant)) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "gray50") +
  geom_pointrange(aes(xmin = CI_low, xmax = CI_high), size = 0.8) +
  scale_color_manual(values = c("Yes" = "#B2182B", "No" = "gray60"), guide = "none") +
  labs(
    title = "Discount Window Borrowing Determinants (Acute Period)",
    subtitle = "Base specification: MTM Loss + Uninsured + Interaction + Controls",
    x = "Coefficient Estimate (95% CI)",
    y = NULL,
    caption = "Note: LPM estimates. Red = significant at 5% level."
  ) +
  theme_minimal(base_size = 12) +
  theme(
    plot.title = element_text(face = "bold"),
    panel.grid.major.y = element_blank()
  )

print(p_dw_coef_base)

save_figure(p_dw_coef_base, "Figure_A3_DW_Acute_Coefficients_Base", width = 10, height = 7)

## Saved: Figure_A3_DW_Acute_Coefficients_Base (PDF + PNG)

# ------------------------------------------------------------------------------
# PLOT A4: DW Acute Period - Risk Category Model
# ------------------------------------------------------------------------------

coef_dw_risk <- extract_all_coefs(m_dw_lpm$`(4) Risk2,3,4`, "DW") %>%
  mutate(Variable = factor(Variable, levels = rev(names(VAR_LABELS)), labels = rev(VAR_LABELS)))

p_dw_coef_risk <- ggplot(coef_dw_risk, aes(x = Coefficient, y = Variable, color = Significant)) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "gray50") +
  geom_pointrange(aes(xmin = CI_low, xmax = CI_high), size = 0.8) +
  scale_color_manual(values = c("Yes" = "#B2182B", "No" = "gray60"), guide = "none") +
  labs(
    title = "Discount Window Borrowing Determinants (Acute Period)",
    subtitle = "Risk category specification (Reference: Risk 1 - Low MTM, Low Uninsured)",
    x = "Coefficient Estimate (95% CI)",
    y = NULL,
    caption = "Note: LPM estimates. Red = significant at 5% level."
  ) +
  theme_minimal(base_size = 12) +
  theme(
    plot.title = element_text(face = "bold"),
    panel.grid.major.y = element_blank()
  )

print(p_dw_coef_risk)

save_figure(p_dw_coef_risk, "Figure_A4_DW_Acute_Coefficients_Risk", width = 10, height = 7)

## Saved: Figure_A4_DW_Acute_Coefficients_Risk (PDF + PNG)

# ------------------------------------------------------------------------------
# PLOT A5: AnyFed Acute Period - All Coefficients (Base Model)
# ------------------------------------------------------------------------------

coef_anyfed_base <- extract_all_coefs(m_all_lpm$`(1) Base`, "AnyFed") %>%
  mutate(Variable = factor(Variable, levels = rev(names(VAR_LABELS)), labels = rev(VAR_LABELS)))

p_anyfed_coef_base <- ggplot(coef_anyfed_base, aes(x = Coefficient, y = Variable, color = Significant)) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "gray50") +
  geom_pointrange(aes(xmin = CI_low, xmax = CI_high), size = 0.8) +
  scale_color_manual(values = c("Yes" = "#4DAF4A", "No" = "gray60"), guide = "none") +
  labs(
    title = "Any Fed Facility Borrowing Determinants (Acute Period)",
    subtitle = "Base specification: MTM Loss + Uninsured + Interaction + Controls",
    x = "Coefficient Estimate (95% CI)",
    y = NULL,
    caption = "Note: LPM estimates. Green = significant at 5% level."
  ) +
  theme_minimal(base_size = 12) +
  theme(
    plot.title = element_text(face = "bold"),
    panel.grid.major.y = element_blank()
  )

print(p_anyfed_coef_base)

save_figure(p_anyfed_coef_base, "Figure_A5_AnyFed_Acute_Coefficients_Base", width = 10, height = 7)

## Saved: Figure_A5_AnyFed_Acute_Coefficients_Base (PDF + PNG)

# ------------------------------------------------------------------------------
# PLOT A6: AnyFed Acute Period - Risk Category Model
# ------------------------------------------------------------------------------

coef_anyfed_risk <- extract_all_coefs(m_all_lpm$`(4) Risk2,3,4`, "AnyFed") %>%
  mutate(Variable = factor(Variable, levels = rev(names(VAR_LABELS)), labels = rev(VAR_LABELS)))

p_anyfed_coef_risk <- ggplot(coef_anyfed_risk, aes(x = Coefficient, y = Variable, color = Significant)) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "gray50") +
  geom_pointrange(aes(xmin = CI_low, xmax = CI_high), size = 0.8) +
  scale_color_manual(values = c("Yes" = "#4DAF4A", "No" = "gray60"), guide = "none") +
  labs(
    title = "Any Fed Facility Borrowing Determinants (Acute Period)",
    subtitle = "Risk category specification (Reference: Risk 1 - Low MTM, Low Uninsured)",
    x = "Coefficient Estimate (95% CI)",
    y = NULL,
    caption = "Note: LPM estimates. Green = significant at 5% level."
  ) +
  theme_minimal(base_size = 12) +
  theme(
    plot.title = element_text(face = "bold"),
    panel.grid.major.y = element_blank()
  )

print(p_anyfed_coef_risk)

save_figure(p_anyfed_coef_risk, "Figure_A6_AnyFed_Acute_Coefficients_Risk", width = 10, height = 7)

## Saved: Figure_A6_AnyFed_Acute_Coefficients_Risk (PDF + PNG)

# ------------------------------------------------------------------------------
# PLOT A7: Combined Acute Period - BTFP vs DW vs AnyFed (Base Model)
# ------------------------------------------------------------------------------

coef_combined_base <- bind_rows(
  extract_all_coefs(m_btfp_lpm$`(1) Base`, "BTFP"),
  extract_all_coefs(m_dw_lpm$`(1) Base`, "DW"),
  extract_all_coefs(m_all_lpm$`(1) Base`, "Any Fed")
) %>%
  mutate(
    Variable = factor(Variable, levels = rev(names(VAR_LABELS)), labels = rev(VAR_LABELS)),
    Model = factor(Model, levels = c("BTFP", "DW", "Any Fed"))
  )

p_combined_base <- ggplot(coef_combined_base, aes(x = Coefficient, y = Variable, color = Model)) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "gray50") +
  geom_pointrange(aes(xmin = CI_low, xmax = CI_high), 
                  position = position_dodge(width = 0.6), size = 0.7) +
  scale_color_manual(values = c("BTFP" = "#2166AC", "DW" = "#B2182B", "Any Fed" = "#4DAF4A")) +
  labs(
    title = "Borrowing Determinants by Facility Type (Acute Period)",
    subtitle = "Base specification: MTM Loss + Uninsured + Interaction + Controls",
    x = "Coefficient Estimate (95% CI)",
    y = NULL,
    color = "Facility",
    caption = "Note: LPM estimates with heteroskedasticity-robust standard errors."
  ) +
  theme_minimal(base_size = 12) +
  theme(
    plot.title = element_text(face = "bold"),
    panel.grid.major.y = element_blank(),
    legend.position = "bottom"
  )

print(p_combined_base)

save_figure(p_combined_base, "Figure_A7_Combined_Acute_Coefficients_Base", width = 12, height = 8)

## Saved: Figure_A7_Combined_Acute_Coefficients_Base (PDF + PNG)

# ------------------------------------------------------------------------------
# PLOT A8: Combined Acute Period - BTFP vs DW vs AnyFed (Risk Model)
# ------------------------------------------------------------------------------

coef_combined_risk <- bind_rows(
  extract_all_coefs(m_btfp_lpm$`(4) Risk2,3,4`, "BTFP"),
  extract_all_coefs(m_dw_lpm$`(4) Risk2,3,4`, "DW"),
  extract_all_coefs(m_all_lpm$`(4) Risk2,3,4`, "Any Fed")
) %>%
  mutate(
    Variable = factor(Variable, levels = rev(names(VAR_LABELS)), labels = rev(VAR_LABELS)),
    Model = factor(Model, levels = c("BTFP", "DW", "Any Fed"))
  )

p_combined_risk <- ggplot(coef_combined_risk, aes(x = Coefficient, y = Variable, color = Model)) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "gray50") +
  geom_pointrange(aes(xmin = CI_low, xmax = CI_high), 
                  position = position_dodge(width = 0.6), size = 0.7) +
  scale_color_manual(values = c("BTFP" = "#2166AC", "DW" = "#B2182B", "Any Fed" = "#4DAF4A")) +
  labs(
    title = "Borrowing Determinants by Facility Type (Acute Period)",
    subtitle = "Risk category specification (Reference: Risk 1 - Low MTM, Low Uninsured)",
    x = "Coefficient Estimate (95% CI)",
    y = NULL,
    color = "Facility",
    caption = "Note: LPM estimates with heteroskedasticity-robust standard errors."
  ) +
  theme_minimal(base_size = 12) +
  theme(
    plot.title = element_text(face = "bold"),
    panel.grid.major.y = element_blank(),
    legend.position = "bottom"
  )

print(p_combined_risk)

save_figure(p_combined_risk, "Figure_A8_Combined_Acute_Coefficients_Risk", width = 12, height = 8)

## Saved: Figure_A8_Combined_Acute_Coefficients_Risk (PDF + PNG)

# ==============================================================================
# TEMPORAL DYNAMICS PLOTS: Coefficient Evolution Across Crisis Phases
# ==============================================================================

# ------------------------------------------------------------------------------
# Function to extract temporal coefficients
# ------------------------------------------------------------------------------

extract_temporal_coefs <- function(models_list, facility_name, key_vars = c("mtm_total", "uninsured_lev", "mtm_x_uninsured")) {
  
  results <- list()
  
  for (period in names(models_list)) {
    model <- models_list[[period]]
    coefs <- coef(model)
    ses <- sqrt(diag(vcov(model)))
    
    for (v in key_vars) {
      if (v %in% names(coefs)) {
        results[[length(results) + 1]] <- tibble(
          Facility = facility_name,
          Period = period,
          Variable = v,
          Coefficient = coefs[v],
          SE = ses[v],
          CI_low = Coefficient - 1.96 * SE,
          CI_high = Coefficient + 1.96 * SE
        )
      }
    }
  }
  
  bind_rows(results)
}

# ------------------------------------------------------------------------------
# BTFP Temporal Models (need to create these first)
# ------------------------------------------------------------------------------

# Create temporal models for BTFP (Base specification)
btfp_temporal_base <- list(
  "Acute" = feols(btfp_acute ~ mtm_total + uninsured_lev + mtm_x_uninsured + 
                    ln_assets + cash_ratio + loan_to_deposit + book_equity_ratio + wholesale + roa,
                  data = df_btfp_acute_s, vcov = "hetero"),
  "Post-Acute" = feols(btfp_post ~ mtm_total + uninsured_lev + mtm_x_uninsured + 
                         ln_assets + cash_ratio + loan_to_deposit + book_equity_ratio + wholesale + roa,
                       data = df_btfp_post_s, vcov = "hetero"),
  "Arbitrage" = feols(btfp_arb ~ mtm_total + uninsured_lev + mtm_x_uninsured + 
                        ln_assets + cash_ratio + loan_to_deposit + book_equity_ratio + wholesale + roa,
                      data = df_btfp_arb_s, vcov = "hetero"),
  "Wind-down" = feols(btfp_wind ~ mtm_total + uninsured_lev + mtm_x_uninsured + 
                        ln_assets + cash_ratio + loan_to_deposit + book_equity_ratio + wholesale + roa,
                      data = df_btfp_wind_s, vcov = "hetero")
)

# Create temporal models for BTFP (Risk specification)
btfp_temporal_risk <- list(
  "Acute" = feols(btfp_acute ~ run_risk_2 + run_risk_3 + run_risk_4 + 
                    ln_assets + cash_ratio + loan_to_deposit + book_equity_ratio + wholesale + roa,
                  data = df_btfp_acute_s, vcov = "hetero"),
  "Post-Acute" = feols(btfp_post ~ run_risk_2 + run_risk_3 + run_risk_4 + 
                         ln_assets + cash_ratio + loan_to_deposit + book_equity_ratio + wholesale + roa,
                       data = df_btfp_post_s, vcov = "hetero"),
  "Arbitrage" = feols(btfp_arb ~ run_risk_2 + run_risk_3 + run_risk_4 + 
                        ln_assets + cash_ratio + loan_to_deposit + book_equity_ratio + wholesale + roa,
                      data = df_btfp_arb_s, vcov = "hetero"),
  "Wind-down" = feols(btfp_wind ~ run_risk_2 + run_risk_3 + run_risk_4 + 
                        ln_assets + cash_ratio + loan_to_deposit + book_equity_ratio + wholesale + roa,
                      data = df_btfp_wind_s, vcov = "hetero")
)

# ------------------------------------------------------------------------------
# DW Temporal Models
# ------------------------------------------------------------------------------

# DW temporal models (Base specification) - only Pre-BTFP and Acute available
dw_temporal_base <- list(
  "Pre-BTFP" = feols(dw_prebtfp ~ mtm_total + uninsured_lev + mtm_x_uninsured + 
                       ln_assets + cash_ratio + loan_to_deposit + book_equity_ratio + wholesale + roa,
                     data = df_dw_pre_s, vcov = "hetero"),
  "Acute" = feols(dw_acute ~ mtm_total + uninsured_lev + mtm_x_uninsured + 
                    ln_assets + cash_ratio + loan_to_deposit + book_equity_ratio + wholesale + roa,
                  data = df_dw_acute, vcov = "hetero")
)

# DW temporal models (Risk specification)
dw_temporal_risk <- list(
  "Pre-BTFP" = feols(dw_prebtfp ~ run_risk_2 + run_risk_3 + run_risk_4 + 
                       ln_assets + cash_ratio + loan_to_deposit + book_equity_ratio + wholesale + roa,
                     data = df_dw_pre_s, vcov = "hetero"),
  "Acute" = feols(dw_acute ~ run_risk_2 + run_risk_3 + run_risk_4 + 
                    ln_assets + cash_ratio + loan_to_deposit + book_equity_ratio + wholesale + roa,
                  data = df_dw_acute, vcov = "hetero")
)

# ------------------------------------------------------------------------------
# PLOT T1: BTFP Temporal Evolution - Base Model (Bar Chart Style)
# ------------------------------------------------------------------------------

btfp_temp_coefs_base <- extract_temporal_coefs(btfp_temporal_base, "BTFP") %>%
  mutate(
    Period = factor(Period, levels = c("Acute", "Post-Acute", "Arbitrage", "Wind-down")),
    Variable = factor(Variable, 
                      levels = c("mtm_x_uninsured", "mtm_total", "uninsured_lev"),
                      labels = c("MTM × Uninsured", "MTM Loss (Total)", "Uninsured Deposits"))
  )

p_btfp_temporal_base <- ggplot(btfp_temp_coefs_base, aes(x = Period, y = Coefficient, fill = Period)) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray30") +
  geom_col(alpha = 0.8, width = 0.7) +
  geom_errorbar(aes(ymin = CI_low, ymax = CI_high), width = 0.2, color = "black") +
  facet_wrap(~ Variable, scales = "free_y", nrow = 1) +
  scale_fill_manual(values = c("Acute" = "#F8766D", "Post-Acute" = "#7CAE00", 
                               "Arbitrage" = "#00BFC4", "Wind-down" = "#C77CFF")) +
  labs(
    title = "BTFP: How Borrowing Determinants Changed Across Crisis Phases",
    subtitle = "Base specification coefficients with 95% confidence intervals",
    x = NULL,
    y = "Coefficient Estimate",
    caption = "Note: LPM estimates. All variables z-standardized."
  ) +
  theme_minimal(base_size = 12) +
  theme(
    plot.title = element_text(face = "bold"),
    legend.position = "none",
    axis.text.x = element_text(angle = 45, hjust = 1),
    strip.text = element_text(face = "bold")
  )

print(p_btfp_temporal_base)

save_figure(p_btfp_temporal_base, "Figure_T1_BTFP_Temporal_Base", width = 14, height = 6)

## Saved: Figure_T1_BTFP_Temporal_Base (PDF + PNG)

# ------------------------------------------------------------------------------
# PLOT T2: BTFP Temporal Evolution - Risk Categories (Bar Chart Style)
# ------------------------------------------------------------------------------

btfp_temp_coefs_risk <- extract_temporal_coefs(btfp_temporal_risk, "BTFP", 
                                                key_vars = c("run_risk_2", "run_risk_3", "run_risk_4")) %>%
  mutate(
    Period = factor(Period, levels = c("Acute", "Post-Acute", "Arbitrage", "Wind-down")),
    Variable = factor(Variable, 
                      levels = c("run_risk_2", "run_risk_3", "run_risk_4"),
                      labels = c("Risk 2\n(Low MTM, High Unins)", 
                                 "Risk 3\n(High MTM, Low Unins)", 
                                 "Risk 4\n(High MTM, High Unins)"))
  )

p_btfp_temporal_risk <- ggplot(btfp_temp_coefs_risk, aes(x = Period, y = Coefficient, fill = Period)) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray30") +
  geom_col(alpha = 0.8, width = 0.7) +
  geom_errorbar(aes(ymin = CI_low, ymax = CI_high), width = 0.2, color = "black") +
  facet_wrap(~ Variable, scales = "free_y", nrow = 1) +
  scale_fill_manual(values = c("Acute" = "#F8766D", "Post-Acute" = "#7CAE00", 
                               "Arbitrage" = "#00BFC4", "Wind-down" = "#C77CFF")) +
  labs(
    title = "BTFP: Risk Category Effects Across Crisis Phases",
    subtitle = "Reference: Risk 1 (Low MTM, Low Uninsured)",
    x = NULL,
    y = "Coefficient Estimate",
    caption = "Note: LPM estimates. Risk categories based on sample median splits."
  ) +
  theme_minimal(base_size = 12) +
  theme(
    plot.title = element_text(face = "bold"),
    legend.position = "none",
    axis.text.x = element_text(angle = 45, hjust = 1),
    strip.text = element_text(face = "bold")
  )

print(p_btfp_temporal_risk)

save_figure(p_btfp_temporal_risk, "Figure_T2_BTFP_Temporal_Risk", width = 14, height = 6)

## Saved: Figure_T2_BTFP_Temporal_Risk (PDF + PNG)

# ------------------------------------------------------------------------------
# PLOT T3: DW Temporal Evolution - Base Model
# ------------------------------------------------------------------------------

dw_temp_coefs_base <- extract_temporal_coefs(dw_temporal_base, "DW") %>%
  mutate(
    Period = factor(Period, levels = c("Pre-BTFP", "Acute")),
    Variable = factor(Variable, 
                      levels = c("mtm_x_uninsured", "mtm_total", "uninsured_lev"),
                      labels = c("MTM × Uninsured", "MTM Loss (Total)", "Uninsured Deposits"))
  )

p_dw_temporal_base <- ggplot(dw_temp_coefs_base, aes(x = Period, y = Coefficient, fill = Period)) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray30") +
  geom_col(alpha = 0.8, width = 0.6) +
  geom_errorbar(aes(ymin = CI_low, ymax = CI_high), width = 0.15, color = "black") +
  facet_wrap(~ Variable, scales = "free_y", nrow = 1) +
  scale_fill_manual(values = c("Pre-BTFP" = "#F8766D", "Acute" = "#00BFC4")) +
  labs(
    title = "Discount Window: How Borrowing Determinants Changed",
    subtitle = "Pre-BTFP (Mar 1-12) vs Acute (Mar 13 - May 1)",
    x = NULL,
    y = "Coefficient Estimate",
    caption = "Note: LPM estimates. All variables z-standardized."
  ) +
  theme_minimal(base_size = 12) +
  theme(
    plot.title = element_text(face = "bold"),
    legend.position = "none",
    strip.text = element_text(face = "bold")
  )

print(p_dw_temporal_base)

save_figure(p_dw_temporal_base, "Figure_T3_DW_Temporal_Base", width = 12, height = 6)

## Saved: Figure_T3_DW_Temporal_Base (PDF + PNG)

# ------------------------------------------------------------------------------
# PLOT T4: DW Temporal Evolution - Risk Categories
# ------------------------------------------------------------------------------

dw_temp_coefs_risk <- extract_temporal_coefs(dw_temporal_risk, "DW", 
                                              key_vars = c("run_risk_2", "run_risk_3", "run_risk_4")) %>%
  mutate(
    Period = factor(Period, levels = c("Pre-BTFP", "Acute")),
    Variable = factor(Variable, 
                      levels = c("run_risk_2", "run_risk_3", "run_risk_4"),
                      labels = c("Risk 2\n(Low MTM, High Unins)", 
                                 "Risk 3\n(High MTM, Low Unins)", 
                                 "Risk 4\n(High MTM, High Unins)"))
  )

p_dw_temporal_risk <- ggplot(dw_temp_coefs_risk, aes(x = Period, y = Coefficient, fill = Period)) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray30") +
  geom_col(alpha = 0.8, width = 0.6) +
  geom_errorbar(aes(ymin = CI_low, ymax = CI_high), width = 0.15, color = "black") +
  facet_wrap(~ Variable, scales = "free_y", nrow = 1) +
  scale_fill_manual(values = c("Pre-BTFP" = "#F8766D", "Acute" = "#00BFC4")) +
  labs(
    title = "Discount Window: Risk Category Effects Over Time",
    subtitle = "Reference: Risk 1 (Low MTM, Low Uninsured)",
    x = NULL,
    y = "Coefficient Estimate",
    caption = "Note: LPM estimates. Risk categories based on sample median splits."
  ) +
  theme_minimal(base_size = 12) +
  theme(
    plot.title = element_text(face = "bold"),
    legend.position = "none",
    strip.text = element_text(face = "bold")
  )

print(p_dw_temporal_risk)

save_figure(p_dw_temporal_risk, "Figure_T4_DW_Temporal_Risk", width = 12, height = 6)

## Saved: Figure_T4_DW_Temporal_Risk (PDF + PNG)

# ------------------------------------------------------------------------------
# PLOT T5: Combined Temporal - BTFP Base (Line Plot Alternative)
# ------------------------------------------------------------------------------

p_btfp_temporal_line <- ggplot(btfp_temp_coefs_base, aes(x = Period, y = Coefficient, 
                                                          group = Variable, color = Variable)) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray50") +
  geom_line(linewidth = 1) +
  geom_pointrange(aes(ymin = CI_low, ymax = CI_high), size = 0.8) +
  scale_color_manual(values = c("MTM × Uninsured" = "#E41A1C", 
                                "MTM Loss (Total)" = "#377EB8", 
                                "Uninsured Deposits" = "#4DAF4A")) +
  labs(
    title = "BTFP: Evolution of Key Determinants Across Crisis Phases",
    subtitle = "How the role of fundamentals and liquidity mismatch changed over time",
    x = NULL,
    y = "Coefficient Estimate (95% CI)",
    color = "Variable",
    caption = "Note: LPM estimates. All variables z-standardized."
  ) +
  theme_minimal(base_size = 12) +
  theme(
    plot.title = element_text(face = "bold"),
    legend.position = "bottom",
    axis.text.x = element_text(angle = 45, hjust = 1)
  )

print(p_btfp_temporal_line)

save_figure(p_btfp_temporal_line, "Figure_T5_BTFP_Temporal_Line", width = 12, height = 7)

## Saved: Figure_T5_BTFP_Temporal_Line (PDF + PNG)

# ------------------------------------------------------------------------------
# PLOT T6: Combined Temporal - Risk 4 Comparison (BTFP vs DW)
# ------------------------------------------------------------------------------

# Combine Risk 4 coefficients from both facilities
risk4_temporal <- bind_rows(
  btfp_temp_coefs_risk %>% filter(Variable == "Risk 4\n(High MTM, High Unins)") %>% 
    mutate(Facility = "BTFP"),
  dw_temp_coefs_risk %>% filter(Variable == "Risk 4\n(High MTM, High Unins)") %>% 
    mutate(Facility = "DW")
) %>%
  mutate(Period = as.character(Period))

p_risk4_comparison <- ggplot(risk4_temporal, aes(x = Period, y = Coefficient, fill = Facility)) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray30") +
  geom_col(position = position_dodge(width = 0.7), alpha = 0.8, width = 0.6) +
  geom_errorbar(aes(ymin = CI_low, ymax = CI_high), 
                position = position_dodge(width = 0.7), width = 0.2, color = "black") +
  scale_fill_manual(values = c("BTFP" = "#2166AC", "DW" = "#B2182B")) +
  labs(
    title = "Risk 4 (High MTM, High Uninsured) Effect: BTFP vs DW",
    subtitle = "Comparing facility choice among highest-risk banks across crisis phases",
    x = NULL,
    y = "Coefficient Estimate (95% CI)",
    fill = "Facility",
    caption = "Note: LPM estimates. Reference: Risk 1 (Low MTM, Low Uninsured)."
  ) +
  theme_minimal(base_size = 12) +
  theme(
    plot.title = element_text(face = "bold"),
    legend.position = "bottom"
  )

print(p_risk4_comparison)

save_figure(p_risk4_comparison, "Figure_T6_Risk4_BTFP_vs_DW", width = 10, height = 6)

## Saved: Figure_T6_Risk4_BTFP_vs_DW (PDF + PNG)

10 APPENDIX: VARIABLE DEFINITIONS

var_defs <- tribble(
  ~Variable, ~Definition,
  "btfp_[period]", "1 if bank borrowed from BTFP during specified period, 0 if pure non-user",
  "dw_[period]", "1 if bank borrowed from DW during specified period, 0 if pure non-user",
  "any_fed", "1 if bank used BTFP or DW during period, 0 if pure non-user",
  "all_user", "1 if bank used BTFP, DW, or had abnormal FHLB increase",
  "non_user", "1 if bank did NOT access any facility (BTFP=0, DW=0, FHLB=0)",
  "user_group", "Categorical: BTFP_Only, DW_Only, Both, Neither",
  "mtm_total", "Z-score: Total MTM loss / Total Assets",
  "uninsured_lev", "Z-score: Uninsured deposits / Total Assets",
  "mtm_x_uninsured", "Interaction: mtm_total × uninsured_lev",
  "run_risk_1", "1 if MTM < median AND uninsured < median (Reference)",
  "run_risk_2", "1 if MTM < median AND uninsured >= median",
  "run_risk_3", "1 if MTM >= median AND uninsured < median",
  "run_risk_4", "1 if MTM >= median AND uninsured >= median",
  "size_cat", "Small (<$1B), Medium ($1B-$100B), Large (>$100B)",
  "ln_assets", "Z-score: Log(Total Assets)",
  "cash_ratio", "Z-score: Cash / Total Assets",
  "loan_to_deposit", "Z-score: Total Loans / Total Deposits",
  "book_equity_ratio", "Z-score: Book Equity / Total Assets",
  "wholesale", "Z-score: Wholesale Funding / Total Liabilities",
  "roa", "Z-score: Return on Assets"
)

var_defs %>%
  kable(format = "html", caption = "Variable Definitions",
        col.names = c("Variable", "Definition")) %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = TRUE, font_size = 10) %>%
  pack_rows("Dependent Variables", 1, 6) %>%
  pack_rows("Key Explanatory (Z-Standardized)", 7, 9) %>%
  pack_rows("Run Risk Categories", 10, 13) %>%
  pack_rows("Size Categories", 14, 14) %>%
  pack_rows("Controls (Z-Standardized)", 15, 20)

Variable Definitions

Variable	Definition
Dependent Variables
btfp_[period]	1 if bank borrowed from BTFP during specified period, 0 if pure non-user
dw_[period]	1 if bank borrowed from DW during specified period, 0 if pure non-user
any_fed	1 if bank used BTFP or DW during period, 0 if pure non-user
all_user	1 if bank used BTFP, DW, or had abnormal FHLB increase
non_user	1 if bank did NOT access any facility (BTFP=0, DW=0, FHLB=0)
user_group	Categorical: BTFP_Only, DW_Only, Both, Neither
Key Explanatory (Z-Standardized)
mtm_total	Z-score: Total MTM loss / Total Assets
uninsured_lev	Z-score: Uninsured deposits / Total Assets
mtm_x_uninsured	Interaction: mtm_total × uninsured_lev
Run Risk Categories
run_risk_1	1 if MTM < median AND uninsured < median (Reference)
run_risk_2	1 if MTM < median AND uninsured >= median
run_risk_3	1 if MTM >= median AND uninsured < median
run_risk_4	1 if MTM >= median AND uninsured >= median
Size Categories
size_cat	Small (<$1B), Medium ($1B-$100B), Large (>$100B)
Controls (Z-Standardized)
ln_assets	Z-score: Log(Total Assets)
cash_ratio	Z-score: Cash / Total Assets
loan_to_deposit	Z-score: Total Loans / Total Deposits
book_equity_ratio	Z-score: Book Equity / Total Assets
wholesale	Z-score: Wholesale Funding / Total Liabilities
roa	Z-score: Return on Assets

# Save Variable Definitions
save_kable_table(var_defs, "Appendix_Variable_Definitions",
                 "Variable Definitions", 
                 "MTM Loss captures bank fundamentals (unrealized securities losses). Uninsured Deposits capture liquidity mismatch (run vulnerability). The interaction tests whether liquidity risk amplifies fundamental weakness. Risk categories based on sample medians. All continuous variables winsorized at 2.5/97.5 percentiles before z-standardization.")

Saved: Appendix_Variable_Definitions (HTML + LaTeX)

11 ============================================================================

12 SESSION INFO

13 ============================================================================

cat("\n============================================\n")

## 
## ============================================

cat("ANALYSIS COMPLETE\n")

## ANALYSIS COMPLETE

cat("============================================\n\n")

## ============================================

cat("Tables saved to:", TABLE_PATH, "\n")

## Tables saved to: C:/Users/mferdo2/OneDrive - Louisiana State University/Finance_PhD/DW_Stigma_paper/Liquidity_project_2025/03_documentation/crisis_borrowing_result/tables

cat("Figures saved to:", FIG_PATH, "\n\n")

## Figures saved to: C:/Users/mferdo2/OneDrive - Louisiana State University/Finance_PhD/DW_Stigma_paper/Liquidity_project_2025/03_documentation/crisis_borrowing_result/figures

cat("Summary of outputs:\n")

## Summary of outputs:

cat("- Tables 1A-1C: Summary Statistics\n")

## - Tables 1A-1C: Summary Statistics

cat("- Tables 2A-2D: User Group Comparisons (with p-values)\n")

## - Tables 2A-2D: User Group Comparisons (with p-values)

cat("- Tables 3-6: Main Regression Tables\n")

## - Tables 3-6: Main Regression Tables

cat("- Tables 4A-4E: BTFP Temporal Analysis\n")

## - Tables 4A-4E: BTFP Temporal Analysis

cat("- Tables 4C-4D: DW Temporal Analysis\n")

## - Tables 4C-4D: DW Temporal Analysis

cat("- Tables 5A-5J: Size Analysis\n")

## - Tables 5A-5J: Size Analysis

cat("- Figures 1-10: Publication-Quality Plots\n")

## - Figures 1-10: Publication-Quality Plots

cat("- Appendix: Variable Definitions\n\n")

## - Appendix: Variable Definitions

sessionInfo()

## R version 4.3.1 (2023-06-16 ucrt)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 11 x64 (build 26100)
## 
## Matrix products: default
## 
## 
## locale:
## [1] LC_COLLATE=English_United States.utf8 
## [2] LC_CTYPE=English_United States.utf8   
## [3] LC_MONETARY=English_United States.utf8
## [4] LC_NUMERIC=C                          
## [5] LC_TIME=English_United States.utf8    
## 
## time zone: America/Chicago
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] patchwork_1.3.0    scales_1.3.0       kableExtra_1.4.0   modelsummary_2.3.0
##  [5] fixest_0.12.1      data.table_1.17.0  lubridate_1.9.4    forcats_1.0.0     
##  [9] stringr_1.5.1      dplyr_1.1.4        purrr_1.0.2        readr_2.1.5       
## [13] tidyr_1.3.1        tibble_3.2.1       ggplot2_3.5.1      tidyverse_2.0.0   
## 
## loaded via a namespace (and not attached):
##  [1] gtable_0.3.6        bayestestR_0.17.0   xfun_0.54          
##  [4] bslib_0.5.1         insight_1.4.3       lattice_0.21-8     
##  [7] tzdb_0.4.0          numDeriv_2016.8-1.1 vctrs_0.6.4        
## [10] tools_4.3.1         generics_0.1.3      datawizard_1.3.0   
## [13] parallel_4.3.1      sandwich_3.1-1      fansi_1.0.5        
## [16] pkgconfig_2.0.3     tinytable_0.15.1    checkmate_2.3.1    
## [19] stringmagic_1.1.2   lifecycle_1.0.4     compiler_4.3.1     
## [22] farver_2.1.1        textshaping_0.3.7   munsell_0.5.0      
## [25] htmltools_0.5.9     sass_0.4.9          yaml_2.3.7         
## [28] Formula_1.2-5       pillar_1.9.0        crayon_1.5.3       
## [31] jquerylib_0.1.4     cachem_1.0.8        nlme_3.1-162       
## [34] tidyselect_1.2.1    digest_0.6.33       performance_0.15.3 
## [37] stringi_1.7.12      labeling_0.4.3      fastmap_1.1.1      
## [40] grid_4.3.1          colorspace_2.1-0    cli_3.6.5          
## [43] magrittr_2.0.3      utf8_1.2.4          withr_2.5.2        
## [46] dreamerr_1.4.0      backports_1.4.1     bit64_4.0.5        
## [49] timechange_0.3.0    rmarkdown_2.29      bit_4.0.5          
## [52] ragg_1.3.0          zoo_1.8-12          hms_1.1.3          
## [55] evaluate_0.23       knitr_1.50          parameters_0.28.3  
## [58] viridisLite_0.4.2   rlang_1.1.6         Rcpp_1.0.12        
## [61] glue_1.8.0          xml2_1.3.6          svglite_2.1.3      
## [64] rstudioapi_0.16.0   vroom_1.6.5         jsonlite_1.8.7     
## [67] R6_2.5.1            tables_0.9.31       systemfonts_1.0.6