1 Theory-to-Code Map

Every equation in theory_standalone_final.pdf is mapped below.

Theory Ref Expression Code Variable Section
Eq. 1 \(C+S+L = D^I + D^U + E\) Balance sheet identity §3
Eq. 2 \(\lambda = (1-\theta_S)S + (1-\theta_L)L\) mtm_total_raw §3
Eq. 3 \(F = (1-\beta^U)r\,D^U/\delta\) f_pp (via DSSW Eq.20) §3
Eq. 4 \(E^{MV} = E - \lambda + F\) emv_pp §3
Eq. 5 \(F^U = (D^U/D)F\) f_u_pp §3
Eq. 6 \(g = w - (C + \theta_S S)\) liquid_buffer_raw §3
Eq. 7 \(v = E^{MV} - F^U\) v_pp §3
Eq. 11 \(\partial^2\Pr(\text{borrow})/\partial\lambda\,\partial F^U > 0\) Interaction test §7
Eq. 12 \(\lambda^\dagger = E + F - F^U\) threshold_run §3
Eq. 13 LPM: \(\Pr(\text{Borrow}_i) = \alpha + c_1\ell + c_2 f^U + c_3(\ell\times f^U) + \gamma'X + \varepsilon\) Main regression §7
Eq. 14 \(\Delta_{\text{BTFP}} = (1-\theta_S)S_{\text{el}}\) mtm_btfp_raw §8
Eq. 15 \(\text{Amount}_i = \alpha + d_1\ell_{S,i} + \gamma'X + \varepsilon\) Intensive margin §8
Falsif. FHLB: \(c_3 \approx 0\); Arb period: \(c_3 \approx 0\) §9

2 SETUP

2.1 Packages

rm(list = ls())
library(data.table); library(dplyr); library(tidyr); library(stringr)
library(lubridate); library(purrr); library(tibble)
library(fixest); library(modelsummary)
library(knitr); library(kableExtra)
library(ggplot2); library(scales); library(patchwork)
library(readr); library(readxl)
cat("All packages loaded.\n")
## All packages loaded.

2.2 Helpers

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])
}

standardize_z <- function(x) {
  if (all(is.na(x))) return(x)
  s <- sd(x, na.rm = TRUE)
  if (is.na(s) || s == 0) return(rep(0, length(x)))
  (x - mean(x, na.rm = TRUE)) / s
}

safe_div <- function(num, denom, default = NA_real_) {
  ifelse(is.na(denom) | denom == 0, default, num / denom)
}

format_pval <- function(p) {
  case_when(is.na(p) ~ "", p < 0.01 ~ "***", p < 0.05 ~ "**", p < 0.10 ~ "*", TRUE ~ "")
}

save_figure <- function(plot_obj, filename, width = 12, height = 8) {
  ggsave(file.path(FIG_PATH, paste0(filename, ".pdf")),
    plot = plot_obj, width = width, height = height, device = "pdf")
  message("Saved: ", filename, ".pdf")
}

dv_summary_rows <- function(model_list) {
  stats <- lapply(model_list, function(m) {
    y <- fitted(m) + resid(m)
    c(N_DV1 = sum(round(y) == 1, na.rm = TRUE),
      Mean_DV = round(mean(y, na.rm = TRUE), 4))
  })
  df <- data.frame(term = c("N (DV = 1)", "Mean (DV)"), stringsAsFactors = FALSE)
  for (nm in names(stats)) df[[nm]] <- as.character(stats[[nm]])
  df
}

save_reg_latex <- function(model_list, filename, ...) {
  tex_path <- file.path(TABLE_PATH, paste0(filename, ".tex"))
  msummary(model_list, output = tex_path,
           stars = c("*" = .10, "**" = .05, "***" = .01),
           gof_omit = "AIC|BIC|Log|RMSE", ...)
  cat(sprintf("Saved: %s.tex\n", filename))
}

save_kbl_latex <- function(df, filename, col.names = NULL, caption = "", align = NULL) {
  tex <- kbl(df, format = "latex", booktabs = TRUE, escape = FALSE,
             col.names = col.names, caption = caption, align = align) %>%
    kable_styling(latex_options = c("hold_position", "scale_down"))
  writeLines(tex, file.path(TABLE_PATH, paste0(filename, ".tex")))
  cat(sprintf("Saved: %s.tex\n", filename))
}

2.3 Paths & 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/new_run_Analysis_revised")
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_MAIN <- "2022Q4"; BASELINE_ARB <- "2023Q3"
CRISIS_START <- as.Date("2023-03-08"); CRISIS_END <- as.Date("2023-05-04")
ARB_START <- as.Date("2023-11-15"); ARB_END <- as.Date("2024-01-24")
DW_DATA_END <- as.Date("2023-12-31")

3 DATA LOADING

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))
dssw_betas <- read_csv(file.path(DATA_PROC, "dssw_deposit_betas.csv"), show_col_types = FALSE) %>%
  mutate(idrssd = as.character(idrssd))
dssw_beta_2022q4 <- dssw_betas %>% filter(estimation_date == "2022Q4") %>%
  select(idrssd, beta_overall, beta_insured, beta_uninsured,
         beta_insured_w, beta_uninsured_w, gamma_hat, alpha_hat)
public_flag <- read_csv(file.path(DATA_PROC, "public_bank_flag.csv"), show_col_types = FALSE) %>%
  mutate(idrssd = as.character(idrssd)) %>% select(idrssd, period, is_public)

deposit_costs_file <- file.path(DATA_PROC, "dssw_deposit_costs.csv")
if (file.exists(deposit_costs_file)) {
  deposit_costs <- read_csv(deposit_costs_file, show_col_types = FALSE) %>%
    mutate(idrssd = as.character(idrssd))
  deposit_costs_2022q4 <- deposit_costs %>%
    filter(period == "2022Q4") %>%
    select(idrssd, deposit_cost_weighted, deposit_cost_insured, deposit_cost_uninsured)
  HAS_DEPOSIT_COSTS <- TRUE
} else {
  deposit_costs_2022q4 <- NULL
  HAS_DEPOSIT_COSTS <- FALSE
}
cat("Call Report:", nrow(call_q), "obs |", n_distinct(call_q$idrssd), "banks\n")
## Call Report: 75989 obs | 5074 banks

3.1 Exclusions & Borrower Indicators

excluded_banks <- call_q %>%
  filter(period == BASELINE_MAIN, failed_bank == 1 | gsib == 1) %>% pull(idrssd)
btfp_loans <- btfp_loans_raw %>% filter(!rssd_id %in% excluded_banks)
dw_loans   <- dw_loans_raw   %>% filter(!rssd_id %in% excluded_banks)

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_crisis <- create_borrower_indicator(btfp_loans,"btfp_loan_date","rssd_id","btfp_loan_amount",CRISIS_START,CRISIS_END,"btfp_crisis")
btfp_arb    <- create_borrower_indicator(btfp_loans,"btfp_loan_date","rssd_id","btfp_loan_amount",ARB_START,ARB_END,"btfp_arb")
dw_crisis   <- create_borrower_indicator(dw_loans,"dw_loan_date","rssd_id","dw_loan_amount",CRISIS_START,min(CRISIS_END,DW_DATA_END),"dw_crisis")
dw_arb      <- create_borrower_indicator(dw_loans,"dw_loan_date","rssd_id","dw_loan_amount",ARB_START,DW_DATA_END,"dw_arb")
cat("BTFP Crisis:", nrow(btfp_crisis), "| DW Crisis:", nrow(dw_crisis), "\n")
## BTFP Crisis: 526 | DW Crisis: 459

4 CALIBRATION & VARIABLE CONSTRUCTION

4.1 Calibration [DSSW Eq. 20]

# ==============================================================================
# DSSW Eq. (20): Decaying perpetuity
#   f_i = pmax( ((1 - β^U_i) * y - c^U_i) / (y + δ) × μ_i,  0 )
#
# Theory Eq. (3) writes F = (1-β)rD^U/δ for exposition.
# Empirically we use the DSSW decaying-perpetuity capitalization:
#   y     = 10-year Treasury yield (discount rate matching 1/δ duration)
#   δ     = 0.10 (deposit decay rate; DSSW baseline)
#   c^U_i = bank-level uninsured deposit operating cost (DSSW hedonic)
# ==============================================================================

y_10yr      <- 0.0392    # 10-year Treasury yield (Feb 2023) — DSSW calibration
delta_decay <- 0.10      # Deposit decay rate (DSSW baseline)
cap_factor  <- 1 / (y_10yr + delta_decay)  # = 7.18

cat("=== CALIBRATION (DSSW Eq. 20) ===\n")
## === CALIBRATION (DSSW Eq. 20) ===
cat(sprintf("  y (10yr Treasury) = %.4f\n", y_10yr))
##   y (10yr Treasury) = 0.0392
cat(sprintf("  δ (decay rate)    = %.2f\n", delta_decay))
##   δ (decay rate)    = 0.10
cat(sprintf("  Cap factor        = 1/(y+δ) = 1/%.4f = %.2f\n", y_10yr + delta_decay, cap_factor))
##   Cap factor        = 1/(y+δ) = 1/0.1392 = 7.18

4.2 Variable Construction [Theory Eq. 1–7, 12, Table 1]

# ==============================================================================
# Maps EVERY variable in theory Table 1 to data:
#   ℓ   = total MTM losses / A                    (Eq. 2)
#   ℓ_S = BTFP eligible securities MTM / A                      (Eq. 2, Eq. 14)
#   ℓ_L = loan MTM / A                            (Eq. 2)
#   f   = pmax(((1-β^U)y - c^U)/(y+δ) × μ, 0)    (Eq. 3, DSSW Eq. 20)
#   e^MV = e - ℓ + f                              (Eq. 4)
#   f^U = (D^U/D) × f                             (Eq. 5) ← KEY NEW VARIABLE
#   v   = e^MV - f^U = e - ℓ + (D^I/D)f           (Eq. 7)
#   λ†  = E + F - F^U = e + (D^I/D)f              (Eq. 12)
# ==============================================================================

construct_analysis_vars <- function(baseline_data) {
  baseline_data %>%
    mutate(
      # ── Eq. 2: MTM loss components ──
      mtm_total_raw  = mtm_loss_to_total_asset,
      mtm_loan_raw   = mtm_loss_total_loan_to_total_asset,           # ℓ_L: Loan losses
      mtm_sec_raw    = mtm_loss_to_total_asset -                     # ℓ_S: Securities losses
                       mtm_loss_total_loan_to_total_asset,           #      = ℓ - ℓ_L
      # BTFP-eligible (OMO) vs non-BTFP decomposition (for P3 par-channel tests)
      mtm_btfp_raw   = mtm_loss_omo_eligible_to_total_asset,        # OMO-eligible securities
      mtm_other_raw  = mtm_loss_non_omo_eligible_to_total_asset,    # non-OMO sec + loans
      # ── Balance sheet ratios ──
      uninsured_lev_raw     = uninsured_deposit_to_total_asset,
      insured_lev_raw       = r_insured_deposit,
      uninsured_share_raw   = uninsured_to_deposit,
      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,
      roa_raw               = roa,
      fhlb_ratio_raw        = fhlb_to_total_asset,
      loan_to_deposit_raw   = loan_to_deposit,
      wholesale_raw = safe_div(
        replace_na(fed_fund_purchase, 0) + replace_na(repo, 0) +
          replace_na(other_borrowed_less_than_1yr, 0),
        total_liability, 0) * 100,

      # ── Deposit beta: clip to [0,1] per theory ──
      uninsured_beta_raw = ifelse(!is.na(beta_uninsured), beta_uninsured, NA_real_),
      beta_u_clipped     = pmin(pmax(uninsured_beta_raw, 0), 1),

      # ── Uninsured deposit cost (DSSW hedonic) ──
      cost_u_raw = ifelse(!is.na(deposit_cost_uninsured), deposit_cost_uninsured, 0),

      # ── Jiang AE: adjusted equity = e - ℓ ──
      adjusted_equity_raw = book_equity_to_total_asset - mtm_loss_to_total_asset,

      # ── Key ratios ──
      mu_decimal        = uninsured_deposit / total_asset,
      insured_share     = safe_div(insured_deposit, insured_deposit + uninsured_deposit, NA_real_),
      uninsured_share_d = safe_div(uninsured_deposit, insured_deposit + uninsured_deposit, NA_real_),

      # ── Eq. 3 (DSSW Eq. 20): Total franchise value f ──
      gross_rent = (1 - beta_u_clipped) * y_10yr,
      net_rent   = gross_rent - cost_u_raw,
      f_decimal  = ifelse(!is.na(beta_u_clipped),
        pmax(net_rent * cap_factor * mu_decimal, 0), NA_real_),
      f_pp = ifelse(!is.na(f_decimal), f_decimal * 100, NA_real_),

      # ── Eq. 4: E^MV = e - ℓ + f ──
      emv_pp = ifelse(!is.na(f_pp),
        book_equity_to_total_asset - mtm_total_raw + f_pp, NA_real_),

      # ── Eq. 5: F^U = (D^U/D) × f  [THE KEY VARIABLE] ──
      f_u_pp = ifelse(!is.na(f_pp) & !is.na(uninsured_share_d),
        uninsured_share_d * f_pp, NA_real_),

      # ── Eq. 7: v = E^MV - F^U  [run value: v < 0 ⟹ run possible] ──
      v_pp = ifelse(!is.na(emv_pp) & !is.na(f_u_pp),
        emv_pp - f_u_pp, NA_real_),

      # ── Eq. 12: λ† = e + f - f^U = e + (D^I/D)f  [run threshold] ──
      threshold_run = ifelse(!is.na(f_pp) & !is.na(insured_share),
        book_equity_to_total_asset + insured_share * f_pp, NA_real_),

      # ── Eq. 6: Liquid buffer = C/A + θ_S·S/A ──
      liquid_buffer_raw = cash_to_total_asset + security_to_total_asset - mtm_btfp_raw ,

      # ── Run indicator: v < 0 ──
      run_possible = as.integer(!is.na(v_pp) & v_pp < 0),

      # ── Collateral (for intensive margin, Eq. 14) ──
      collateral_capacity_raw = omo_eligible_to_total_asset,

      # ── Winsorized ──
     
      mtm_total_w=winsorize(mtm_total_raw), mtm_sec_w=winsorize(mtm_sec_raw),
      mtm_loan_w=winsorize(mtm_loan_raw),
      mtm_btfp_w=winsorize(mtm_btfp_raw),   # BTFP-eligible (OMO)
      mtm_other_w=winsorize(mtm_other_raw),  # non-BTFP (non-OMO sec + loans)
      uninsured_lev_w=winsorize(uninsured_lev_raw),
      insured_lev_w=winsorize(r_insured_deposit), adjusted_equity_w=winsorize(adjusted_equity_raw),
      f_pp_w=winsorize(f_pp), f_u_pp_w=winsorize(f_u_pp), emv_pp_w=winsorize(emv_pp),
      uninsured_beta_w=winsorize(uninsured_beta_raw), ln_assets_w=winsorize(ln_assets_raw),
      cash_ratio_w=winsorize(cash_ratio_raw), book_equity_ratio_w=winsorize(book_equity_ratio_raw),
      roa_w=winsorize(roa_raw), loan_to_deposit_w=winsorize(loan_to_deposit_raw),
      wholesale_w=winsorize(wholesale_raw),
      deposit_cost_w=winsorize(cost_u_raw),

      # ── Z-standardized (full-sample z-scores) ──
      mtm_total=standardize_z(mtm_total_w), mtm_sec=standardize_z(mtm_sec_w),
      mtm_loan=standardize_z(mtm_loan_w),
      mtm_btfp=standardize_z(mtm_btfp_w),   # BTFP-eligible (OMO) z-score
      mtm_other=standardize_z(mtm_other_w),  # non-BTFP z-score
      uninsured_lev=standardize_z(uninsured_lev_w),
      insured_lev=standardize_z(insured_lev_w), adjusted_equity=standardize_z(adjusted_equity_w),
      franchise_value=standardize_z(f_pp_w), uninsured_franchise=standardize_z(f_u_pp_w), emv=standardize_z(emv_pp_w),
      uninsured_beta=standardize_z(uninsured_beta_w), ln_assets=standardize_z(ln_assets_w),
      cash_ratio=standardize_z(cash_ratio_w), book_equity_ratio=standardize_z(book_equity_ratio_w),
      roa=standardize_z(roa_w), loan_to_deposit=standardize_z(loan_to_deposit_w),
      wholesale=standardize_z(wholesale_w),

      # ── Eq. 11/13: THE interaction ──
      # ℓ × f^U  (NOT ℓ × f — the theory's cross-partial is w.r.t. F^U)
      mtm_x_f_u = mtm_total * uninsured_franchise
    )
}

5 BUILD DATASETS

df_2022q4 <- call_q %>%
  filter(period == BASELINE_MAIN, !idrssd %in% excluded_banks,
         !is.na(omo_eligible) & omo_eligible > 0) %>%
  left_join(dssw_beta_2022q4, by = "idrssd") %>%
  { if (HAS_DEPOSIT_COSTS) left_join(., deposit_costs_2022q4, by = "idrssd") else . } %>%
  left_join(public_flag %>% filter(period == "2022Q4") %>% select(idrssd, is_public), by = "idrssd") %>%
  mutate(is_public = replace_na(is_public, 0L)) %>%
  construct_analysis_vars()

df_2023q3 <- call_q %>%
  filter(period == BASELINE_ARB, !idrssd %in% excluded_banks,
         !is.na(omo_eligible) & omo_eligible > 0) %>%
  left_join(dssw_beta_2022q4, by = "idrssd") %>%
  { if (HAS_DEPOSIT_COSTS) left_join(., deposit_costs_2022q4, by = "idrssd") else . } %>%
  left_join(public_flag %>% filter(period == "2023Q3") %>% select(idrssd, is_public), by = "idrssd") %>%
  mutate(is_public = replace_na(is_public, 0L)) %>%
  construct_analysis_vars()

cat("2022Q4:", nrow(df_2022q4), "| β^U available:", sum(!is.na(df_2022q4$f_pp)), "\n")
## 2022Q4: 4292 | β^U available: 4226
join_all_borrowers <- function(df_base, btfp_df, dw_df, btfp_var, dw_var) {
  df_base %>%
    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("{btfp_var}" := replace_na(!!sym(btfp_var), 0L),
           "{dw_var}" := replace_na(!!sym(dw_var), 0L),
           fhlb_user = as.integer(abnormal_fhlb_borrowing_10pct == 1),
           any_fed = as.integer(!!sym(btfp_var) == 1 | !!sym(dw_var) == 1),
           both_fed = as.integer(!!sym(btfp_var) == 1 & !!sym(dw_var) == 1),
           # Total Fed amount borrowed (for intensive margin)
           fed_amt = replace_na(!!sym(paste0(btfp_var, "_amt")), 0) +
                     replace_na(!!sym(paste0(dw_var, "_amt")), 0))
}

df_crisis <- join_all_borrowers(df_2022q4, btfp_crisis, dw_crisis, "btfp_crisis", "dw_crisis")
df_arb <- df_2023q3 %>%
  left_join(btfp_arb %>% select(idrssd, btfp_arb, btfp_arb_amt), by = "idrssd") %>%
  left_join(dw_arb %>% select(idrssd, dw_arb), by = "idrssd") %>%
  mutate(btfp_arb = replace_na(btfp_arb, 0L), dw_arb = replace_na(dw_arb, 0L),
         fhlb_user = as.integer(abnormal_fhlb_borrowing_10pct == 1),
         any_fed = as.integer(btfp_arb == 1 | dw_arb == 1))

cat("Crisis:", nrow(df_crisis), "| Any Fed:", sum(df_crisis$any_fed), "\n")
## Crisis: 4292 | Any Fed: 828

5.1 Sample Cleaning

core_vars <- c("book_equity_to_total_asset", "total_asset", "total_liability",
               "cash_to_total_asset", "security_to_total_asset",
               "total_loan_to_total_asset", "loan_to_deposit", "roa",
               "uninsured_deposit_to_total_asset", "wholesale_raw")

df_crisis_clean <- df_crisis %>%
  filter(!is.na(mtm_loss_to_total_asset)) %>%
  filter(if_all(all_of(core_vars), ~ !is.na(.))) %>%
  filter(!is.na(uninsured_deposit) & !is.na(insured_deposit) & insured_deposit > 0) %>%
  filter(!is.na(abnormal_fhlb_borrowing_10pct)) %>%
  filter(total_asset > 0 & is.finite(ln_assets_raw))

# Theory model sample: require β^U (and hence f, f^U, v)
df_model <- df_crisis_clean %>% filter(!is.na(f_u_pp))

df_arb_clean <- df_arb %>%
  filter(!is.na(mtm_loss_to_total_asset)) %>%
  filter(if_all(all_of(core_vars), ~ !is.na(.))) %>%
  filter(!is.na(uninsured_deposit) & !is.na(insured_deposit) & insured_deposit > 0) %>%
  filter(!is.na(abnormal_fhlb_borrowing_10pct)) %>%
  filter(total_asset > 0)

cat(sprintf("Clean: %d | Theory model (β^U available): %d | Arb: %d\n",
    nrow(df_crisis_clean), nrow(df_model), nrow(df_arb_clean)))
## Clean: 4251 | Theory model (β^U available): 4226 | Arb: 4168
cat(sprintf("Fed borrowers in model sample: %d (%.1f%%)\n",
    sum(df_model$any_fed), 100 * mean(df_model$any_fed)))
## Fed borrowers in model sample: 822 (19.5%)

6 DESCRIPTIVE STATISTICS

6.1 Franchise Diagnostics

cat("=== FRANCHISE & RUN VALUE DIAGNOSTICS ===\n\n")
## === FRANCHISE & RUN VALUE DIAGNOSTICS ===
cat(sprintf("Banks: %d | Borrowers: %d (%.1f%%)\n",
    nrow(df_model), sum(df_model$any_fed), 100 * mean(df_model$any_fed)))
## Banks: 4226 | Borrowers: 822 (19.5%)
cat(sprintf("\nTotal franchise f: mean=%.2f, med=%.2f, sd=%.2f\n",
    mean(df_model$f_pp), median(df_model$f_pp), sd(df_model$f_pp)))
## 
## Total franchise f: mean=2.66, med=2.38, sd=1.61
cat(sprintf("Uninsured franchise f^U: mean=%.2f, med=%.2f, sd=%.2f\n",
    mean(df_model$f_u_pp, na.rm=T), median(df_model$f_u_pp, na.rm=T), sd(df_model$f_u_pp, na.rm=T)))
## Uninsured franchise f^U: mean=0.91, med=0.62, sd=1.05
cat(sprintf("Run value v: mean=%.2f, med=%.2f, sd=%.2f\n",
    mean(df_model$v_pp, na.rm=T), median(df_model$v_pp, na.rm=T), sd(df_model$v_pp, na.rm=T)))
## Run value v: mean=5.86, med=5.33, sd=6.31
cat(sprintf("\nInputs: β^U(clipped): mean=%.3f | μ: mean=%.3f | D^U/D: mean=%.3f\n",
    mean(df_model$beta_u_clipped, na.rm=T), mean(df_model$mu_decimal, na.rm=T),
    mean(df_model$uninsured_share_d, na.rm=T)))
## 
## Inputs: β^U(clipped): mean=0.331 | μ: mean=0.238 | D^U/D: mean=0.277
cat(sprintf("  c^U: mean=%.4f | Net rent ((1-β^U)y - c^U): mean=%.4f\n",
    mean(df_model$cost_u_raw, na.rm=T),
    mean(df_model$net_rent, na.rm=T)))
##   c^U: mean=0.0108 | Net rent ((1-β^U)y - c^U): mean=0.0154
cat(sprintf("  Cap factor = 1/(y+δ) = %.2f | y=%.2f%%, δ=%.0f%%\n",
    cap_factor, y_10yr*100, delta_decay*100))
##   Cap factor = 1/(y+δ) = 7.18 | y=3.92%, δ=10%

6.2 Borrower vs Non-Borrower

desc_vars <- c("ln_assets_raw", "mtm_total_raw", "mtm_btfp_raw", "mtm_other_raw",
  "book_equity_ratio_raw", "cash_ratio_raw", "loan_to_deposit_raw",
  "wholesale_raw", "roa_raw",
  "uninsured_lev_raw", "uninsured_share_raw", "uninsured_beta_raw",
  "cost_u_raw", "f_pp", "f_u_pp", "emv_pp", "v_pp")
desc_labels <- c("Log(Assets)", "Total MTM (ℓ)", "Securities MTM (ℓ_S)", "Loan MTM (ℓ_L)",
  "Book Equity (e)", "Cash/TA", "Loan/Deposit",
  "Wholesale (%)", "ROA",
  "Uninsured/TA (μ)", "D^U/D", "β^U",
  "Cost c^U", "Franchise (f)", "Uninsured Franchise (f^U)", "MV Equity (E^MV)", "Run Value (v)")

desc_tbl <- map_dfr(seq_along(desc_vars), function(i) {
  v <- desc_vars[i]
  b <- df_model[[v]][df_model$any_fed == 1]
  n <- df_model[[v]][df_model$any_fed == 0]
  tt <- tryCatch(t.test(b, n), error = function(e) NULL)
  tibble(Variable = desc_labels[i],
         Borrower = sprintf("%.3f (%.3f)", mean(b, na.rm=T), sd(b, na.rm=T)),
         NonBorrower = sprintf("%.3f (%.3f)", mean(n, na.rm=T), sd(n, na.rm=T)),
         Diff = round(mean(b, na.rm=T) - mean(n, na.rm=T), 3),
         p = if (!is.null(tt)) tt$p.value else NA_real_,
         Stars = format_pval(p))
}) %>% mutate(Difference = sprintf("%.3f%s", Diff, Stars))

kbl(desc_tbl %>% select(Variable, Borrower, NonBorrower, Difference),
    format = "html", escape = FALSE,
    col.names = c("Variable", "Borrower Mean (SD)", "Non-Borrower Mean (SD)", "Difference"),
    caption = sprintf("Borrower vs Non-Borrower (N_B=%d, N_NB=%d)",
      sum(df_model$any_fed), sum(df_model$any_fed == 0))) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE) %>%
  footnote(general = "*** p<0.01, ** p<0.05, * p<0.10 (Welch t-test).")
Borrower vs Non-Borrower (N_B=822, N_NB=3404)
Variable Borrower Mean (SD) Non-Borrower Mean (SD) Difference
Log(Assets) 13.771 (1.525) 12.701 (1.367) 1.071***
Total MTM (ℓ) 5.936 (1.978) 5.384 (2.223) 0.551***
Securities MTM (ℓ_S) 0.812 (0.883) 0.649 (0.819) 0.163***
Loan MTM (ℓ_L) 4.984 (1.847) 4.542 (2.074) 0.443***
Book Equity (e) 8.585 (3.110) 9.851 (6.096) -1.266***
Cash/TA 5.251 (5.521) 8.619 (9.161) -3.367***
Loan/Deposit 74.216 (19.811) 70.645 (24.382) 3.572***
Wholesale (%) 1.529 (3.480) 0.868 (3.004) 0.661***
ROA 1.094 (0.586) 1.060 (1.919) 0.034
Uninsured/TA (μ) 27.133 (12.061) 22.944 (11.821) 4.189***
D^U/D 31.816 (14.201) 26.750 (13.805) 5.066***
β^U 0.342 (0.115) 0.329 (0.104) 0.014***
Cost c^U 0.011 (0.000) 0.011 (0.000) 0.000
Franchise (f) 2.953 (1.660) 2.583 (1.588) 0.369***
Uninsured Franchise (f^U) 1.106 (1.103) 0.866 (1.032) 0.240***
MV Equity (E^MV) 5.602 (4.347) 7.050 (6.838) -1.448***
Run Value (v) 4.496 (4.061) 6.184 (6.704) -1.688***
Note:
*** p<0.01, ** p<0.05, * p<0.10 (Welch t-test). 
save_kbl_latex(desc_tbl %>% select(Variable, Borrower, NonBorrower, Difference),
  "Table_Desc_BorrowerVsNon",
  col.names = c("Variable", "Borrower", "Non-Borrower", "Difference"),
  caption = "Borrower vs Non-Borrower: Descriptive Statistics")
## Saved: Table_Desc_BorrowerVsNon.tex

7 EQUITY WATERFALL: Jiang → Franchise → Run Value

Jiang AE says ~825 banks insolvent. Adding franchise value rescues most. F^U shows what’s at stake in a run.

wf <- df_model %>%
  mutate(B = ifelse(any_fed == 1, "Fed Borrower", "Non-Borrower")) %>%
  group_by(B) %>%
  summarise(
    `Book Equity (e)` = mean(book_equity_ratio_raw, na.rm=T),
    `− MTM Loss (ℓ)` = -mean(mtm_total_raw, na.rm=T),
    `= Jiang AE` = mean(adjusted_equity_raw, na.rm=T),
    `+ Franchise (f)` = mean(f_pp, na.rm=T),
    `= E^MV` = mean(emv_pp, na.rm=T),
    `− Unins. Franchise (f^U)` = -mean(f_u_pp, na.rm=T),
    `= Run Value (v)` = mean(v_pp, na.rm=T),
    .groups = "drop") %>%
  pivot_longer(-B, names_to = "Component", values_to = "Value") %>%
  mutate(
    Component = factor(Component, levels = c("Book Equity (e)", "− MTM Loss (ℓ)",
      "= Jiang AE", "+ Franchise (f)", "= E^MV",
      "− Unins. Franchise (f^U)", "= Run Value (v)")),
    Type = case_when(grepl("^=", Component) ~ "Sub", grepl("^−", Component) ~ "Neg", TRUE ~ "Pos"))

p_wf <- ggplot(wf, aes(x = Component, y = Value, fill = Type)) +
  geom_col(alpha = 0.85, width = 0.65) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  geom_text(aes(label = sprintf("%.2f", Value)),
            vjust = ifelse(wf$Value >= 0, -0.3, 1.3), size = 3) +
  facet_wrap(~ B) +
  scale_fill_manual(values = c("Pos" = "#43A047", "Neg" = "#E53935", "Sub" = "#37474F"), guide = "none") +
  scale_x_discrete(labels = function(x) str_wrap(x, 12)) +
  labs(title = "Equity Waterfall: Book → MTM → Franchise → Run Value",
       subtitle = "Mean (pp of TA). v = E^MV − F^U: negative means a run equilibrium is possible.",
       x = NULL, y = "pp of TA") +
  theme_minimal(base_size = 13) +
  theme(plot.title = element_text(face = "bold"), axis.text.x = element_text(size = 9))
print(p_wf)

8 PREDICTION 1: EXTENSIVE MARGIN [Eq. 13]

Theory Eq. (11/13): \(c_3 > 0\). The sensitivity of borrowing probability to MTM losses is increasing in the uninsured franchise \(f^U\). The interaction \(\ell \times f^U\) is THE coordination signature.

8.1 7.1 Main Table: Eq. (13)

controls <- c("ln_assets", "cash_ratio", "loan_to_deposit", "book_equity_ratio", "wholesale", "roa")
ctrl_str <- paste(controls, collapse = " + ")

# Eq. 13: Pr(Borrow) = α + c₁ℓ + c₂f^U + c₃(ℓ×f^U) + γ'X + ε
fml_main <- as.formula(paste0("any_fed ~ mtm_total + uninsured_franchise + mtm_x_f_u + ", ctrl_str))

coef_map <- c("mtm_total" = "ℓ (MTM Loss)", "uninsured_franchise" = "f^U (Uninsured Franchise)",
  "mtm_x_f_u" = "ℓ × f^U (KEY: c₃>0)",
  "ln_assets" = "Log(Assets)", "cash_ratio" = "Cash/TA",
  "loan_to_deposit" = "Loan/Dep", "book_equity_ratio" = "Equity/TA",
  "wholesale" = "Wholesale", "roa" = "ROA")

reg_main <- list(
  "Any Fed"  = feols(fml_main, data = df_model, vcov = "hetero"),
  "BTFP"     = feols(update(fml_main, btfp_crisis ~ .), data = df_model, vcov = "hetero"),
  "DW"       = feols(update(fml_main, dw_crisis ~ .), data = df_model, vcov = "hetero")
)

msummary(reg_main, stars = c("*"=.10, "**"=.05, "***"=.01), gof_omit = "AIC|BIC|Log|RMSE",
  add_rows = dv_summary_rows(reg_main), coef_rename = coef_map,
  title = "Eq. (13): Extensive Margin — Probability of Emergency Borrowing")
Eq. (13): Extensive Margin — Probability of Emergency Borrowing
Any Fed BTFP DW
* p < 0.1, ** p < 0.05, *** p < 0.01
(Intercept) 0.194*** 0.117*** 0.102***
(0.006) (0.005) (0.005)
ℓ (MTM Loss) 0.027*** 0.018*** 0.011**
(0.006) (0.005) (0.005)
f^U (Uninsured Franchise) 0.019*** 0.017*** 0.010*
(0.007) (0.006) (0.006)
ℓ × f^U (KEY: c₃>0) 0.021*** 0.013** 0.017***
(0.006) (0.005) (0.005)
Log(Assets) 0.103*** 0.054*** 0.075***
(0.008) (0.006) (0.007)
Cash/TA -0.021*** -0.026*** 0.001
(0.006) (0.005) (0.005)
Loan/Dep -0.005 -0.010* 0.001
(0.007) (0.006) (0.006)
Equity/TA -0.019*** -0.018*** -0.002
(0.006) (0.005) (0.004)
Wholesale 0.029*** 0.022*** 0.015***
(0.007) (0.006) (0.005)
ROA -0.002 -0.005 -0.002
(0.006) (0.005) (0.005)
Num.Obs. 4226 4226 4226
R2 0.110 0.068 0.074
R2 Adj. 0.109 0.066 0.072
Std.Errors Heteroskedasticity-robust Heteroskedasticity-robust Heteroskedasticity-robust
N (DV = 1) 822 496 430
Mean (DV) 0.1945 0.1174 0.1018 
save_reg_latex(reg_main, "Table_Eq13_Main",
  add_rows = dv_summary_rows(reg_main), coef_rename = coef_map,
  title = "Eq.\\ (13): Extensive Margin")
## Saved: Table_Eq13_Main.tex

8.2 7.2 Progressive Build-Up

Shows the reader what each ingredient adds. Column 1 = Jiang (losses only). Column 2 = add f^U. Column 3 = add the interaction (full theory).

reg_prog <- list(
  "Losses only"   = feols(as.formula(paste0("any_fed ~ mtm_total + ", ctrl_str)),
                          data = df_model, vcov = "hetero"),
  "+ f^U"         = feols(as.formula(paste0("any_fed ~ mtm_total + uninsured_franchise + ", ctrl_str)),
                          data = df_model, vcov = "hetero"),
  "+ ℓ×f^U (Eq.13)" = feols(fml_main, data = df_model, vcov = "hetero")
)

msummary(reg_prog, stars = c("*"=.10, "**"=.05, "***"=.01), gof_omit = "AIC|BIC|Log|RMSE",
  add_rows = dv_summary_rows(reg_prog), coef_rename = coef_map,
  title = "Progressive Build-Up: Jiang → + Uninsured Franchise → + Interaction")
Progressive Build-Up: Jiang → + Uninsured Franchise → + Interaction
Losses only + f^U + ℓ×f^U (Eq.13)
* p < 0.1, ** p < 0.05, *** p < 0.01
(Intercept) 0.192*** 0.192*** 0.194***
(0.006) (0.006) (0.006)
ℓ (MTM Loss) 0.022*** 0.024*** 0.027***
(0.006) (0.006) (0.006)
Log(Assets) 0.109*** 0.104*** 0.103***
(0.007) (0.008) (0.008)
Cash/TA -0.022*** -0.023*** -0.021***
(0.006) (0.006) (0.006)
Loan/Dep -0.010 -0.007 -0.005
(0.007) (0.007) (0.007)
Equity/TA -0.021*** -0.021*** -0.019***
(0.006) (0.006) (0.006)
Wholesale 0.028*** 0.029*** 0.029***
(0.007) (0.007) (0.007)
ROA 0.000 -0.001 -0.002
(0.006) (0.006) (0.006)
f^U (Uninsured Franchise) 0.013* 0.019***
(0.007) (0.007)
ℓ × f^U (KEY: c₃>0) 0.021***
(0.006)
Num.Obs. 4226 4226 4226
R2 0.107 0.108 0.110
R2 Adj. 0.105 0.106 0.109
Std.Errors Heteroskedasticity-robust Heteroskedasticity-robust Heteroskedasticity-robust
N (DV = 1) 822 822 822
Mean (DV) 0.1945 0.1945 0.1945 
save_reg_latex(reg_prog, "Table_Eq13_Progressive",
  add_rows = dv_summary_rows(reg_prog), coef_rename = coef_map,
  title = "Progressive Build-Up: Losses $\\to$ $+f^U$ $\\to$ $+\\ell \\times f^U$")
## Saved: Table_Eq13_Progressive.tex

8.3 7.3 Observable vs Unobservable Signal

Theory Eq. (14): Securities losses ℓ_S are publicly observable and determine BTFP capacity. Loan losses ℓ_L are unobservable. Coordination should amplify through ℓ_S × f^U, not ℓ_L × f^U.

df_signal <- df_model %>%
  mutate(sec_x_f_u  = mtm_sec * uninsured_franchise,
         loan_x_f_u = mtm_loan * uninsured_franchise)

fml_signal <- as.formula(paste0(
  "any_fed ~ mtm_sec + mtm_loan + uninsured_franchise + sec_x_f_u + loan_x_f_u + ", ctrl_str))

reg_signal <- list(
  "Any Fed" = feols(fml_signal, data = df_signal, vcov = "hetero"),
  "BTFP"    = feols(update(fml_signal, btfp_crisis ~ .), data = df_signal, vcov = "hetero"),
  "DW"      = feols(update(fml_signal, dw_crisis ~ .), data = df_signal, vcov = "hetero")
)

msummary(reg_signal, stars = c("*"=.10, "**"=.05, "***"=.01), gof_omit = "AIC|BIC|Log|RMSE",
  add_rows = dv_summary_rows(reg_signal),
  coef_rename = c("mtm_sec" = "ℓ_S (Securities)", "mtm_loan" = "ℓ_L (Loans)",
    "uninsured_franchise" = "f^U", "sec_x_f_u" = "ℓ_S × f^U (observable; expect >0)",
    "loan_x_f_u" = "ℓ_L × f^U (unobservable; expect ≈0)"),
  title = "Observable vs Unobservable Signal: ℓ_S × f^U vs ℓ_L × f^U")
Observable vs Unobservable Signal: ℓ_S × f^U vs ℓ_L × f^U
Any Fed BTFP DW
* p < 0.1, ** p < 0.05, *** p < 0.01
(Intercept) 0.195*** 0.117*** 0.103***
(0.006) (0.005) (0.005)
ℓ_S (Securities) 0.029*** 0.022*** 0.010
(0.009) (0.008) (0.007)
ℓ_L (Loans) 0.022*** 0.011** 0.009
(0.007) (0.006) (0.006)
f^U 0.020*** 0.017*** 0.011*
(0.007) (0.006) (0.006)
ℓ_S × f^U (observable; expect >0) 0.012** 0.009* 0.011**
(0.006) (0.005) (0.005)
ℓ_L × f^U (unobservable; expect ≈0) 0.023*** 0.010* 0.016***
(0.006) (0.006) (0.006)
ln_assets 0.102*** 0.054*** 0.075***
(0.008) (0.006) (0.007)
cash_ratio -0.018** -0.022*** 0.002
(0.007) (0.005) (0.006)
loan_to_deposit 0.001 -0.001 0.002
(0.010) (0.009) (0.008)
book_equity_ratio -0.018*** -0.017*** -0.002
(0.006) (0.005) (0.004)
wholesale 0.028*** 0.021*** 0.015***
(0.007) (0.006) (0.005)
roa -0.002 -0.005 -0.002
(0.006) (0.005) (0.005)
Num.Obs. 4226 4226 4226
R2 0.111 0.069 0.075
R2 Adj. 0.109 0.066 0.072
Std.Errors Heteroskedasticity-robust Heteroskedasticity-robust Heteroskedasticity-robust
N (DV = 1) 822 496 430
Mean (DV) 0.1945 0.1174 0.1018 
save_reg_latex(reg_signal, "Table_Eq13_Signal",
  add_rows = dv_summary_rows(reg_signal),
  coef_rename = c("mtm_sec"="$\\ell_S$", "mtm_loan"="$\\ell_L$",
    "uninsured_franchise"="$f^U$",
    "sec_x_f_u"="$\\ell_S \\times f^U$", "loan_x_f_u"="$\\ell_L \\times f^U$"),
  title = "Observable vs.\\ Unobservable Signal")
## Saved: Table_Eq13_Signal.tex

8.4 7.4 By Facility: ℓ_S vs ℓ_L

Theory Eq. (14): Securities losses should predict BTFP (par channel), not DW.

fml_sep <- function(dv) as.formula(paste0(dv, " ~ mtm_sec + mtm_loan + uninsured_franchise + ", ctrl_str))

reg_facility <- list(
  "BTFP"  = feols(fml_sep("btfp_crisis"), data = df_model, vcov = "hetero"),
  "DW"    = feols(fml_sep("dw_crisis"),   data = df_model, vcov = "hetero"),
  "FHLB"  = feols(fml_sep("fhlb_user"),   data = df_model, vcov = "hetero")
)

msummary(reg_facility, stars = c("*"=.10, "**"=.05, "***"=.01), gof_omit = "AIC|BIC|Log|RMSE",
  add_rows = dv_summary_rows(reg_facility),
  coef_rename = c("mtm_sec" = "ℓ_S (Securities)", "mtm_loan" = "ℓ_L (Loans)",
    "uninsured_franchise" = "f^U"),
  title = "Par-Value Channel: ℓ_S Should Predict BTFP, Not DW or FHLB (Eq. 14)")
Par-Value Channel: ℓ_S Should Predict BTFP, Not DW or FHLB (Eq. 14)
BTFP DW FHLB
* p < 0.1, ** p < 0.05, *** p < 0.01
(Intercept) 0.115*** 0.100*** 0.070***
(0.005) (0.004) (0.004)
ℓ_S (Securities) 0.022*** 0.010 -0.008
(0.008) (0.007) (0.006)
ℓ_L (Loans) 0.008 0.005 0.007
(0.005) (0.005) (0.005)
f^U 0.013** 0.005 -0.002
(0.006) (0.006) (0.004)
ln_assets 0.054*** 0.076*** 0.016***
(0.006) (0.007) (0.005)
cash_ratio -0.022*** 0.001 -0.020***
(0.005) (0.006) (0.004)
loan_to_deposit -0.001 0.002 0.013*
(0.009) (0.008) (0.007)
book_equity_ratio -0.018*** -0.004 0.008*
(0.004) (0.004) (0.004)
wholesale 0.021*** 0.015*** -0.001
(0.006) (0.005) (0.004)
roa -0.005 -0.001 -0.007*
(0.005) (0.005) (0.004)
Num.Obs. 4226 4226 4226
R2 0.067 0.072 0.025
R2 Adj. 0.065 0.070 0.023
Std.Errors Heteroskedasticity-robust Heteroskedasticity-robust Heteroskedasticity-robust
N (DV = 1) 496 430 299
Mean (DV) 0.1174 0.1018 0.0708 
save_reg_latex(reg_facility, "Table_Eq14_Facility",
  add_rows = dv_summary_rows(reg_facility),
  coef_rename = c("mtm_sec"="$\\ell_S$", "mtm_loan"="$\\ell_L$", "uninsured_franchise"="$f^U$"),
  title = "Par-Value Channel: $\\ell_S$ Predicts BTFP (Eq.\\ 14)")
## Saved: Table_Eq14_Facility.tex

9 PREDICTION 2: INTENSIVE MARGIN [Eq. 15]

Theory Eq. (15): Among borrowers, amounts should increase in securities losses ℓ_S (which determine BTFP capacity Δ_BTFP), but not necessarily in total losses or loan losses.

9.1 8.1 Amount Regressions

# Among Fed borrowers only
df_borr <- df_model %>%
  filter(any_fed == 1) %>%
  mutate(
    # Log amount (in thousands)
    ln_fed_amt = log(fed_amt + 1),
    ln_btfp_amt = log(replace_na(btfp_crisis_amt, 0) + 1),
    # Normalized: amount / total assets
    fed_amt_ta = fed_amt / (total_asset * 1000),  # total_asset likely in thousands
    btfp_amt_ta = replace_na(btfp_crisis_amt, 0) / (total_asset * 1000)
  )

cat(sprintf("Borrowers: %d | Mean Fed amount: $%.0fM | Mean BTFP amount: $%.0fM\n",
    nrow(df_borr),
    mean(df_borr$fed_amt, na.rm=T) / 1e6,
    mean(df_borr$btfp_crisis_amt, na.rm=T) / 1e6))
## Borrowers: 822 | Mean Fed amount: $813M | Mean BTFP amount: $263M
if (nrow(df_borr) >= 30) {
  # Eq. 15: Amount = α + d₁ ℓ_S + γ'X + ε
  fml_amt <- as.formula(paste0("ln_fed_amt ~ mtm_sec + mtm_loan + uninsured_franchise + ", ctrl_str))
  fml_btfp_amt <- as.formula(paste0("ln_btfp_amt ~ mtm_sec + mtm_loan + uninsured_franchise + ", ctrl_str))

  reg_intensive <- list(
    "ln(Fed Amount)"  = feols(fml_amt, data = df_borr, vcov = "hetero"),
    "ln(BTFP Amount)" = feols(fml_btfp_amt, data = df_borr %>% filter(btfp_crisis == 1), vcov = "hetero")
  )

  msummary(reg_intensive, stars = c("*"=.10, "**"=.05, "***"=.01), gof_omit = "AIC|BIC|Log|RMSE",
    coef_rename = c("mtm_sec" = "ℓ_S (Securities; d₁>0)",
      "mtm_loan" = "ℓ_L (Loans)", "uninsured_franchise" = "f^U"),
    title = "Eq. (15): Intensive Margin — Amount Borrowed (Among Borrowers Only)",
    notes = "Theory: d₁ > 0. Securities losses determine BTFP capacity (Δ_BTFP).")

  save_reg_latex(reg_intensive, "Table_Eq15_Intensive",
    coef_rename = c("mtm_sec"="$\\ell_S$ ($d_1>0$)", "mtm_loan"="$\\ell_L$",
      "uninsured_franchise"="$f^U$"),
    title = "Eq.\\ (15): Intensive Margin --- Amount Borrowed")
}
## Saved: Table_Eq15_Intensive.tex

10 FALSIFICATION TESTS

10.1 9.1 FHLB Falsification

Theory scope condition: FHLB borrowing is routine, not run-driven. \(c_3 \approx 0\).

fml_fhlb <- as.formula(paste0("fhlb_user ~ mtm_total + uninsured_franchise + mtm_x_f_u + ", ctrl_str))

reg_falsif <- list(
  "Any Fed (crisis)" = feols(fml_main, data = df_model, vcov = "hetero"),
  "FHLB (crisis)"    = feols(fml_fhlb, data = df_model, vcov = "hetero")
)

msummary(reg_falsif, stars = c("*"=.10, "**"=.05, "***"=.01), gof_omit = "AIC|BIC|Log|RMSE",
  add_rows = dv_summary_rows(reg_falsif), coef_rename = coef_map,
  title = "Falsification: FHLB Should NOT Show ℓ × f^U > 0",
  notes = "FHLB = routine advance. Significant ℓ × f^U for FHLB would undermine coordination interpretation.")
Falsification: FHLB Should NOT Show ℓ × f^U > 0
Any Fed (crisis) FHLB (crisis)
* p < 0.1, ** p < 0.05, *** p < 0.01
FHLB = routine advance. Significant ℓ × f^U for FHLB would undermine coordination interpretation.
(Intercept) 0.194*** 0.069***
(0.006) (0.004)
ℓ (MTM Loss) 0.027*** -0.000
(0.006) (0.005)
f^U (Uninsured Franchise) 0.019*** -0.004
(0.007) (0.004)
ℓ × f^U (KEY: c₃>0) 0.021*** -0.007**
(0.006) (0.003)
Log(Assets) 0.103*** 0.016***
(0.008) (0.005)
Cash/TA -0.021*** -0.016***
(0.006) (0.004)
Loan/Dep -0.005 0.022***
(0.007) (0.005)
Equity/TA -0.019*** 0.007*
(0.006) (0.004)
Wholesale 0.029*** -0.002
(0.007) (0.004)
ROA -0.002 -0.007*
(0.006) (0.004)
Num.Obs. 4226 4226
R2 0.110 0.025
R2 Adj. 0.109 0.023
Std.Errors Heteroskedasticity-robust Heteroskedasticity-robust
N (DV = 1) 822 299
Mean (DV) 0.1945 0.0708 
save_reg_latex(reg_falsif, "Table_Falsif_FHLB",
  add_rows = dv_summary_rows(reg_falsif), coef_rename = coef_map,
  title = "Falsification: FHLB vs.\\ Fed Borrowing")
## Saved: Table_Falsif_FHLB.tex

10.2 9.2 Temporal Falsification: Arb Period

Theory scope condition: After Nov 2023, BTFP rate fell below IOR → carry trade. \(c_3 \approx 0\).

df_arb_model <- df_arb_clean %>% filter(!is.na(f_u_pp))

df_arb_reg <- df_arb_model %>%
  mutate(
    mtm_total_a = standardize_z(winsorize(mtm_total_raw)),
    f_u_a       = standardize_z(winsorize(f_u_pp)),
    mtm_x_f_u_a = mtm_total_a * f_u_a,
    ln_assets_a = standardize_z(winsorize(ln_assets_raw)),
    cash_ratio_a = standardize_z(winsorize(cash_ratio_raw)),
    loan_to_deposit_a = standardize_z(winsorize(loan_to_deposit_raw)),
    book_equity_ratio_a = standardize_z(winsorize(book_equity_ratio_raw)),
    wholesale_a = standardize_z(winsorize(wholesale_raw)),
    roa_a = standardize_z(winsorize(roa_raw))
  )

arb_ctrl <- c("ln_assets_a", "cash_ratio_a", "loan_to_deposit_a",
              "book_equity_ratio_a", "wholesale_a", "roa_a")
fml_arb <- as.formula(paste0("any_fed ~ mtm_total_a + f_u_a + mtm_x_f_u_a + ",
                              paste(arb_ctrl, collapse = " + ")))

reg_arb <- list(
  "Crisis: Any Fed" = feols(fml_main, data = df_model, vcov = "hetero")
)
if (sum(df_arb_reg$any_fed) >= 10 & sum(df_arb_reg$any_fed == 0) >= 10)
  reg_arb[["Arb: Any Fed"]] <- feols(fml_arb, data = df_arb_reg, vcov = "hetero")

arb_map <- c("mtm_total" = "ℓ", "uninsured_franchise" = "f^U", "mtm_x_f_u" = "ℓ × f^U",
  "mtm_total_a" = "ℓ", "f_u_a" = "f^U", "mtm_x_f_u_a" = "ℓ × f^U",
  "ln_assets" = "Log(Assets)", "ln_assets_a" = "Log(Assets)",
  "cash_ratio" = "Cash/TA", "cash_ratio_a" = "Cash/TA",
  "loan_to_deposit" = "Loan/Dep", "loan_to_deposit_a" = "Loan/Dep",
  "book_equity_ratio" = "Equity/TA", "book_equity_ratio_a" = "Equity/TA",
  "wholesale" = "Wholesale", "wholesale_a" = "Wholesale",
  "roa" = "ROA", "roa_a" = "ROA")

msummary(reg_arb, stars = c("*"=.10, "**"=.05, "***"=.01), gof_omit = "AIC|BIC|Log|RMSE",
  coef_rename = arb_map,
  title = "Temporal Falsification: Crisis vs Arb Period",
  notes = "ℓ × f^U significant in crisis, ≈ 0 in arb period (carry trade, not coordination).")
Temporal Falsification: Crisis vs Arb Period
Crisis: Any Fed Arb: Any Fed
* p < 0.1, ** p < 0.05, *** p < 0.01
ℓ × f^U significant in crisis, ≈ 0 in arb period (carry trade, not coordination).
(Intercept) 0.194*** 0.241***
(0.006) (0.006)
ℓ 0.027*** 0.014*
(0.006) (0.007)
f^U 0.019*** 0.007
(0.007) (0.007)
ℓ × f^U 0.021*** 0.006
(0.006) (0.006)
Log(Assets) 0.103*** 0.092***
(0.008) (0.008)
Cash/TA -0.021*** -0.033***
(0.006) (0.007)
Loan/Dep -0.005 -0.001
(0.007) (0.007)
Equity/TA -0.019*** -0.012**
(0.006) (0.006)
Wholesale 0.029*** 0.099***
(0.007) (0.008)
ROA -0.002 -0.024***
(0.006) (0.006)
Num.Obs. 4226 4151
R2 0.110 0.150
R2 Adj. 0.109 0.148
Std.Errors Heteroskedasticity-robust Heteroskedasticity-robust 
save_reg_latex(reg_arb, "Table_Falsif_Arb",
  coef_rename = arb_map,
  title = "Temporal Falsification: Crisis vs.\\ Arb Period")
## Saved: Table_Falsif_Arb.tex

11 ROBUSTNESS

11.1 10.1 Logit

fml_logit <- as.formula(paste0("any_fed ~ mtm_total + uninsured_franchise + mtm_x_f_u + ", ctrl_str))

reg_logit <- list(
  "LPM (baseline)" = feols(fml_main, data = df_model, vcov = "hetero"),
  "Logit"          = feglm(fml_logit, data = df_model, family = binomial("logit"), vcov = "hetero")
)

msummary(reg_logit, stars = c("*"=.10, "**"=.05, "***"=.01), gof_omit = "AIC|BIC|Log|RMSE",
  coef_rename = coef_map,
  title = "Robustness: LPM vs Logit")
Robustness: LPM vs Logit
LPM (baseline) Logit
* p < 0.1, ** p < 0.05, *** p < 0.01
(Intercept) 0.194*** -1.684***
(0.006) (0.048)
ℓ (MTM Loss) 0.027*** 0.188***
(0.006) (0.048)
f^U (Uninsured Franchise) 0.019*** 0.109**
(0.007) (0.045)
ℓ × f^U (KEY: c₃>0) 0.021*** 0.085**
(0.006) (0.042)
Log(Assets) 0.103*** 0.654***
(0.008) (0.048)
Cash/TA -0.021*** -0.268***
(0.006) (0.071)
Loan/Dep -0.005 0.020
(0.007) (0.055)
Equity/TA -0.019*** -0.257***
(0.006) (0.061)
Wholesale 0.029*** 0.159***
(0.007) (0.039)
ROA -0.002 0.033
(0.006) (0.048)
Num.Obs. 4226 4226
R2 0.110 0.117
R2 Adj. 0.109 0.112
Std.Errors Heteroskedasticity-robust Heteroskedasticity-robust 
save_reg_latex(reg_logit, "Table_Robust_Logit",
  coef_rename = coef_map,
  title = "Robustness: LPM vs.\\ Logit")
## Saved: Table_Robust_Logit.tex

11.2 10.2 Decay Rate Sensitivity

DSSW baseline: δ = 10%. Lower δ → slower deposit decay → higher f and f^U.

delta_grid <- c(0.05, 0.10, 0.15)
reg_delta <- list()

for (dd in delta_grid) {
  cf_dd <- 1 / (y_10yr + dd)
  lab <- sprintf("δ=%.0f%% (1/(y+δ)=%.2f)%s", dd*100, cf_dd,
                 ifelse(dd == delta_decay, " [baseline]", ""))

  df_d <- df_model %>%
    mutate(
      c_u = ifelse(!is.na(cost_u_raw), cost_u_raw, 0),
      f_d = pmax(((1 - beta_u_clipped) * y_10yr - c_u) * cf_dd * mu_decimal, 0) * 100,
      f_u_d = uninsured_share_d * f_d,
      f_u_d_w = winsorize(f_u_d),
      fu_d = standardize_z(f_u_d_w),
      mtm_x_fu_d = mtm_total * fu_d
    )

  fml_d <- as.formula(paste0("any_fed ~ mtm_total + fu_d + mtm_x_fu_d + ", ctrl_str))
  reg_delta[[lab]] <- feols(fml_d, data = df_d, vcov = "hetero")
}

msummary(reg_delta, stars = c("*"=.10, "**"=.05, "***"=.01), gof_omit = "AIC|BIC|Log|RMSE",
  coef_rename = c("mtm_total" = "ℓ", "fu_d" = "f^U", "mtm_x_fu_d" = "ℓ × f^U"),
  title = "Robustness: Eq. (13) Under Alternative Deposit Decay Rates",
  notes = "DSSW baseline: δ=10%. ℓ × f^U should be stable across decay rates.")
Robustness: Eq. (13) Under Alternative Deposit Decay Rates
δ=5% (1/(y+δ)=11.21) δ=10% (1/(y+δ)=7.18) [baseline] δ=15% (1/(y+δ)=5.29)
* p < 0.1, ** p < 0.05, *** p < 0.01
DSSW baseline: δ=10%. ℓ × f^U should be stable across decay rates.
(Intercept) 0.194*** 0.194*** 0.194***
(0.006) (0.006) (0.006)
ℓ 0.027*** 0.027*** 0.027***
(0.006) (0.006) (0.006)
f^U 0.019*** 0.019*** 0.019***
(0.007) (0.007) (0.007)
ℓ × f^U 0.021*** 0.021*** 0.021***
(0.006) (0.006) (0.006)
ln_assets 0.103*** 0.103*** 0.103***
(0.008) (0.008) (0.008)
cash_ratio -0.021*** -0.021*** -0.021***
(0.006) (0.006) (0.006)
loan_to_deposit -0.005 -0.005 -0.005
(0.007) (0.007) (0.007)
book_equity_ratio -0.019*** -0.019*** -0.019***
(0.006) (0.006) (0.006)
wholesale 0.029*** 0.029*** 0.029***
(0.007) (0.007) (0.007)
roa -0.002 -0.002 -0.002
(0.006) (0.006) (0.006)
Num.Obs. 4226 4226 4226
R2 0.110 0.110 0.110
R2 Adj. 0.109 0.109 0.109
Std.Errors Heteroskedasticity-robust Heteroskedasticity-robust Heteroskedasticity-robust 
save_reg_latex(reg_delta, "Table_Robust_Delta",
  coef_rename = c("mtm_total"="$\\ell$", "fu_d"="$f^U$", "mtm_x_fu_d"="$\\ell \\times f^U$"),
  title = "Robustness: Alternative Deposit Decay Rates")
## Saved: Table_Robust_Delta.tex
delta_grid <- c(0.05, 0.10, 0.15)
reg_delta <- list()

for (dd in delta_grid) {
  cf_dd <- 1 / (y_10yr + dd)
  lab <- sprintf("δ=%.0f%% (1/(y+δ)=%.2f)%s", dd*100, cf_dd,
                 ifelse(dd == delta_decay, " [baseline]", ""))

  df_d <- df_model %>%
    mutate(
      c_u = ifelse(!is.na(cost_u_raw), cost_u_raw, 0),
      f_d = pmax(((1 - beta_u_clipped) * y_10yr - c_u) * cf_dd * mu_decimal, 0) * 100,
      f_u_d = uninsured_share_d * f_d,
      f_u_d_w = winsorize(f_u_d),
      
      # MEAN-CENTER THE RAW VARIABLES
      mtm_c = mtm_total_w - mean(mtm_total_w, na.rm = TRUE),
      f_u_c = f_u_d_w - mean(f_u_d_w, na.rm = TRUE),
      
      # INTERACT THE CENTERED VARIABLES
      mtm_x_fu_c = mtm_c * f_u_c
    )

  # Use the centered variables in the regression
  fml_d <- as.formula(paste0("any_fed ~ mtm_c + f_u_c + mtm_x_fu_c + ", ctrl_str))
  reg_delta[[lab]] <- feols(fml_d, data = df_d, vcov = "hetero")
}

# 1. Print to RMarkdown HTML document
msummary(reg_delta, stars = c("*"=.10, "**"=.05, "***"=.01), gof_omit = "AIC|BIC|Log|RMSE",
  coef_rename = c("mtm_c" = "ℓ (Centered)", "f_u_c" = "f^U (Centered)", "mtm_x_fu_c" = "ℓ × f^U (Centered)"),
  title = "Robustness: Eq. (13) Under Alternative Deposit Decay Rates (Mean-Centered)",
  notes = "Explanatory variables are mean-centered raw values. Controls are Z-standardized.")
Robustness: Eq. (13) Under Alternative Deposit Decay Rates (Mean-Centered)
δ=5% (1/(y+δ)=11.21) δ=10% (1/(y+δ)=7.18) [baseline] δ=15% (1/(y+δ)=5.29)
* p < 0.1, ** p < 0.05, *** p < 0.01
Explanatory variables are mean-centered raw values. Controls are Z-standardized.
(Intercept) 0.195*** 0.195*** 0.195***
(0.006) (0.006) (0.006)
ℓ (Centered) 0.013*** 0.013*** 0.013***
(0.003) (0.003) (0.003)
f^U (Centered) 0.015*** 0.023*** 0.031***
(0.006) (0.009) (0.012)
ℓ × f^U (Centered) 0.008*** 0.012*** 0.016***
(0.002) (0.003) (0.004)
ln_assets 0.103*** 0.103*** 0.103***
(0.008) (0.008) (0.008)
cash_ratio -0.021*** -0.021*** -0.021***
(0.006) (0.006) (0.006)
loan_to_deposit -0.005 -0.005 -0.005
(0.007) (0.007) (0.007)
book_equity_ratio -0.019*** -0.019*** -0.019***
(0.006) (0.006) (0.006)
wholesale 0.029*** 0.029*** 0.029***
(0.007) (0.007) (0.007)
roa -0.002 -0.002 -0.002
(0.006) (0.006) (0.006)
Num.Obs. 4226 4226 4226
R2 0.110 0.110 0.110
R2 Adj. 0.109 0.109 0.109
Std.Errors Heteroskedasticity-robust Heteroskedasticity-robust Heteroskedasticity-robust 
# 2. Save to LaTeX folder
save_reg_latex(reg_delta, "Table_Robust_Delta_Centered",
  coef_rename = c("mtm_c" = "$\\ell$ (Centered)", "f_u_c" = "$f^U$ (Centered)", "mtm_x_fu_c" = "$\\ell \\times f^U$ (Centered)"),
  title = "Robustness: Eq. (13) Under Alternative Deposit Decay Rates (Mean-Centered)",
  notes = "DSSW baseline: $\\delta=10\\%$. Explanatory variables are mean-centered raw values to preserve magnitude shifts. Controls are Z-standardized.")
## Saved: Table_Robust_Delta_Centered.tex

12 PLOTS

borr_colors <- c("Fed Borrower" = "#C62828", "Non-Borrower" = "#1565C0")
theme_gp <- theme_minimal(base_size = 13) +
  theme(plot.title = element_text(face = "bold", size = 14),
        plot.subtitle = element_text(color = "grey40", size = 11),
        plot.caption = element_text(color = "grey50", size = 9, hjust = 0),
        legend.position = "bottom", panel.grid.minor = element_blank())

12.1 Scatter: ℓ vs f^U

p1 <- ggplot(df_model, aes(x = mtm_total_raw, y = f_u_pp)) +
  geom_point(aes(color = ifelse(any_fed == 1, "Fed Borrower", "Non-Borrower"),
                 shape = ifelse(any_fed == 1, "Fed Borrower", "Non-Borrower")),
             alpha = 0.5, size = 1.8) +
  scale_color_manual(values = borr_colors, name = NULL) +
  scale_shape_manual(values = c("Fed Borrower" = 17, "Non-Borrower" = 1), name = NULL) +
  labs(title = "MTM Losses vs Uninsured Franchise Value",
       subtitle = "Triangles = Fed borrowers. Theory: borrowing concentrates in the upper-right (high ℓ, high f^U)",
       x = "Total MTM Loss / TA (pp)", y = "Uninsured Franchise f^U / TA (pp)") +
  theme_gp
print(p1)

12.2 Run Value Distribution

p2 <- ggplot(df_model %>% mutate(B = ifelse(any_fed == 1, "Fed Borrower", "Non-Borrower")),
       aes(x = v_pp, fill = B, color = B)) +
  geom_density(alpha = 0.3, linewidth = 0.7) +
  geom_vline(xintercept = 0, linetype = "dashed", linewidth = 0.8) +
  scale_fill_manual(values = borr_colors) + scale_color_manual(values = borr_colors) +
  labs(title = "Run Value (v = E^MV − F^U): Borrowers vs Non-Borrowers",
       subtitle = "Dashed line: v=0. Left of line: run equilibrium possible (F^U > E^MV).",
       x = "Run Value v / TA (pp)", y = "Density") +
  theme_gp
print(p2)

13 New Date and period classifications

14 NEW REGRESSIONS: FINE-GRAINED DATES & SIZE SUBSAMPLES

14.1 Data Preparation for New Timeframes & Size

# 1. Create specific date indicators from the raw loan data
dw_m1_7  <- create_borrower_indicator(dw_loans, "dw_loan_date", "rssd_id", "dw_loan_amount", 
                                      as.Date("2023-03-01"), as.Date("2023-03-07"), "dw_m1_7")
dw_m8_12 <- create_borrower_indicator(dw_loans, "dw_loan_date", "rssd_id", "dw_loan_amount", 
                                      as.Date("2023-03-08"), as.Date("2023-03-12"), "dw_m8_12")

btfp_r1  <- create_borrower_indicator(btfp_loans, "btfp_loan_date", "rssd_id", "btfp_loan_amount", 
                                      as.Date("2023-03-13"), as.Date("2023-04-28"), "btfp_r1")
dw_r1    <- create_borrower_indicator(dw_loans, "dw_loan_date", "rssd_id", "dw_loan_amount", 
                                      as.Date("2023-03-13"), as.Date("2023-04-28"), "dw_r1")

btfp_r2  <- create_borrower_indicator(btfp_loans, "btfp_loan_date", "rssd_id", "btfp_loan_amount", 
                                      as.Date("2023-04-29"), as.Date("2023-05-04"), "btfp_r2")
dw_r2    <- create_borrower_indicator(dw_loans, "dw_loan_date", "rssd_id", "dw_loan_amount", 
                                      as.Date("2023-04-29"), as.Date("2023-05-04"), "dw_r2")

# 2. Join these new indicators into our main df_model
# Note: total_asset in Call Reports is in thousands. 
# So 1 Billion = 1,000,000 thousands; 10 Billion = 10,000,000 thousands.
df_model_extended <- df_model %>%
  left_join(dw_m1_7 %>% select(idrssd, dw_m1_7), by = "idrssd") %>%
  left_join(dw_m8_12 %>% select(idrssd, dw_m8_12), by = "idrssd") %>%
  left_join(btfp_r1 %>% select(idrssd, btfp_r1), by = "idrssd") %>%
  left_join(dw_r1 %>% select(idrssd, dw_r1), by = "idrssd") %>%
  left_join(btfp_r2 %>% select(idrssd, btfp_r2), by = "idrssd") %>%
  left_join(dw_r2 %>% select(idrssd, dw_r2), by = "idrssd") %>%
  mutate(
    dw_m1_7    = replace_na(dw_m1_7, 0L),
    dw_m8_12   = replace_na(dw_m8_12, 0L),
    btfp_r1    = replace_na(btfp_r1, 0L),
    dw_r1      = replace_na(dw_r1, 0L),
    any_fed_r1 = as.integer(btfp_r1 == 1 | dw_r1 == 1),
    btfp_r2    = replace_na(btfp_r2, 0L),
    dw_r2      = replace_na(dw_r2, 0L),
    any_fed_r2 = as.integer(btfp_r2 == 1 | dw_r2 == 1),
    
    # Create Size Buckets
    size_bucket = case_when(
      total_asset < 1000000 ~ "< 1B",
      total_asset >= 1000000 & total_asset < 10000000 ~ "1B - 10B",
      total_asset >= 10000000 ~ "> 10B"
    )
  )

# Base formula RHS for all the new models
fml_rhs_base <- paste0("mtm_total + uninsured_franchise + mtm_x_f_u + ", ctrl_str)

15 Early DW Usage: Mar 1-7 vs Mar 8-12

reg_dw_early <- list(
  "DW (Mar 1-7)"  = feols(as.formula(paste("dw_m1_7 ~", fml_rhs_base)), data = df_model_extended, vcov = "hetero"),
  "DW (Mar 8-12)" = feols(as.formula(paste("dw_m8_12 ~", fml_rhs_base)), data = df_model_extended, vcov = "hetero")
)

msummary(reg_dw_early, stars = c("*"=.10, "**"=.05, "***"=.01), gof_omit = "AIC|BIC|Log|RMSE",
  add_rows = dv_summary_rows(reg_dw_early), coef_rename = coef_map,
  title = "Table DW: DW Usage Pre- and Early-Crisis")
Table DW: DW Usage Pre- and Early-Crisis
DW (Mar 1-7) DW (Mar 8-12)
* p < 0.1, ** p < 0.05, *** p < 0.01
(Intercept) 0.016*** 0.016***
(0.002) (0.002)
ℓ (MTM Loss) 0.001 0.000
(0.002) (0.002)
f^U (Uninsured Franchise) 0.001 -0.000
(0.003) (0.003)
ℓ × f^U (KEY: c₃>0) 0.003 0.004*
(0.002) (0.002)
Log(Assets) 0.012*** 0.016***
(0.003) (0.003)
Cash/TA -0.002 -0.002
(0.002) (0.002)
Loan/Dep -0.001 -0.001
(0.002) (0.002)
Equity/TA 0.001 -0.002
(0.002) (0.002)
Wholesale 0.011*** 0.010***
(0.003) (0.003)
ROA -0.002 -0.003*
(0.002) (0.002)
Num.Obs. 4226 4226
R2 0.021 0.027
R2 Adj. 0.019 0.025
Std.Errors Heteroskedasticity-robust Heteroskedasticity-robust
N (DV = 1) 67 66
Mean (DV) 0.0159 0.0156 
save_reg_latex(reg_dw_early, "Table_DW_Early",
  add_rows = dv_summary_rows(reg_dw_early), coef_rename = coef_map,
  title = "Table DW: DW Usage Pre- and Early-Crisis")
## Saved: Table_DW_Early.tex
# ==============================================================================
# 1. DESCRIPTIVE SPLIT: What characteristics define a Repeat Borrower?
# Comparing Means: Single vs Repeat (March 8-12 panic window)
# ==============================================================================

# 1. Calculate the exact number of DW loans per bank for the early periods
dw_counts_m1_7 <- dw_loans %>%
  filter(dw_loan_date >= as.Date("2023-03-01") & dw_loan_date <= as.Date("2023-03-07")) %>%
  group_by(rssd_id) %>%
  summarise(dw_m1_7_count = n(), .groups = "drop") %>%
  rename(idrssd = rssd_id)

dw_counts_m8_12 <- dw_loans %>%
  filter(dw_loan_date >= as.Date("2023-03-08") & dw_loan_date <= as.Date("2023-03-12")) %>%
  group_by(rssd_id) %>%
  summarise(dw_m8_12_count = n(), .groups = "drop") %>%
  rename(idrssd = rssd_id)

# 2. Merge counts AND amounts into the main dataset, create Repeat Dummies
df_model_extended <- df_model_extended %>%
  select(-any_of(c("dw_m1_7_count", "dw_m8_12_count", "dw_m1_7_repeat", "dw_m8_12_repeat", "dw_m1_7_amt", "dw_m8_12_amt"))) %>%
  left_join(dw_counts_m1_7, by = "idrssd") %>%
  left_join(dw_counts_m8_12, by = "idrssd") %>%
  left_join(dw_m1_7 %>% select(idrssd, dw_m1_7_amt), by = "idrssd") %>%
  left_join(dw_m8_12 %>% select(idrssd, dw_m8_12_amt), by = "idrssd") %>%
  mutate(
    dw_m1_7_count = replace_na(dw_m1_7_count, 0L),
    dw_m8_12_count = replace_na(dw_m8_12_count, 0L),
    
    # Dummies for Repeat Borrowers
    dw_m1_7_repeat = as.integer(dw_m1_7_count > 1),
    dw_m8_12_repeat = as.integer(dw_m8_12_count > 1),
    
    # Log Amounts
    ln_dw_amt_m1_7 = log(replace_na(dw_m1_7_amt, 0) + 1),
    ln_dw_amt_m8_12 = log(replace_na(dw_m8_12_amt, 0) + 1)
  )

# 3. Create subsets for BORROWERS ONLY in each period
df_borr_1_7 <- df_model_extended %>% filter(dw_m1_7 == 1)
df_borr_8_12 <- df_model_extended %>% filter(dw_m8_12 == 1)

# 4. Define variables: Theory variables (to interact) vs Controls (standalone)
theory_vars <- c("mtm_total", "uninsured_franchise", "mtm_x_f_u")
# Note: ctrl_str is already defined globally in your script as "ln_assets + cash_ratio + ..."

# 5. Define the specific formulas: (Theory) * repeat_dummy + Controls
fml_interact_1_7 <- as.formula(
  paste0("ln_dw_amt_m1_7 ~ (", paste(theory_vars, collapse = " + "), ") * dw_m1_7_repeat + ", ctrl_str)
)

fml_interact_8_12 <- as.formula(
  paste0("ln_dw_amt_m8_12 ~ (", paste(theory_vars, collapse = " + "), ") * dw_m8_12_repeat + ", ctrl_str)
)

# 6. Run Regressions (Checking degrees of freedom first)
reg_interact <- list()

if(nrow(df_borr_1_7) > 15) {
  reg_interact[["Mar 1-7 (Interacted)"]] <- feols(fml_interact_1_7, data = df_borr_1_7, vcov = "hetero")
} else {
  cat("Not enough borrowers in Mar 1-7 to run the interaction.\n")
}

if(nrow(df_borr_8_12) > 15) {
  reg_interact[["Mar 8-12 (Interacted)"]] <- feols(fml_interact_8_12, data = df_borr_8_12, vcov = "hetero")
} else {
  cat("Not enough borrowers in Mar 8-12 to run the interaction.\n")
}

# 7. Print and Save
if(length(reg_interact) > 0) {
  
  save_reg_latex(reg_interact, "Table_DW_Early_Repeat_Interacted",
    title = "Intensive Margin: Repeat Borrower Sensitivities (Theory Interacted, Controls Standalone)")
    
  interact_table <- msummary(reg_interact, stars = c("*"=.10, "**"=.05, "***"=.01), gof_omit = "AIC|BIC|Log|RMSE",
    title = "Intensive Margin: Repeat Borrower Sensitivities (Theory Interacted, Controls Standalone)")
    
  interact_table
}
## Saved: Table_DW_Early_Repeat_Interacted.tex
Intensive Margin: Repeat Borrower Sensitivities (Theory Interacted, Controls Standalone)
Mar 1-7 (Interacted) Mar 8-12 (Interacted)
* p < 0.1, ** p < 0.05, *** p < 0.01
(Intercept) 12.362*** 12.579***
(0.562) (0.807)
mtm_total 0.597 -0.384
(0.537) (0.724)
uninsured_franchise 1.440*** -0.328
(0.353) (0.759)
mtm_x_f_u 0.590* 4.714***
(0.328) (0.830)
dw_m1_7_repeat 3.764***
(0.694)
ln_assets 0.338 0.964**
(0.534) (0.448)
cash_ratio -0.843* -0.446
(0.451) (0.673)
loan_to_deposit 0.615 -0.242
(0.451) (0.572)
book_equity_ratio -0.114 0.250
(0.400) (0.807)
wholesale 0.651*** -0.008
(0.174) (0.343)
roa -0.257 -1.018
(0.566) (0.644)
mtm_total × dw_m1_7_repeat -0.659
(0.619)
uninsured_franchise × dw_m1_7_repeat -1.486***
(0.434)
mtm_x_f_u × dw_m1_7_repeat 0.293
(0.509)
dw_m8_12_repeat 3.653***
(0.840)
mtm_total × dw_m8_12_repeat 0.278
(0.837)
uninsured_franchise × dw_m8_12_repeat 0.394
(0.785)
mtm_x_f_u × dw_m8_12_repeat -3.997***
(0.974)
Num.Obs. 67 66
R2 0.640 0.551
R2 Adj. 0.552 0.439
Std.Errors Heteroskedasticity-robust Heteroskedasticity-robust

16 Comparing aggregate, BTFP, and DW borrowing probabilities across two distinct periods: the heavy initial panic (Mar 13 - Apr 28) and the subsequent tail (Apr 29 - May 4).

# Range 1: March 13 - April 28
reg_full_r1 <- list(
  "Any Fed (Mar 13-Apr 28)" = feols(as.formula(paste("any_fed_r1 ~", fml_rhs_base)), data = df_model_extended, vcov = "hetero"),
  "BTFP (Mar 13-Apr 28)"    = feols(as.formula(paste("btfp_r1 ~", fml_rhs_base)), data = df_model_extended, vcov = "hetero"),
  "DW (Mar 13-Apr 28)"      = feols(as.formula(paste("dw_r1 ~", fml_rhs_base)), data = df_model_extended, vcov = "hetero")
)

# Range 2: April 29 - May 4
reg_full_r2 <- list(
  "Any Fed (Apr 29-May 4)" = feols(as.formula(paste("any_fed_r2 ~", fml_rhs_base)), data = df_model_extended, vcov = "hetero"),
  "BTFP (Apr 29-May 4)"    = feols(as.formula(paste("btfp_r2 ~", fml_rhs_base)), data = df_model_extended, vcov = "hetero"),
  "DW (Apr 29-May 4)"      = feols(as.formula(paste("dw_r2 ~", fml_rhs_base)), data = df_model_extended, vcov = "hetero")
)

msummary(reg_full_r1, stars = c("*"=.10, "**"=.05, "***"=.01), gof_omit = "AIC|BIC|Log|RMSE",
  add_rows = dv_summary_rows(reg_full_r1), coef_rename = coef_map,
  title = "New Full Regression 1: March 13 - April 28")
New Full Regression 1: March 13 - April 28
Any Fed (Mar 13-Apr 28) BTFP (Mar 13-Apr 28) DW (Mar 13-Apr 28)
* p < 0.1, ** p < 0.05, *** p < 0.01
(Intercept) 0.176*** 0.106*** 0.091***
(0.006) (0.005) (0.004)
ℓ (MTM Loss) 0.028*** 0.018*** 0.012**
(0.006) (0.005) (0.005)
f^U (Uninsured Franchise) 0.019*** 0.015** 0.013**
(0.007) (0.006) (0.006)
ℓ × f^U (KEY: c₃>0) 0.021*** 0.014*** 0.016***
(0.006) (0.005) (0.005)
Log(Assets) 0.096*** 0.052*** 0.067***
(0.007) (0.006) (0.006)
Cash/TA -0.016*** -0.021*** 0.002
(0.006) (0.004) (0.005)
Loan/Dep -0.003 -0.006 0.001
(0.007) (0.005) (0.005)
Equity/TA -0.018*** -0.016*** -0.002
(0.005) (0.004) (0.004)
Wholesale 0.026*** 0.019*** 0.014***
(0.007) (0.006) (0.005)
ROA -0.001 -0.003 -0.001
(0.006) (0.004) (0.004)
Num.Obs. 4226 4226 4226
R2 0.104 0.064 0.068
R2 Adj. 0.102 0.062 0.066
Std.Errors Heteroskedasticity-robust Heteroskedasticity-robust Heteroskedasticity-robust
N (DV = 1) 742 448 382
Mean (DV) 0.1756 0.106 0.0904 
save_reg_latex(reg_full_r1, "Table_Full_Reg_R1",
  add_rows = dv_summary_rows(reg_full_r1), coef_rename = coef_map,
  title = "Emergency Facility Borrowing: March 13 - April 28")
## Saved: Table_Full_Reg_R1.tex
msummary(reg_full_r2, stars = c("*"=.10, "**"=.05, "***"=.01), gof_omit = "AIC|BIC|Log|RMSE",
  add_rows = dv_summary_rows(reg_full_r2), coef_rename = coef_map,
  title = "New Full Regression 2: April 29 - May 4")
New Full Regression 2: April 29 - May 4
Any Fed (Apr 29-May 4) BTFP (Apr 29-May 4) DW (Apr 29-May 4)
* p < 0.1, ** p < 0.05, *** p < 0.01
(Intercept) 0.058*** 0.039*** 0.021***
(0.004) (0.003) (0.002)
ℓ (MTM Loss) 0.007* 0.007** 0.001
(0.004) (0.003) (0.002)
f^U (Uninsured Franchise) 0.008* 0.009** 0.000
(0.005) (0.004) (0.003)
ℓ × f^U (KEY: c₃>0) 0.008** 0.004 0.006**
(0.004) (0.003) (0.003)
Log(Assets) 0.026*** 0.016*** 0.012***
(0.005) (0.004) (0.003)
Cash/TA -0.015*** -0.013*** -0.002
(0.003) (0.003) (0.002)
Loan/Dep -0.005 -0.007** 0.001
(0.004) (0.003) (0.003)
Equity/TA -0.004 -0.005 -0.000
(0.004) (0.003) (0.002)
Wholesale 0.016*** 0.011*** 0.008***
(0.005) (0.004) (0.003)
ROA -0.007** -0.006** -0.002
(0.003) (0.003) (0.002)
Num.Obs. 4226 4226 4226
R2 0.033 0.027 0.014
R2 Adj. 0.031 0.025 0.012
Std.Errors Heteroskedasticity-robust Heteroskedasticity-robust Heteroskedasticity-robust
N (DV = 1) 243 167 86
Mean (DV) 0.0575 0.0395 0.0204 
save_reg_latex(reg_full_r2, "Table_Full_Reg_R2",
  add_rows = dv_summary_rows(reg_full_r2), coef_rename = coef_map,
  title = "Emergency Facility Borrowing: April 29 - May 4")
## Saved: Table_Full_Reg_R2.tex

17 Subsamples by Size: Range 1 (Mar 13 - Apr 28)

run_size_regs <- function(dv, base_df) {
  fml <- as.formula(paste(dv, "~", fml_rhs_base))
  list(
    "< 1B"     = feols(fml, data = base_df %>% filter(size_bucket == "< 1B"), vcov = "hetero"),
    "1B - 10B" = feols(fml, data = base_df %>% filter(size_bucket == "1B - 10B"), vcov = "hetero"),
    "> 10B"    = feols(fml, data = base_df %>% filter(size_bucket == "> 10B"), vcov = "hetero")
  )
}

reg_r1_any  <- run_size_regs("any_fed_r1", df_model_extended)
reg_r1_btfp <- run_size_regs("btfp_r1", df_model_extended)
reg_r1_dw   <- run_size_regs("dw_r1", df_model_extended)

msummary(reg_r1_any, stars = c("*"=.10, "**"=.05, "***"=.01), gof_omit = "AIC|BIC|Log|RMSE",
  add_rows = dv_summary_rows(reg_r1_any), coef_rename = coef_map,
  title = "Size Subsamples (Mar 13 - Apr 28): ANY FED")
Size Subsamples (Mar 13 - Apr 28): ANY FED
< 1B 1B - 10B > 10B
* p < 0.1, ** p < 0.05, *** p < 0.01
(Intercept) 0.168*** 0.233*** 4.415**
(0.008) (0.050) (2.044)
ℓ (MTM Loss) 0.022*** 0.042** -0.029
(0.007) (0.020) (0.060)
f^U (Uninsured Franchise) 0.012 0.009 0.117***
(0.008) (0.015) (0.032)
ℓ × f^U (KEY: c₃>0) 0.011* 0.033** 0.042
(0.006) (0.013) (0.035)
Log(Assets) 0.090*** 0.058 -1.579**
(0.010) (0.037) (0.792)
Cash/TA -0.011* -0.047* -0.116*
(0.006) (0.028) (0.068)
Loan/Dep -0.001 -0.026 0.027
(0.007) (0.023) (0.066)
Equity/TA -0.018*** -0.032 -0.107*
(0.005) (0.025) (0.061)
Wholesale 0.027*** 0.022 -0.014
(0.007) (0.017) (0.036)
ROA 0.000 0.019 -0.005
(0.005) (0.019) (0.057)
Num.Obs. 3309 783 134
R2 0.060 0.047 0.151
R2 Adj. 0.057 0.036 0.089
Std.Errors Heteroskedasticity-robust Heteroskedasticity-robust Heteroskedasticity-robust
N (DV = 1) 429 250 63
Mean (DV) 0.1296 0.3193 0.4701 
msummary(reg_r1_btfp, stars = c("*"=.10, "**"=.05, "***"=.01), gof_omit = "AIC|BIC|Log|RMSE",
  add_rows = dv_summary_rows(reg_r1_btfp), coef_rename = coef_map,
  title = "Size Subsamples (Mar 13 - Apr 28): BTFP")
Size Subsamples (Mar 13 - Apr 28): BTFP
< 1B 1B - 10B > 10B
* p < 0.1, ** p < 0.05, *** p < 0.01
(Intercept) 0.104*** 0.141*** 0.790
(0.006) (0.040) (2.142)
ℓ (MTM Loss) 0.016*** 0.028* -0.046
(0.005) (0.016) (0.054)
f^U (Uninsured Franchise) 0.014* 0.000 0.076**
(0.007) (0.013) (0.034)
ℓ × f^U (KEY: c₃>0) 0.012** 0.010 0.023
(0.005) (0.013) (0.033)
Log(Assets) 0.048*** 0.025 -0.229
(0.008) (0.030) (0.830)
Cash/TA -0.018*** -0.036* -0.067
(0.004) (0.021) (0.061)
Loan/Dep -0.005 -0.017 0.024
(0.005) (0.019) (0.062)
Equity/TA -0.012*** -0.053*** -0.094*
(0.004) (0.020) (0.055)
Wholesale 0.015*** 0.027* 0.010
(0.006) (0.016) (0.035)
ROA 0.000 -0.010 0.015
(0.004) (0.015) (0.051)
Num.Obs. 3309 783 134
R2 0.042 0.051 0.067
R2 Adj. 0.039 0.040 -0.001
Std.Errors Heteroskedasticity-robust Heteroskedasticity-robust Heteroskedasticity-robust
N (DV = 1) 268 140 40
Mean (DV) 0.081 0.1788 0.2985 
msummary(reg_r1_dw, stars = c("*"=.10, "**"=.05, "***"=.01), gof_omit = "AIC|BIC|Log|RMSE",
  add_rows = dv_summary_rows(reg_r1_dw), coef_rename = coef_map,
  title = "Size Subsamples (Mar 13 - Apr 28): DW")
Size Subsamples (Mar 13 - Apr 28): DW
< 1B 1B - 10B > 10B
* p < 0.1, ** p < 0.05, *** p < 0.01
(Intercept) 0.079*** 0.089** 5.199**
(0.006) (0.042) (1.994)
ℓ (MTM Loss) 0.008* 0.013 0.011
(0.005) (0.017) (0.046)
f^U (Uninsured Franchise) 0.004 0.012 0.085**
(0.006) (0.013) (0.035)
ℓ × f^U (KEY: c₃>0) 0.005 0.028** 0.061**
(0.005) (0.012) (0.029)
Log(Assets) 0.049*** 0.069** -1.930**
(0.008) (0.032) (0.772)
Cash/TA 0.005 -0.016 -0.082
(0.005) (0.025) (0.062)
Loan/Dep 0.002 -0.009 0.022
(0.005) (0.020) (0.044)
Equity/TA -0.007** 0.002 -0.046
(0.004) (0.021) (0.055)
Wholesale 0.013** 0.006 0.022
(0.005) (0.014) (0.035)
ROA 0.000 0.015 0.003
(0.004) (0.017) (0.050)
Num.Obs. 3309 783 134
R2 0.026 0.020 0.144
R2 Adj. 0.023 0.009 0.082
Std.Errors Heteroskedasticity-robust Heteroskedasticity-robust Heteroskedasticity-robust
N (DV = 1) 194 145 43
Mean (DV) 0.0586 0.1852 0.3209 
# Save combined latex (or separate as preferred)
save_reg_latex(reg_r1_any, "Table_Size_R1_AnyFed", add_rows = dv_summary_rows(reg_r1_any), coef_rename = coef_map, title = "Any Fed by Size (Mar 1-Apr 28)")
## Saved: Table_Size_R1_AnyFed.tex
save_reg_latex(reg_r1_btfp, "Table_Size_R1_BTFP", add_rows = dv_summary_rows(reg_r1_btfp), coef_rename = coef_map, title = "BTFP by Size (Mar 1-Apr 28)")
## Saved: Table_Size_R1_BTFP.tex
save_reg_latex(reg_r1_dw, "Table_Size_R1_DW", add_rows = dv_summary_rows(reg_r1_dw), coef_rename = coef_map, title = "DW by Size (Mar 1-Apr 28)")
## Saved: Table_Size_R1_DW.tex

18 Subsamples by Size: Range 2 (Apr 29 - May 4)

reg_r2_any  <- run_size_regs("any_fed_r2", df_model_extended)
reg_r2_btfp <- run_size_regs("btfp_r2", df_model_extended)
reg_r2_dw   <- run_size_regs("dw_r2", df_model_extended)

msummary(reg_r2_any, stars = c("*"=.10, "**"=.05, "***"=.01), gof_omit = "AIC|BIC|Log|RMSE",
  add_rows = dv_summary_rows(reg_r2_any), coef_rename = coef_map,
  title = "Size Subsamples (Apr 29 - May 4): ANY FED")
Size Subsamples (Apr 29 - May 4): ANY FED
< 1B 1B - 10B > 10B
* p < 0.1, ** p < 0.05, *** p < 0.01
(Intercept) 0.056*** 0.097*** 0.702
(0.005) (0.037) (1.528)
ℓ (MTM Loss) -0.006 0.037*** 0.019
(0.004) (0.014) (0.035)
f^U (Uninsured Franchise) 0.005 0.001 0.035
(0.005) (0.011) (0.024)
ℓ × f^U (KEY: c₃>0) -0.003 0.012 0.054**
(0.004) (0.010) (0.022)
Log(Assets) 0.030*** 0.008 -0.247
(0.006) (0.026) (0.594)
Cash/TA -0.014*** -0.028* -0.060*
(0.003) (0.017) (0.032)
Loan/Dep -0.004 -0.023 -0.020
(0.004) (0.015) (0.030)
Equity/TA -0.006* 0.003 -0.000
(0.004) (0.017) (0.037)
Wholesale 0.011** 0.034** 0.003
(0.005) (0.014) (0.021)
ROA -0.006* -0.010 0.012
(0.003) (0.011) (0.021)
Num.Obs. 3309 783 134
R2 0.021 0.055 0.095
R2 Adj. 0.019 0.044 0.029
Std.Errors Heteroskedasticity-robust Heteroskedasticity-robust Heteroskedasticity-robust
N (DV = 1) 145 84 14
Mean (DV) 0.0438 0.1073 0.1045 
msummary(reg_r2_btfp, stars = c("*"=.10, "**"=.05, "***"=.01), gof_omit = "AIC|BIC|Log|RMSE",
  add_rows = dv_summary_rows(reg_r2_btfp), coef_rename = coef_map,
  title = "Size Subsamples (Apr 29 - May 4): BTFP")
Size Subsamples (Apr 29 - May 4): BTFP
< 1B 1B - 10B > 10B
* p < 0.1, ** p < 0.05, *** p < 0.01
(Intercept) 0.040*** 0.082*** 0.355
(0.004) (0.031) (1.240)
ℓ (MTM Loss) -0.003 0.028** 0.025
(0.003) (0.012) (0.023)
f^U (Uninsured Franchise) 0.005 0.005 0.042*
(0.004) (0.009) (0.022)
ℓ × f^U (KEY: c₃>0) -0.003 0.004 0.039**
(0.003) (0.009) (0.019)
Log(Assets) 0.022*** -0.009 -0.122
(0.005) (0.022) (0.482)
Cash/TA -0.012*** -0.025** -0.023
(0.003) (0.011) (0.019)
Loan/Dep -0.006* -0.022* -0.022
(0.003) (0.013) (0.027)
Equity/TA -0.005 -0.005 -0.027
(0.003) (0.013) (0.029)
Wholesale 0.008** 0.018 0.009
(0.004) (0.011) (0.019)
ROA -0.004* -0.016* 0.018
(0.003) (0.009) (0.015)
Num.Obs. 3309 783 134
R2 0.019 0.049 0.118
R2 Adj. 0.016 0.038 0.054
Std.Errors Heteroskedasticity-robust Heteroskedasticity-robust Heteroskedasticity-robust
N (DV = 1) 103 55 9
Mean (DV) 0.0311 0.0702 0.0672 
msummary(reg_r2_dw, stars = c("*"=.10, "**"=.05, "***"=.01), gof_omit = "AIC|BIC|Log|RMSE",
  add_rows = dv_summary_rows(reg_r2_dw), coef_rename = coef_map,
  title = "Size Subsamples (Apr 29 - May 4): DW")
Size Subsamples (Apr 29 - May 4): DW
< 1B 1B - 10B > 10B
* p < 0.1, ** p < 0.05, *** p < 0.01
(Intercept) 0.019*** 0.021 0.556
(0.003) (0.025) (0.948)
ℓ (MTM Loss) -0.004* 0.010 0.002
(0.002) (0.009) (0.029)
f^U (Uninsured Franchise) -0.000 -0.001 -0.005
(0.003) (0.008) (0.012)
ℓ × f^U (KEY: c₃>0) 0.000 0.014* 0.022
(0.002) (0.008) (0.015)
Log(Assets) 0.010** 0.015 -0.205
(0.004) (0.018) (0.368)
Cash/TA -0.002 -0.001 -0.045
(0.002) (0.013) (0.029)
Loan/Dep 0.002 -0.004 -0.007
(0.003) (0.008) (0.017)
Equity/TA -0.002 0.002 0.026
(0.002) (0.012) (0.024)
Wholesale 0.005* 0.020** 0.005
(0.003) (0.010) (0.015)
ROA -0.002 0.003 0.006
(0.002) (0.008) (0.019)
Num.Obs. 3309 783 134
R2 0.007 0.026 0.053
R2 Adj. 0.004 0.015 -0.016
Std.Errors Heteroskedasticity-robust Heteroskedasticity-robust Heteroskedasticity-robust
N (DV = 1) 48 32 6
Mean (DV) 0.0145 0.0409 0.0448 
# Save outputs
save_reg_latex(reg_r2_any, "Table_Size_R2_AnyFed", add_rows = dv_summary_rows(reg_r2_any), coef_rename = coef_map, title = "Any Fed by Size (Apr 29-May 4)")
## Saved: Table_Size_R2_AnyFed.tex
save_reg_latex(reg_r2_btfp, "Table_Size_R2_BTFP", add_rows = dv_summary_rows(reg_r2_btfp), coef_rename = coef_map, title = "BTFP by Size (Apr 29-May 4)")
## Saved: Table_Size_R2_BTFP.tex
save_reg_latex(reg_r2_dw, "Table_Size_R2_DW", add_rows = dv_summary_rows(reg_r2_dw), coef_rename = coef_map, title = "DW by Size (Apr 29-May 4)")
## Saved: Table_Size_R2_DW.tex

18.1 Daily Borrowers: March 1 - May 4

# 1. Aggregate daily unique borrowers for BTFP
btfp_daily <- btfp_loans %>%
  filter(btfp_loan_date >= as.Date("2023-03-01") & btfp_loan_date <= as.Date("2023-05-04")) %>%
  group_by(btfp_loan_date) %>%
  summarise(Borrowers = n_distinct(rssd_id), .groups = "drop") %>%
  rename(Date = btfp_loan_date) %>%
  mutate(Facility = "BTFP")

# 2. Aggregate daily unique borrowers for DW
dw_daily <- dw_loans %>%
  filter(dw_loan_date >= as.Date("2023-03-01") & dw_loan_date <= as.Date("2023-05-04")) %>%
  group_by(dw_loan_date) %>%
  summarise(Borrowers = n_distinct(rssd_id), .groups = "drop") %>%
  rename(Date = dw_loan_date) %>%
  mutate(Facility = "Discount Window")

# 3. Combine into a single dataframe
daily_counts <- bind_rows(btfp_daily, dw_daily)

# 4. Create the plot
# Using geom_col (bar plot) works best to show the episodic spikes in DW vs steady BTFP usage
p_daily <- ggplot(daily_counts, aes(x = Date, y = Borrowers, fill = Facility)) +
  geom_col(position = position_dodge(width = 0.8), width = 0.8, alpha = 0.85) +
  # Add vertical lines to delineate key date ranges used in regressions
  geom_vline(xintercept = as.numeric(as.Date("2023-03-08")), linetype = "dashed", color = "grey40", linewidth = 0.6) +
  geom_vline(xintercept = as.numeric(as.Date("2023-04-28")), linetype = "dashed", color = "grey40", linewidth = 0.6) +
  scale_fill_manual(values = c("BTFP" = "#1565C0", "Discount Window" = "#C62828")) +
  scale_x_date(date_breaks = "1 week", date_labels = "%b %d", expand = c(0.02, 0)) +
  labs(
    title = "Daily Number of Borrowers: BTFP vs. Discount Window",
    subtitle = "March 1, 2023 – May 4, 2023. Dashed lines denote Mar 8 (Initial panic) and Apr 28.",
    x = "Date",
    y = "Number of Unique Borrowers"
  ) +
  theme_gp +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1),
    legend.title = element_blank()
  )

print(p_daily)

# Optionally save the figure
save_figure(p_daily, "Plot_Daily_Borrowers_Mar_May", width = 12, height = 6)
# 1. Aggregate daily unique borrowers for BTFP
btfp_daily <- btfp_loans %>%
  filter(btfp_loan_date >= as.Date("2023-03-01") & btfp_loan_date <= as.Date("2023-05-04")) %>%
  group_by(btfp_loan_date) %>%
  summarise(Borrowers = n_distinct(rssd_id), .groups = "drop") %>%
  rename(Date = btfp_loan_date) %>%
  mutate(Facility = "BTFP")

# 2. Aggregate daily unique borrowers for DW
dw_daily <- dw_loans %>%
  filter(dw_loan_date >= as.Date("2023-03-01") & dw_loan_date <= as.Date("2023-05-04")) %>%
  group_by(dw_loan_date) %>%
  summarise(Borrowers = n_distinct(rssd_id), .groups = "drop") %>%
  rename(Date = dw_loan_date) %>%
  mutate(Facility = "Discount Window")

# 3. Combine into a single dataframe
daily_counts <- bind_rows(btfp_daily, dw_daily)

# 4. Create the plot
p_daily <- ggplot(daily_counts, aes(x = Date, y = Borrowers, fill = Facility)) +
  geom_col(position = position_dodge(width = 0.8), width = 0.8, alpha = 0.85) +
  
  # --- NEW: Add the actual numbers on top of the bars ---
  # ifelse hides '0' values to prevent clutter on days with no borrowing
  geom_text(aes(label = ifelse(Borrowers > 0, Borrowers, "")), 
            position = position_dodge(width = 0.8), 
            vjust = -0.5,   # Nudges the text just above the bar
            size = 3.5,     # Adjust size if it looks too crowded
            color = "black") +
            
  # Add vertical lines to delineate key date ranges
  geom_vline(xintercept = as.numeric(as.Date("2023-03-08")), linetype = "dashed", color = "grey40", linewidth = 0.6) +
  geom_vline(xintercept = as.numeric(as.Date("2023-04-28")), linetype = "dashed", color = "grey40", linewidth = 0.6) +
  scale_fill_manual(values = c("BTFP" = "#1565C0", "Discount Window" = "#C62828")) +
  scale_x_date(date_breaks = "1 week", date_labels = "%b %d", expand = c(0.02, 0)) +
  
  # --- NEW: Expand the top of the y-axis so high numbers don't get cut off ---
  scale_y_continuous(expand = expansion(mult = c(0, 0.15))) +
  
  labs(
    title = "Daily Number of Borrowers: BTFP vs. Discount Window",
    subtitle = "March 1, 2023 – May 4, 2023. Actual borrower counts displayed on bars.",
    x = "Date",
    y = "Number of Unique Borrowers"
  ) +
  theme_gp +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1),
    legend.title = element_blank()
  )

print(p_daily)
## Warning in scale_x_date(date_breaks = "1 week", date_labels = "%b %d", expand = c(0.02, : A <numeric> value was passed to a Date scale.
## ℹ The value was converted to a <Date> object.
## A <numeric> value was passed to a Date scale.
## ℹ The value was converted to a <Date> object.

# Optionally save the figure
save_figure(p_daily, "Plot_Daily_Borrowers_Mar_May_Labeled", width = 18, height = 8)
## Warning in scale_x_date(date_breaks = "1 week", date_labels = "%b %d", expand = c(0.02, : A <numeric> value was passed to a Date scale.
## ℹ The value was converted to a <Date> object.
## A <numeric> value was passed to a Date scale.
## ℹ The value was converted to a <Date> object.
## Saved: Plot_Daily_Borrowers_Mar_May_Labeled.pdf