Summary

This analysis examines who used the Bank Term Funding Program (BTFP) during the 2023 banking crisis and what determined that choice. Building on Jiang et al. (2024) Market-to-market (MTM) Loss, we test whether BTFP’s unique par valuation feature attracted banks facing insolvency risk (not just liquidity risk).

Key Hypotheses:

Banks with higher insolvency risk AND more OMO-eligible assets preferentially chose BTFP
BTFP usage is driven more by insolvency concerns than pure liquidity needs
The interaction of run risk and OMO-eligible capacity predicts BTFP usage (Borrowing)

Conceptual Framework

BTFP Design and Economic Incentives

The Bank Term Funding Program (BTFP) was announced on March 12, 2023, with a critical design feature:

Collateral valued at PAR rather than market value

This created differential incentives across banks:

Bank Type	BTFP Benefit	Expected Behavior
High MTM losses, High OMO-eligible	Large valuation subsidy	Strong BTFP preference
High MTM losses, Low OMO-eligible	Cannot fully benefit	May use DW or FHLB
Low MTM losses	Little subsidy value	Indifferent or prefer FHLB

Separating Liquidity Risk from Insolvency Risk

This is a key conceptual contribution. We distinguish:

Liquidity Risk (Run Risk)

\[\text{Run Risk} = \frac{\text{Uninsured Deposits}}{\text{Total Deposits}} \times \frac{\text{MTM Loss}}{\text{Total Assets}}\]

Interpretation: Banks with high uninsured deposits AND large MTM losses face the greatest liquidity pressure from potential runs.

Insolvency Risk (Insured Deposit Coverage Ratio)

\[\text{Insured Deposit Coverage Ratio} = \frac{\text{MV Assets} - s \times \text{Uninsured Deposits} - \text{Insured Deposits}}{\text{Insured Deposits}}\]

Where \(s \in [0,1]\) is the assumed run rate on uninsured deposits.

Interpretation: After paying off running depositors, can remaining assets cover FDIC-insured deposits? Negative values indicate insolvency.

Causal Identification Strategy

Identification Challenge

Banks that chose BTFP may differ from non-users in unobservable ways: - Management quality - Risk appetite - Regulatory relationships

Our Strategy

Pre-treatment conditioning: Use 2022Q4 data (before crisis)
Sample restriction: Focus on banks with omo_eligible > 0
Rich controls: Size, wholesale funding, FHLB usage
Interaction tests: Does subsidy value matter MORE for distressed banks?
Progressive model building: Start simple, add complexity

Key Identifying Assumption

Conditional on observable 2022Q4 characteristics, BTFP usage is driven by economic fundamentals (insolvency risk, collateral availability), not unobservable factors.

Variable Definitions

Dependent Variables

Variable	Definition	Construction
`btfp_user`	Used BTFP during acute crisis	=1 if any BTFP loan Mar 13 - Apr 30, 2023
`dw_user`	Used Discount Window	=1 if any DW loan during crisis
`btfp_amount_z`	BTFP borrowing intensity	log(1 + BTFP amount / Assets), standardized

Key Independent Variables

Variable	Definition	Formula
`pct_uninsured_z`	Uninsured deposit share	Uninsured / Total Deposits, standardized
`pct_mtm_loss_z`	MTM loss ratio	MTM Loss / Total Assets, standardized
`run_risk_z`	Liquidity risk	pct_uninsured × pct_mtm_loss, standardized
`insured_coverage_z`	Insolvency risk	(MV Assets - s × Uninsured - Insured) / Insured, standardized
`valuation_subsidy_z`	BTFP benefit	MTM Loss on OMO-Eligible / Assets, standardized
`omo_eligible_z`	Eligible as Par-value collateral (treasury, agency_debt, agency_mbs)	OMO-Eligible / Assets, standardized

Control Variables

Variable	Definition
`log_assets`	log(Total Assets)
`fhlb_ratio_z`	FHLB Advances / Total Assets
`wholesale_ratio_z`	(Repo + Fed Funds + Other Borrowing) / Assets
`cash_ratio_z`	Cash / Uninsured Deposits
`book_equity_ratio_z`	Book Equity / Total Assets

Note: All continuous variables are standardized (z-scores) for interpretability.

Setup and Data Loading

################################################################################
# SECTION 0: SETUP AND CONFIGURATION
################################################################################

# Required packages
required_packages <- c(
  "data.table", "dplyr", "tidyr", "tibble", "lubridate", "stringr", "readr",
  "fixest", "sandwich", "lmtest", "broom",
  "ggplot2", "ggthemes", "scales", "patchwork",
  "knitr", "kableExtra",
  "DescTools", "moments"
)

# Install missing packages
missing_packages <- required_packages[!required_packages %in% installed.packages()[,"Package"]]
if (length(missing_packages) > 0) {
  cat("Installing:", paste(missing_packages, collapse = ", "), "\n")
  install.packages(missing_packages, repos = "https://cloud.r-project.org")
}

# Load packages
suppressPackageStartupMessages({
  library(data.table)
  library(dplyr)
  library(tidyr)
  library(tibble)
  library(lubridate)
  library(stringr)
  library(readr)
  library(fixest)
  library(sandwich)
  library(lmtest)
  library(broom)
  library(ggplot2)
  library(ggthemes)
  library(scales)
  library(patchwork)
  library(knitr)
  library(kableExtra)
  library(DescTools)
  library(moments)
})

options(scipen = 999)
set.seed(20230313)  # SVB failure date as seed

################################################################################
# Configure Paths (MODIFY FOR YOUR ENVIRONMENT)
################################################################################

BASE_PATH <- "C:/Users/mferdo2/OneDrive - Louisiana State University/Finance_PhD/DW_Stigma_paper/Liquidity_project_2025"

# Only set working directory if path exists (for portability)
if (dir.exists(BASE_PATH)) {
  setwd(BASE_PATH)
}

DATA_RAW    <- file.path(BASE_PATH, "01_data/raw")
DATA_PROC   <- file.path(BASE_PATH, "01_data/processed")
DOC_PATH    <- file.path(BASE_PATH, "03_documentation")
TABLE_PATH  <- file.path(DOC_PATH, "regression_tables/btfp_analysis")
FIG_PATH    <- file.path(DOC_PATH, "figures/btfp_analysis")

# Create output directories
dir.create(TABLE_PATH, recursive = TRUE, showWarnings = FALSE)
dir.create(FIG_PATH, recursive = TRUE, showWarnings = FALSE)

################################################################################
# Crisis Parameters
################################################################################

CRISIS_START      <- as.Date("2023-03-13")  # Monday after SVB weekend
CRISIS_END_ACUTE  <- as.Date("2023-04-30")  # End of acute phase
CRISIS_END_BTFP   <- as.Date("2024-03-11")  # BTFP program closure
MTM_START         <- as.Date("2022-03-16")  # First FFR hike
MTM_END           <- as.Date("2023-03-10")  # SVB failure

cat("\n")

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

#> ================================================================

cat("BTFP USAGE ANALYSIS - 2023 BANKING CRISIS\n")

#> BTFP USAGE ANALYSIS - 2023 BANKING CRISIS

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

#> ================================================================

cat("Crisis Window:     ", as.character(CRISIS_START), " to ", as.character(CRISIS_END_ACUTE), "\n")

#> Crisis Window:      2023-03-13  to  2023-04-30

cat("MTM Period:        ", as.character(MTM_START), " to ", as.character(MTM_END), "\n")

#> MTM Period:         2022-03-16  to  2023-03-10

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

#> ================================================================

################################################################################
# Helper Functions
################################################################################

#' Winsorize at specified quantiles
winsorize <- function(x, probs = c(0.01, 0.99)) {
  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-score transformation)
standardize <- function(x) {
  (x - mean(x, na.rm = TRUE)) / sd(x, na.rm = TRUE)
}

#' Safe division (returns NA if denominator is 0 or NA)
safe_div <- function(num, denom) {
  ifelse(is.na(denom) | denom == 0, NA_real_, num / denom)
}

#' Parse dates flexibly
parse_date_safe <- function(x) {
  if (inherits(x, "Date")) return(x)
  x_char <- as.character(x)
  result <- suppressWarnings(mdy(x_char))
  if (all(is.na(result))) result <- suppressWarnings(ymd(x_char))
  if (all(is.na(result))) result <- suppressWarnings(dmy(x_char))
  return(result)
}

#' Create summary statistics table
create_summary_table <- function(df, vars, digits = 3) {
  df %>%
    select(all_of(vars)) %>%
    pivot_longer(everything(), names_to = "Variable", values_to = "Value") %>%
    group_by(Variable) %>%
    summarise(
      N = sum(!is.na(Value)),
      Mean = round(mean(Value, na.rm = TRUE), digits),
      SD = round(sd(Value, na.rm = TRUE), digits),
      Min = round(min(Value, na.rm = TRUE), digits),
      P25 = round(quantile(Value, 0.25, na.rm = TRUE), digits),
      Median = round(quantile(Value, 0.50, na.rm = TRUE), digits),
      P75 = round(quantile(Value, 0.75, na.rm = TRUE), digits),
      Max = round(max(Value, na.rm = TRUE), digits),
      .groups = "drop"
    )
}

Data Loading

################################################################################
# SECTION 1: LOAD DATA
################################################################################

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

#> 
#> === LOADING DATA ===

# 1.1 Load Call Report Data (2022Q4 baseline)
call_q <- read_csv(
  file.path(DATA_PROC, "final_df_fhlb_gsib.csv"), 
  show_col_types = FALSE
) %>%
  mutate(idrssd = as.character(idrssd))

cat("✓ Call reports loaded:", nrow(call_q), "obs,", n_distinct(call_q$idrssd), "banks\n")

#> ✓ Call reports loaded: 61002 obs, 4926 banks

# 1.2 Load BTFP Loan Data
btfp_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 = parse_date_safe(btfp_loan_date)
  )

cat("✓ BTFP loans loaded:", nrow(btfp_raw), "obs\n")

#> ✓ BTFP loans loaded: 6734 obs

cat("  Date range:", as.character(min(btfp_raw$btfp_loan_date, na.rm = TRUE)), 
    "to", as.character(max(btfp_raw$btfp_loan_date, na.rm = TRUE)), "\n")

#>   Date range: 2023-03-13 to 2024-03-11

# 1.3 Load DW Loan Data
dw_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 = parse_date_safe(dw_loan_date)
  )

cat("✓ DW loans loaded:", nrow(dw_raw), "obs\n")

#> ✓ DW loans loaded: 7724 obs

cat("  Date range:", as.character(min(dw_raw$dw_loan_date, na.rm = TRUE)), 
    "to", as.character(max(dw_raw$dw_loan_date, na.rm = TRUE)), "\n")

#>   Date range: 2023-01-03 to 2023-09-29

Data Preparation

Aggregate Facility Usage

################################################################################
# SECTION 2: AGGREGATE FACILITY USAGE
################################################################################

cat("\n=== AGGREGATING FACILITY USAGE ===\n\n")

#> 
#> === AGGREGATING FACILITY USAGE ===

# 2.1 BTFP usage during acute crisis
btfp_crisis <- btfp_raw %>%
  filter(btfp_loan_date >= CRISIS_START, btfp_loan_date <= CRISIS_END_ACUTE) %>%
  group_by(rssd_id) %>%
  summarise(
    btfp_user = 1L,
    btfp_total_amt = sum(btfp_loan_amount, na.rm = TRUE),
    btfp_n_loans = n(),
    btfp_first_date = min(btfp_loan_date),
    .groups = "drop"
  ) %>%
  rename(idrssd = rssd_id)

cat("✓ BTFP users (acute crisis):", nrow(btfp_crisis), "\n")

#> ✓ BTFP users (acute crisis): 481

# 2.2 DW usage during acute crisis
dw_crisis <- dw_raw %>%
  filter(dw_loan_date >= CRISIS_START, dw_loan_date <= CRISIS_END_ACUTE) %>%
  group_by(rssd_id) %>%
  summarise(
    dw_user = 1L,
    dw_total_amt = sum(dw_loan_amount, na.rm = TRUE),
    dw_n_loans = n(),
    dw_first_date = min(dw_loan_date),
    .groups = "drop"
  ) %>%
  rename(idrssd = rssd_id)

cat("✓ DW users (acute crisis):", nrow(dw_crisis), "\n")

#> ✓ DW users (acute crisis): 415

# 2.3 Extended BTFP usage (full program)
btfp_full <- btfp_raw %>%
  filter(btfp_loan_date >= CRISIS_START, btfp_loan_date <= CRISIS_END_BTFP) %>%
  group_by(rssd_id) %>%
  summarise(
    btfp_user_full = 1L,
    btfp_total_amt_full = sum(btfp_loan_amount, na.rm = TRUE),
    .groups = "drop"
  ) %>%
  rename(idrssd = rssd_id)

cat("✓ BTFP users (full program):", nrow(btfp_full), "\n")

#> ✓ BTFP users (full program): 1327

Create 2022Q4 Baseline Variables

################################################################################
# SECTION 3: CREATE 2022Q4 BASELINE WITH KEY VARIABLES
################################################################################

cat("\n=== CREATING 2022Q4 BASELINE ===\n\n")

#> 
#> === CREATING 2022Q4 BASELINE ===

# Check available columns
has_failed_bank <- "failed_bank" %in% names(call_q)
has_mv_asset <- "mv_asset" %in% names(call_q)

cat("Data availability checks:\n")

#> Data availability checks:

cat("  - failed_bank column:", ifelse(has_failed_bank, "YES", "NO"), "\n")

#>   - failed_bank column: YES

cat("  - mv_asset column:", ifelse(has_mv_asset, "YES", "NO"), "\n\n")

#>   - mv_asset column: NO

# Create baseline dataset
baseline <- call_q %>%
  filter(quarter == "2022Q4") %>%
  transmute(
    idrssd = as.character(idrssd),
    size_bin = size_bin,
    
    # Failed bank flag
    failed_bank = if (has_failed_bank) failed_bank else 0L,
    
    # =========================================================================
    # A. BALANCE SHEET FUNDAMENTALS
    # =========================================================================
    total_asset,
    total_liability,
    total_equity,
    total_deposit,
    insured_deposit = pmax(insured_deposit, 1),
    uninsured_deposit = pmax(uninsured_deposit, 0),
    
    # Size measures
    log_assets = log(pmax(total_asset, 1)),
    
    # =========================================================================
    # B. MARK-TO-MARKET VARIABLES
    # =========================================================================
    mtm_loss_to_total_asset,
    
    # MV Equity (ALLOW NEGATIVE - captures truly insolvent banks)
    mv_equity = equity_after_mtm,
    
    # Market-value assets
    mv_asset = if (has_mv_asset) {
      mm_asset
    } else {
      total_asset * (1 - mtm_loss_to_total_asset)
    },
    
    # MTM insolvency indicator
    mtm_insolvent = (equity_after_mtm < 0),
    
    # =========================================================================
    # C. OMO-ELIGIBLE SECURITIES & VALUATION SUBSIDY
    # =========================================================================
    omo_eligible,
    mtm_loss_omo_eligible,
    mtm_loss_omo_eligible_to_total_asset,
    omo_eligible_to_total_asset,
    
    # =========================================================================
    # D. WHOLESALE FUNDING & FHLB
    # =========================================================================
    repo,
    fed_fund_purchase,
    other_borr,
    fhlb_adv,
    fhlb_less_than_1yr,
    
    # =========================================================================
    # E. LIQUIDITY & CONTROL VARIABLES
    # =========================================================================
    cash_to_uninsured_deposit,
    book_equity_to_total_asset
  )

cat("✓ Initial baseline:", nrow(baseline), "banks\n")

#> ✓ Initial baseline: 4737 banks

Construct Key Risk Variables

################################################################################
# SECTION 4: CONSTRUCT KEY RISK VARIABLES
################################################################################

cat("\n=== CONSTRUCTING KEY RISK VARIABLES ===\n\n")

#> 
#> === CONSTRUCTING KEY RISK VARIABLES ===

baseline <- baseline %>%
  mutate(
    # =========================================================================
    # KEY VARIABLE 1: PERCENT UNINSURED DEPOSITS
    # =========================================================================
    pct_uninsured = safe_div(uninsured_deposit, total_deposit)*100,
    
    # =========================================================================
    # KEY VARIABLE 2: PERCENT MTM LOSS
    # =========================================================================
    pct_mtm_loss = mtm_loss_to_total_asset,
    
    # =========================================================================
    # KEY VARIABLE 3: RUN RISK (Liquidity Risk)
    # Joint vulnerability: high uninsured AND high MTM loss
    # =========================================================================
    run_risk = pct_uninsured * pct_mtm_loss,
    
    # =========================================================================
    # KEY VARIABLE 4: INSURED DEPOSIT COVERAGE RATIO (Insolvency Risk)
    # After 100% uninsured run, can MV assets cover insured deposits?
    # Negative = insolvent
    # =========================================================================
    insured_coverage = safe_div(
      mv_asset - uninsured_deposit - insured_deposit,
      insured_deposit
    ),
    
    # Alternative: 50% uninsured run scenario
    insured_coverage_50pct = safe_div(
      mv_asset - 0.5 * uninsured_deposit - insured_deposit,
      insured_deposit
    ),
    
    # =========================================================================
    # KEY VARIABLE 5: VALUATION SUBSIDY (BTFP Benefit)
    # MTM loss on OMO-eligible = benefit from par valuation
    # =========================================================================
    valuation_subsidy = mtm_loss_omo_eligible_to_total_asset,
    
    # =========================================================================
    # KEY VARIABLE 6: OMO-ELIGIBLE RATIO
    # Capacity to use BTFP
    # =========================================================================
    omo_eligible_ratio = omo_eligible_to_total_asset,
    
    # Has any OMO-eligible collateral
    has_omo_eligible = (omo_eligible > 0),
    
    # =========================================================================
    # CONTROL VARIABLES
    # =========================================================================
    fhlb_ratio = safe_div(fhlb_adv, total_asset),
    fhlb_st_ratio = safe_div(fhlb_less_than_1yr, total_asset),
    wholesale_ratio = safe_div(repo + fed_fund_purchase + other_borr, total_asset),
    cash_ratio = cash_to_uninsured_deposit,
    book_equity_ratio = book_equity_to_total_asset
  )

# Filter valid observations
baseline <- baseline %>%
  filter(
    is.finite(pct_uninsured),
    is.finite(pct_mtm_loss),
    is.finite(run_risk),
    is.finite(insured_coverage)
  )

cat("✓ After constructing risk variables:", nrow(baseline), "banks\n")

#> ✓ After constructing risk variables: 4669 banks

# Summary of key variables
cat("\nKey Variable Summary:\n")

#> 
#> Key Variable Summary:

cat("  - Run Risk: Mean =", round(mean(baseline$run_risk, na.rm = TRUE), 4), 
    ", SD =", round(sd(baseline$run_risk, na.rm = TRUE), 4), "\n")

#>   - Run Risk: Mean = 148.6125 , SD = 294.2765

cat("  - Insured Coverage: Mean =", round(mean(baseline$insured_coverage, na.rm = TRUE), 4), 
    ", SD =", round(sd(baseline$insured_coverage, na.rm = TRUE), 4), "\n")

#>   - Insured Coverage: Mean = -20.727 , SD = 566.1702

cat("  - Valuation Subsidy: Mean =", round(mean(baseline$valuation_subsidy, na.rm = TRUE), 4), 
    ", SD =", round(sd(baseline$valuation_subsidy, na.rm = TRUE), 4), "\n")

#>   - Valuation Subsidy: Mean = 0.6314 , SD = 0.8385

cat("  - MTM Insolvent:", sum(baseline$mtm_insolvent, na.rm = TRUE), 
    "(", round(100 * mean(baseline$mtm_insolvent, na.rm = TRUE), 2), "%)\n")

#>   - MTM Insolvent: 866 ( 18.55 %)

cat("  - Has OMO-Eligible:", sum(baseline$has_omo_eligible, na.rm = TRUE), 
    "(", round(100 * mean(baseline$has_omo_eligible, na.rm = TRUE), 2), "%)\n")

#>   - Has OMO-Eligible: 4283 ( 91.73 %)

Winsorize and Standardize Variables

################################################################################
# SECTION 5: WINSORIZE AND STANDARDIZE VARIABLES
################################################################################

cat("\n=== WINSORIZING AND STANDARDIZING ===\n\n")

#> 
#> === WINSORIZING AND STANDARDIZING ===

# Variables to winsorize
vars_to_winsorize <- c(
  "pct_uninsured", "pct_mtm_loss", "run_risk", 
  "insured_coverage", "insured_coverage_50pct",
  "valuation_subsidy", "omo_eligible_ratio",
  "fhlb_ratio", "fhlb_st_ratio", "wholesale_ratio",
  "cash_ratio", "book_equity_ratio"
)

# Winsorize at 1st and 99th percentiles
baseline <- baseline %>%
  mutate(across(all_of(vars_to_winsorize), ~winsorize(.x, c(0.05, 0.95))))

cat("✓ Winsorized", length(vars_to_winsorize), "variables at 1/99 percentiles\n")

#> ✓ Winsorized 12 variables at 1/99 percentiles

# Create standardized (z-score) versions
baseline <- baseline %>%
  mutate(
    # Key independent variables (standardized)
    pct_uninsured_z = standardize(pct_uninsured),
    pct_mtm_loss_z = standardize(pct_mtm_loss),
    run_risk_z = standardize(run_risk),
    insured_coverage_z = standardize(insured_coverage),
    valuation_subsidy_z = standardize(valuation_subsidy),
    omo_eligible_z = standardize(omo_eligible_ratio),
    
    # Controls (standardized)
    fhlb_ratio_z = standardize(fhlb_ratio),
    fhlb_st_ratio_z = standardize(fhlb_st_ratio),
    wholesale_ratio_z = standardize(wholesale_ratio),
    cash_ratio_z = standardize(cash_ratio),
    book_equity_ratio_z = standardize(book_equity_ratio),
    log_assets_z = standardize(log_assets)
  )

cat("✓ Created standardized (z-score) versions of all variables\n")

#> ✓ Created standardized (z-score) versions of all variables

Merge Facility Usage and Create Analysis Dataset

################################################################################
# SECTION 6: MERGE AND CREATE FINAL ANALYSIS DATASET
################################################################################

cat("\n=== CREATING ANALYSIS DATASET ===\n\n")

#> 
#> === CREATING ANALYSIS DATASET ===

# Merge facility usage
analysis <- baseline %>%
  left_join(btfp_crisis, by = "idrssd") %>%
  left_join(dw_crisis, by = "idrssd") %>%
  left_join(btfp_full, by = "idrssd") %>%
  mutate(
    # Fill NA with 0 for non-users
    btfp_user = coalesce(btfp_user, 0L),
    dw_user = coalesce(dw_user, 0L),
    btfp_user_full = coalesce(btfp_user_full, 0L),
    btfp_total_amt = coalesce(btfp_total_amt, 0),
    dw_total_amt = coalesce(dw_total_amt, 0),
    
    # Scale amounts by assets
    btfp_amt_scaled = btfp_total_amt / pmax(total_asset * 1000, 1),
    dw_amt_scaled = dw_total_amt / pmax(total_asset * 1000, 1),
    
    # Log-transformed amounts (for intensive margin)
    btfp_amt_log = log1p(btfp_amt_scaled * 1e6),
    dw_amt_log = log1p(dw_amt_scaled * 1e6),
    
    # Standardized amounts
    btfp_amt_z = standardize(btfp_amt_scaled),
    
    # Usage patterns
    used_any_fed = pmax(btfp_user, dw_user),
    used_both = (btfp_user == 1 & dw_user == 1),
    used_only_btfp = (btfp_user == 1 & dw_user == 0),
    used_only_dw = (btfp_user == 0 & dw_user == 1),
    
    # For conditional choice model
    chose_btfp = ifelse(used_only_btfp | used_only_dw, 
                        as.integer(used_only_btfp), NA_integer_)
  )

# Exclude failed banks
n_before <- nrow(analysis)
failed_ids <- analysis %>% filter(failed_bank == 1) %>% pull(idrssd)

analysis_clean <- analysis %>%
  filter(failed_bank == 0)

cat("✓ Excluded", n_before - nrow(analysis_clean), "failed banks\n")

#> ✓ Excluded 8 failed banks

if (length(failed_ids) > 0) {
  cat("  Failed bank IDs:", paste(failed_ids, collapse = ", "), "\n")
}

#>   Failed bank IDs: 802866, 1216826, 2942690, 3076604, 3350724, 3390627, 3437483, 4114567

# Final sample
cat("\n✓ FINAL ANALYSIS SAMPLE:", nrow(analysis_clean), "banks\n")

#> 
#> ✓ FINAL ANALYSIS SAMPLE: 4661 banks

cat("  - BTFP users:", sum(analysis_clean$btfp_user), 
    "(", round(100*mean(analysis_clean$btfp_user), 2), "%)\n")

#>   - BTFP users: 473 ( 10.15 %)

cat("  - DW users:", sum(analysis_clean$dw_user), 
    "(", round(100*mean(analysis_clean$dw_user), 2), "%)\n")

#>   - DW users: 407 ( 8.73 %)

cat("  - Both facilities:", sum(analysis_clean$used_both), "\n")

#>   - Both facilities: 94

cat("  - Only BTFP:", sum(analysis_clean$used_only_btfp), "\n")

#>   - Only BTFP: 379

cat("  - Only DW:", sum(analysis_clean$used_only_dw), "\n")

#>   - Only DW: 313

cat("  - MTM Insolvent:", sum(analysis_clean$mtm_insolvent), 
    "(", round(100*mean(analysis_clean$mtm_insolvent), 2), "%)\n")

#>   - MTM Insolvent: 864 ( 18.54 %)

cat("  - Has OMO-Eligible:", sum(analysis_clean$has_omo_eligible), 
    "(", round(100*mean(analysis_clean$has_omo_eligible), 2), "%)\n")

#>   - Has OMO-Eligible: 4275 ( 91.72 %)

Summary Statistics

Full Sample Summary

################################################################################
# SECTION 7: SUMMARY STATISTICS
################################################################################

cat("\n=== SUMMARY STATISTICS ===\n\n")

#> 
#> === SUMMARY STATISTICS ===

# Variables for summary table
vars_summary <- c(
  "pct_uninsured", "pct_mtm_loss", "run_risk",
  "insured_coverage", "valuation_subsidy", "omo_eligible_ratio",
  "fhlb_ratio", "wholesale_ratio", "cash_ratio", "book_equity_ratio",
  "log_assets"
)

# Full sample summary
summary_full <- create_summary_table(analysis_clean, vars_summary)

# Display table
kable(summary_full, 
      caption = "Table 1: Summary Statistics - Full Sample",
      digits = 3) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
                full_width = FALSE)

Table 1: Summary Statistics - Full Sample
Variable	N	Mean	SD	Min	P25	Median	P75	Max
book_equity_ratio	4661	9.491	3.343	4.496	7.223	8.959	11.066	17.848
cash_ratio	4651	0.449	0.509	0.051	0.119	0.243	0.540	2.000
fhlb_ratio	4661	0.025	0.035	0.000	0.000	0.004	0.042	0.115
insured_coverage	4661	-8.683	3.454	-15.467	-11.150	-8.412	-6.000	-3.002
log_assets	4661	12.869	1.515	8.033	11.883	12.695	13.614	21.887
omo_eligible_ratio	4661	9.245	8.529	0.000	2.302	6.742	14.056	29.680
pct_mtm_loss	4661	5.433	2.060	1.977	3.818	5.294	6.966	9.233
pct_uninsured	4661	27.297	12.719	7.910	17.749	25.729	35.124	55.227
run_risk	4661	141.438	76.653	29.626	80.376	130.306	191.580	306.629
valuation_subsidy	4661	0.580	0.633	0.000	0.074	0.346	0.876	2.179
wholesale_ratio	4661	0.007	0.014	0.000	0.000	0.000	0.005	0.051

Summary by BTFP Usage

# Summary by BTFP usage
summary_by_btfp <- analysis_clean %>%
  group_by(btfp_user) %>%
  summarise(
    N = n(),
    `% Uninsured` = mean(pct_uninsured, na.rm = TRUE),
    `% MTM Loss` = mean(pct_mtm_loss, na.rm = TRUE),
    `Run Risk` = mean(run_risk, na.rm = TRUE),
    `Insured Deposit Coverage Ratio` = mean(insured_coverage, na.rm = TRUE),
    `Valuation Subsidy` = mean(valuation_subsidy, na.rm = TRUE),
    `OMO Eligible Ratio` = mean(omo_eligible_ratio, na.rm = TRUE),
    `FHLB Ratio` = mean(fhlb_ratio, na.rm = TRUE),
    `Log Assets` = mean(log_assets, na.rm = TRUE),
    `% MTM Insolvent` = 100 * mean(mtm_insolvent, na.rm = TRUE),
    .groups = "drop"
  ) %>%
  mutate(btfp_user = ifelse(btfp_user == 1, "BTFP User", "Non-User"))

kable(summary_by_btfp, 
      caption = "Table 2: Summary Statistics by BTFP Usage",
      digits = 3) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
                full_width = FALSE)

Table 2: Summary Statistics by BTFP Usage
btfp_user	N	% Uninsured	% MTM Loss	Run Risk	Insured Deposit Coverage Ratio	Valuation Subsidy	OMO Eligible Ratio	FHLB Ratio	Log Assets	% MTM Insolvent
Non-User	4188	26.847	5.362	136.932	-8.499	0.558	9.125	0.024	12.767	17.407
BTFP User	473	31.280	6.064	181.331	-10.314	0.770	10.306	0.036	13.772	28.541

Summary by OMO-Eligible Status

# Summary by OMO-eligible status
summary_by_omo <- analysis_clean %>%
  group_by(has_omo_eligible) %>%
  summarise(
    N = n(),
    `% Used BTFP` = 100 * mean(btfp_user, na.rm = TRUE),
    `% Used DW` = 100 * mean(dw_user, na.rm = TRUE),
    `% MTM Insolvent` = 100 * mean(mtm_insolvent, na.rm = TRUE),
    `Run Risk` = mean(run_risk, na.rm = TRUE),
    `Insured Coverage` = mean(insured_coverage, na.rm = TRUE),
    `Valuation Subsidy` = mean(valuation_subsidy, na.rm = TRUE),
    .groups = "drop"
  ) %>%
  mutate(has_omo_eligible = ifelse(has_omo_eligible, "Has OMO-Eligible", "No OMO-Eligible"))

kable(summary_by_omo, 
      caption = "Table 3: Summary Statistics by OMO-Eligible Status",
      digits = 3) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
                full_width = FALSE)

Table 3: Summary Statistics by OMO-Eligible Status
has_omo_eligible	N	% Used BTFP	% Used DW	% MTM Insolvent	Run Risk	Insured Coverage	Valuation Subsidy
No OMO-Eligible	386	3.886	4.663	10.104	114.245	-8.013	0.006
Has OMO-Eligible	4275	10.713	9.099	19.298	143.893	-8.744	0.631

Summary by Insolvency Status

# Summary by MTM insolvency status
summary_by_insolvency <- analysis_clean %>%
  group_by(mtm_insolvent) %>%
  summarise(
    N = n(),
    `% Used BTFP` = 100 * mean(btfp_user, na.rm = TRUE),
    `% Used DW` = 100 * mean(dw_user, na.rm = TRUE),
    `% Has OMO-Eligible` = 100 * mean(has_omo_eligible, na.rm = TRUE),
    `Run Risk` = mean(run_risk, na.rm = TRUE),
    `Valuation Subsidy` = mean(valuation_subsidy, na.rm = TRUE),
    `% Uninsured` = mean(pct_uninsured, na.rm = TRUE),
    .groups = "drop"
  ) %>%
  mutate(mtm_insolvent = ifelse(mtm_insolvent, "MTM Insolvent", "MTM Solvent"))

kable(summary_by_insolvency, 
      caption = "Table 4: Summary Statistics by MTM Insolvency Status",
      digits = 3) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
                full_width = FALSE)

Table 4: Summary Statistics by MTM Insolvency Status
mtm_insolvent	N	% Used BTFP	% Used DW	% Has OMO-Eligible	Run Risk	Valuation Subsidy	% Uninsured
MTM Solvent	3797	8.902	8.665	90.861	128.703	0.512	27.489
MTM Insolvent	864	15.625	9.028	95.486	197.403	0.877	26.454

Main Regression Analysis

Model Specification

Our baseline model building strategy:

Model	Specification	Purpose
1	Individual controls only	Baseline effects
2	Add Run Risk (liquidity)	Test liquidity channel
3	Add Insured Deposit Coverage Ratio (insolvency)	Test insolvency channel
4	Add Valuation Subsidy	Test BTFP incentive
5	Add Interaction: Insolvency × Subsidy	Key test of hypothesis
6	Full model with all controls	Robustness

################################################################################
# SECTION 8: MAIN REGRESSION ANALYSIS
################################################################################

cat("\n=== ESTIMATING MAIN MODELS ===\n\n")

#> 
#> === ESTIMATING MAIN MODELS ===

# Set up variable dictionary for nice labels
setFixest_dict(c(
  # Key variables
  pct_uninsured_z = "% Uninsured (z)",
  pct_mtm_loss_z = "% MTM Loss (z)",
  run_risk_z = "Run Risk (z)",
  insured_coverage_z = "Insured Coverage (z)",
  valuation_subsidy_z = "Valuation Subsidy (z)",
  omo_eligible_z = "OMO-Eligible Ratio (z)",
  
  # Interactions
  "insured_coverage_z:valuation_subsidy_z" = "Insolvency × Subsidy",
  "run_risk_z:valuation_subsidy_z" = "Run Risk × Subsidy",
  "run_risk_z:omo_eligible_z" = "Run Risk × OMO-Eligible",
  "mtm_insolventTRUE" = "MTM Insolvent",
  "mtm_insolventTRUE:valuation_subsidy_z" = "Insolvent × Subsidy",
  
  # Controls
  fhlb_ratio_z = "FHLB Ratio (z)",
  fhlb_st_ratio_z = "FHLB Short-term (z)",
  wholesale_ratio_z = "Wholesale Funding (z)",
  cash_ratio_z = "Cash Ratio (z)",
  book_equity_ratio_z = "Book Equity Ratio (z)",
  log_assets_z = "Log Assets (z)",
  "has_omo_eligibleTRUE" = "Has OMO-Eligible"
))

Extensive Margin Models (Logit)

Model Set 1: BTFP Usage

################################################################################
# MODEL SET 1: BTFP USAGE (LOGIT) - PROGRESSIVE BUILDING
################################################################################

# Restrict to banks with OMO-eligible assets (treatment pool)
analysis_omo <- analysis_clean %>% filter(has_omo_eligible)

cat("Analysis sample (OMO-eligible banks):", nrow(analysis_omo), "\n")

#> Analysis sample (OMO-eligible banks): 4275

cat("BTFP users in sample:", sum(analysis_omo$btfp_user), "\n\n")

#> BTFP users in sample: 458

# Model 1: Individual effects only
m1_btfp <- feglm(
  btfp_user ~ pct_uninsured_z + pct_mtm_loss_z + 
    log_assets_z,
  data = analysis_omo, 
  family = binomial("logit"), 
  vcov = "hetero"
)

# Model 2: Add Run Risk (liquidity channel)
m2_btfp <- feglm(
  btfp_user ~ run_risk_z + 
    log_assets_z,
  data = analysis_omo, 
  family = binomial("logit"), 
  vcov = "hetero"
)

# Model 3: Add Insured Coverage (insolvency channel)
m3_btfp <- feglm(
  btfp_user ~ run_risk_z + insured_coverage_z +
    log_assets_z,
  data = analysis_omo, 
  family = binomial("logit"), 
  vcov = "hetero"
)

# Model 4: Add Valuation Subsidy
m4_btfp <- feglm(
  btfp_user ~ run_risk_z + insured_coverage_z + valuation_subsidy_z +
    log_assets_z,
  data = analysis_omo, 
  family = binomial("logit"), 
  vcov = "hetero"
)

# Model 5: KEY - Interaction of Insolvency × Subsidy
m5_btfp <- feglm(
  btfp_user ~ run_risk_z + insured_coverage_z * valuation_subsidy_z +
    log_assets_z,
  data = analysis_omo, 
  family = binomial("logit"), 
  vcov = "hetero"
)

# Model 6: Full model with all controls
m6_btfp <- feglm(
  btfp_user ~ run_risk_z + insured_coverage_z * valuation_subsidy_z +
    fhlb_ratio_z + wholesale_ratio_z + cash_ratio_z + book_equity_ratio_z +
    log_assets_z,
  data = analysis_omo, 
  family = binomial("logit"), 
  vcov = "hetero"
)

cat("✓ BTFP logit models estimated\n")

#> ✓ BTFP logit models estimated

# Display results
etable(m1_btfp, m2_btfp, m3_btfp, m4_btfp, m5_btfp, m6_btfp,
       title = "Table 5: BTFP Usage - Logit Models (OMO-Eligible Banks Only)",
       headers = c("(1) Base", "(2) +Run Risk", "(3) +Coverage", 
                   "(4) +Subsidy", "(5) +Interaction", "(6) Full"),
       notes = "Heteroskedasticity-robust standard errors. Sample restricted to banks with OMO-eligible assets.",
       fitstat = ~ n + pr2)

#>                                  m1_btfp            m2_btfp             m3_btfp
#>                                 (1) Base      (2) +Run Risk       (3) +Coverage
#> Dependent Var.:                btfp_user          btfp_user           btfp_user
#>                                                                                
#> Constant              -2.345*** (0.0568) -2.330*** (0.0559)  -2.337*** (0.0560)
#> % Uninsured (z)       0.2053*** (0.0566)                                       
#> % MTM Loss (z)        0.4095*** (0.0514)                                       
#> Log Assets (z)        0.4959*** (0.0499) 0.4386*** (0.0460)  0.4535*** (0.0471)
#> Run Risk (z)                             0.3805*** (0.0492)     0.1395 (0.0860)
#> Insured Coverage (z)                                        -0.2915*** (0.0859)
#> Valuation Subsidy (z)                                                          
#> Insolvency × Subsidy                                                          
#> FHLB Ratio (z)                                                                 
#> Wholesale Funding (z)                                                          
#> Cash Ratio (z)                                                                 
#> Book Equity Ratio (z)                                                          
#> _____________________ __________________ __________________ ___________________
#> S.E. type             Heteroskedas.-rob. Heteroskedas.-rob. Heteroskedast.-rob.
#> Observations                       4,275              4,275               4,275
#> Pseudo R2                        0.07243            0.06987             0.07312
#> 
#>                                  m4_btfp            m5_btfp             m6_btfp
#>                             (4) +Subsidy   (5) +Interaction            (6) Full
#> Dependent Var.:                btfp_user          btfp_user           btfp_user
#>                                                                                
#> Constant              -2.347*** (0.0561) -2.344*** (0.0566)  -2.575*** (0.0773)
#> % Uninsured (z)                                                                
#> % MTM Loss (z)                                                                 
#> Log Assets (z)        0.4369*** (0.0479) 0.4424*** (0.0479)  0.4141*** (0.0556)
#> Run Risk (z)            0.1435. (0.0868)   0.1489. (0.0870)    -0.0189 (0.1122)
#> Insured Coverage (z)  -0.2641** (0.0878) -0.2642** (0.0879)   -0.2195* (0.1089)
#> Valuation Subsidy (z)   0.1191* (0.0482)  0.1498** (0.0539)    0.1041. (0.0587)
#> Insolvency × Subsidy                       0.0526 (0.0474)     0.0383 (0.0502)
#> FHLB Ratio (z)                                               0.2080*** (0.0539)
#> Wholesale Funding (z)                                         0.1441** (0.0462)
#> Cash Ratio (z)                                              -0.5452*** (0.1463)
#> Book Equity Ratio (z)                                       -0.3349*** (0.0735)
#> _____________________ __________________ __________________ ___________________
#> S.E. type             Heteroskedas.-rob. Heteroskedas.-rob. Heteroskedast.-rob.
#> Observations                       4,275              4,275               4,270
#> Pseudo R2                        0.07514            0.07551             0.10978
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Model Set 2: DW Usage (Comparison)

################################################################################
# MODEL SET 2: DW USAGE (LOGIT) - FOR COMPARISON
################################################################################

# Model 5 equivalent for DW
m5_dw <- feglm(
  dw_user ~ run_risk_z + insured_coverage_z * valuation_subsidy_z +
    log_assets_z,
  data = analysis_omo, 
  family = binomial("logit"), 
  vcov = "hetero"
)

# Full model for DW
m6_dw <- feglm(
  dw_user ~ run_risk_z + insured_coverage_z * valuation_subsidy_z +
    fhlb_ratio_z + wholesale_ratio_z + cash_ratio_z + book_equity_ratio_z +
    log_assets_z,
  data = analysis_omo, 
  family = binomial("logit"), 
  vcov = "hetero"
)

cat("✓ DW logit models estimated\n")

#> ✓ DW logit models estimated

# Compare BTFP vs DW
etable(m5_btfp, m6_btfp, m5_dw, m6_dw,
       title = "Table 6: BTFP vs DW Usage - Logit Models",
       headers = c("BTFP (5)", "BTFP (6)", "DW (5)", "DW (6)"),
       notes = "Heteroskedasticity-robust standard errors.",
       fitstat = ~ n + pr2)

#>                                  m5_btfp             m6_btfp              m5_dw
#>                                 BTFP (5)            BTFP (6)             DW (5)
#> Dependent Var.:                btfp_user           btfp_user            dw_user
#>                                                                                
#> Constant              -2.344*** (0.0566)  -2.575*** (0.0773) -2.553*** (0.0623)
#> Run Risk (z)            0.1489. (0.0870)    -0.0189 (0.1122)   0.2413* (0.1004)
#> Insured Coverage (z)  -0.2642** (0.0879)   -0.2195* (0.1089)   -0.0448 (0.1038)
#> Valuation Subsidy (z)  0.1498** (0.0539)    0.1041. (0.0587)    0.0083 (0.0655)
#> Log Assets (z)        0.4424*** (0.0479)  0.4141*** (0.0556) 0.6088*** (0.0492)
#> Insolvency × Subsidy    0.0526 (0.0474)     0.0383 (0.0502)    0.0667 (0.0579)
#> FHLB Ratio (z)                            0.2080*** (0.0539)                   
#> Wholesale Funding (z)                      0.1441** (0.0462)                   
#> Cash Ratio (z)                           -0.5452*** (0.1463)                   
#> Book Equity Ratio (z)                    -0.3349*** (0.0735)                   
#> _____________________ __________________ ___________________ __________________
#> S.E. type             Heteroskedas.-rob. Heteroskedast.-rob. Heteroskedas.-rob.
#> Observations                       4,275               4,270              4,275
#> Pseudo R2                        0.07551             0.10978            0.08741
#> 
#>                                    m6_dw
#>                                   DW (6)
#> Dependent Var.:                  dw_user
#>                                         
#> Constant              -2.602*** (0.0675)
#> Run Risk (z)            0.1926. (0.1159)
#> Insured Coverage (z)    -0.0138 (0.1130)
#> Valuation Subsidy (z)   -0.0169 (0.0685)
#> Log Assets (z)        0.5882*** (0.0529)
#> Insolvency × Subsidy    0.0612 (0.0589)
#> FHLB Ratio (z)           0.0701 (0.0570)
#> Wholesale Funding (z)  0.1373** (0.0485)
#> Cash Ratio (z)          -0.1612 (0.1019)
#> Book Equity Ratio (z)   -0.0908 (0.0715)
#> _____________________ __________________
#> S.E. type             Heteroskedas.-rob.
#> Observations                       4,270
#> Pseudo R2                        0.09342
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Model Set 3: Alternative Insolvency Measure (MTM Insolvent Indicator)

################################################################################
# MODEL SET 3: MTM INSOLVENCY INDICATOR
################################################################################

# Model with binary insolvency indicator
m_btfp_insolvent <- feglm(
  btfp_user ~ run_risk_z + mtm_insolvent * valuation_subsidy_z +
    fhlb_ratio_z + wholesale_ratio_z + cash_ratio_z + book_equity_ratio_z +
    log_assets_z,
  data = analysis_omo, 
  family = binomial("logit"), 
  vcov = "hetero"
)

m_dw_insolvent <- feglm(
  dw_user ~ run_risk_z + mtm_insolvent * valuation_subsidy_z +
    fhlb_ratio_z + wholesale_ratio_z + cash_ratio_z + book_equity_ratio_z +
    log_assets_z,
  data = analysis_omo, 
  family = binomial("logit"), 
  vcov = "hetero"
)

cat("✓ MTM Insolvency indicator models estimated\n")

#> ✓ MTM Insolvency indicator models estimated

etable(m_btfp_insolvent, m_dw_insolvent,
       title = "Table 7: Facility Usage with MTM Insolvency Indicator",
       headers = c("BTFP", "DW"),
       notes = "mtm_insolvent = 1 if MV Equity < 0.",
       fitstat = ~ n + pr2)

#>                          m_btfp_insolvent     m_dw_insolvent
#>                                      BTFP                 DW
#> Dependent Var.:                 btfp_user            dw_user
#>                                                             
#> Constant               -2.580*** (0.0782) -2.563*** (0.0725)
#> Run Risk (z)             0.1569* (0.0640) 0.2306*** (0.0687)
#> MTM Insolvent             0.0578 (0.1635)   -0.2800 (0.1847)
#> Valuation Subsidy (z)     0.0958 (0.0621)   -0.0408 (0.0694)
#> FHLB Ratio (z)         0.2442*** (0.0505)    0.0751 (0.0548)
#> Wholesale Funding (z)  0.1527*** (0.0457)  0.1349** (0.0483)
#> Cash Ratio (z)        -0.5132*** (0.1401)   -0.1513 (0.1027)
#> Book Equity Ratio (z) -0.3194*** (0.0860)  -0.1588. (0.0846)
#> Log Assets (z)         0.3972*** (0.0562) 0.5655*** (0.0539)
#> Insolvent × Subsidy     -0.0022 (0.1014)   -0.0008 (0.1250)
#> _____________________ ___________________ __________________
#> S.E. type             Heteroskedast.-rob. Heteroskedas.-rob.
#> Observations                        4,270              4,270
#> Pseudo R2                         0.10830            0.09398
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Linear Probability Models (LPM) for Robustness

################################################################################
# MODEL SET 4: LINEAR PROBABILITY MODELS (ROBUSTNESS)
################################################################################

# LPM version of key model
m_btfp_lpm <- feols(
  btfp_user ~ run_risk_z + insured_coverage_z * valuation_subsidy_z +
    fhlb_ratio_z + wholesale_ratio_z + cash_ratio_z + book_equity_ratio_z +
    log_assets_z,
  data = analysis_omo, 
  vcov = "hetero"
)

m_dw_lpm <- feols(
  dw_user ~ run_risk_z + insured_coverage_z * valuation_subsidy_z +
    fhlb_ratio_z + wholesale_ratio_z + cash_ratio_z + book_equity_ratio_z +
    log_assets_z,
  data = analysis_omo, 
  vcov = "hetero"
)

cat("✓ LPM models estimated\n")

#> ✓ LPM models estimated

etable(m_btfp_lpm, m_dw_lpm, m6_btfp, m6_dw,
       title = "Table 8: LPM vs Logit Comparison",
       headers = c("BTFP LPM", "DW LPM", "BTFP Logit", "DW Logit"),
       notes = "LPM coefficients are directly interpretable as probability changes.",
       fitstat = ~ n + r2 + pr2)

#>                                m_btfp_lpm           m_dw_lpm
#>                                  BTFP LPM             DW LPM
#> Dependent Var.:                 btfp_user            dw_user
#>                                                             
#> Constant               0.1009*** (0.0046) 0.0883*** (0.0043)
#> Run Risk (z)              0.0117 (0.0096)   0.0189* (0.0084)
#> Insured Coverage (z)     -0.0132 (0.0083)    0.0007 (0.0073)
#> Valuation Subsidy (z)    0.0120* (0.0048)    0.0026 (0.0045)
#> FHLB Ratio (z)         0.0253*** (0.0058)    0.0073 (0.0050)
#> Wholesale Funding (z)   0.0169** (0.0053)  0.0127** (0.0048)
#> Cash Ratio (z)          -0.0079. (0.0045)    0.0047 (0.0045)
#> Book Equity Ratio (z) -0.0200*** (0.0047)   -0.0029 (0.0045)
#> Log Assets (z)         0.0425*** (0.0062) 0.0603*** (0.0061)
#> Insolvency × Subsidy    -0.0057 (0.0053)    0.0028 (0.0049)
#> _____________________ ___________________ __________________
#> Family                                OLS                OLS
#> S.E. type             Heteroskedast.-rob. Heteroskedas.-rob.
#> Observations                        4,270              4,270
#> R2                                0.07048            0.06243
#> Pseudo R2                         0.14917            0.18602
#> 
#>                                   m6_btfp              m6_dw
#>                                BTFP Logit           DW Logit
#> Dependent Var.:                 btfp_user            dw_user
#>                                                             
#> Constant               -2.575*** (0.0773) -2.602*** (0.0675)
#> Run Risk (z)             -0.0189 (0.1122)   0.1926. (0.1159)
#> Insured Coverage (z)    -0.2195* (0.1089)   -0.0138 (0.1130)
#> Valuation Subsidy (z)    0.1041. (0.0587)   -0.0169 (0.0685)
#> FHLB Ratio (z)         0.2080*** (0.0539)    0.0701 (0.0570)
#> Wholesale Funding (z)   0.1441** (0.0462)  0.1373** (0.0485)
#> Cash Ratio (z)        -0.5452*** (0.1463)   -0.1612 (0.1019)
#> Book Equity Ratio (z) -0.3349*** (0.0735)   -0.0908 (0.0715)
#> Log Assets (z)         0.4141*** (0.0556) 0.5882*** (0.0529)
#> Insolvency × Subsidy     0.0383 (0.0502)    0.0612 (0.0589)
#> _____________________ ___________________ __________________
#> Family                              Logit              Logit
#> S.E. type             Heteroskedast.-rob. Heteroskedas.-rob.
#> Observations                        4,270              4,270
#> R2                                     --                 --
#> Pseudo R2                         0.10978            0.09342
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Intensive Margin Models

################################################################################
# MODEL SET 5: INTENSIVE MARGIN (AMOUNT BORROWED)
################################################################################

# BTFP amount (users only)
btfp_users <- analysis_omo %>% filter(btfp_user == 1)

if (nrow(btfp_users) >= 30) {
  m_btfp_intensive <- feols(
    btfp_amt_log ~ run_risk_z + insured_coverage_z * valuation_subsidy_z +
      fhlb_ratio_z + wholesale_ratio_z + cash_ratio_z + book_equity_ratio_z +
      log_assets_z,
    data = btfp_users, 
    vcov = "hetero"
  )
  
  cat("✓ BTFP intensive margin model estimated (n =", nrow(btfp_users), ")\n")
  
  etable(m_btfp_intensive,
         title = "Table 9: BTFP Intensive Margin (Users Only)",
         notes = "Dependent variable: log(1 + BTFP Amount / Assets)",
         fitstat = ~ n + r2)
} else {
  cat("⚠ Insufficient BTFP users for intensive margin analysis\n")
}

#> ✓ BTFP intensive margin model estimated (n = 458 )

#>                          m_btfp_intensive
#> Dependent Var.:              btfp_amt_log
#>                                          
#> Constant                9.642*** (0.2279)
#> Run Risk (z)             -0.3990 (0.2458)
#> Insured Coverage (z)    -0.5207* (0.2584)
#> Valuation Subsidy (z)     0.0335 (0.1301)
#> FHLB Ratio (z)          0.3728** (0.1279)
#> Wholesale Funding (z)  0.3132*** (0.0837)
#> Cash Ratio (z)           -0.6327 (0.3850)
#> Book Equity Ratio (z)    -0.2278 (0.1643)
#> Log Assets (z)        -0.9522*** (0.1637)
#> Insolvency × Subsidy    -0.0669 (0.1340)
#> _____________________ ___________________
#> S.E. type             Heteroskedast.-rob.
#> Observations                          457
#> R2                                0.19765
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Conditional Choice Model (BTFP vs DW)

################################################################################
# MODEL SET 6: CONDITIONAL CHOICE (BTFP vs DW among single-facility users)
################################################################################

# Sample: Banks that used exactly one facility
choice_sample <- analysis_omo %>%
  filter(used_only_btfp | used_only_dw)

cat("Conditional choice sample size:", nrow(choice_sample), "\n")

#> Conditional choice sample size: 663

cat("  Chose BTFP:", sum(choice_sample$chose_btfp, na.rm = TRUE), "\n")

#>   Chose BTFP: 366

cat("  Chose DW:", sum(choice_sample$chose_btfp == 0, na.rm = TRUE), "\n\n")

#>   Chose DW: 297

if (nrow(choice_sample) >= 50) {
  m_choice_logit <- feglm(
    chose_btfp ~ run_risk_z + insured_coverage_z * valuation_subsidy_z +
      fhlb_ratio_z + wholesale_ratio_z + cash_ratio_z + book_equity_ratio_z +
      log_assets_z,
    data = choice_sample, 
    family = binomial("logit"), 
    vcov = "hetero"
  )
  
  m_choice_lpm <- feols(
    chose_btfp ~ run_risk_z + insured_coverage_z * valuation_subsidy_z +
      fhlb_ratio_z + wholesale_ratio_z + cash_ratio_z + book_equity_ratio_z +
      log_assets_z,
    data = choice_sample, 
    vcov = "hetero"
  )
  
  cat("✓ Conditional choice models estimated\n")
  
  etable(m_choice_logit, m_choice_lpm,
         title = "Table 10: Conditional Choice - BTFP vs DW",
         headers = c("Logit", "LPM"),
         notes = "Sample: Banks that used exactly one facility. Dependent variable: chose_btfp = 1 if used only BTFP.",
         fitstat = ~ n + pr2 + r2)
} else {
  cat("⚠ Insufficient observations for conditional choice analysis\n")
}

#> ✓ Conditional choice models estimated

#>                          m_choice_logit       m_choice_lpm
#>                                   Logit                LPM
#> Dependent Var.:              chose_btfp         chose_btfp
#>                                                           
#> Constant                0.0653 (0.1096) 0.5173*** (0.0254)
#> Run Risk (z)           -0.2175 (0.1857)   -0.0505 (0.0424)
#> Insured Coverage (z)   -0.2115 (0.1798)   -0.0492 (0.0411)
#> Valuation Subsidy (z)   0.0954 (0.1017)    0.0220 (0.0236)
#> FHLB Ratio (z)         0.1646. (0.0872)   0.0385. (0.0198)
#> Wholesale Funding (z)   0.0064 (0.0723)    0.0017 (0.0168)
#> Cash Ratio (z)        -0.3382. (0.1768)  -0.0768* (0.0380)
#> Book Equity Ratio (z) -0.2698* (0.1059) -0.0637** (0.0246)
#> Log Assets (z)        -0.2129* (0.0970)  -0.0491* (0.0218)
#> Insolvency × Subsidy  -0.0540 (0.0910)   -0.0120 (0.0206)
#> _____________________ _________________ __________________
#> Family                            Logit                OLS
#> S.E. type             Heteroskeda.-rob. Heteroskedas.-rob.
#> Observations                        662                662
#> Pseudo R2                       0.03708            0.03523
#> R2                                   --            0.04950
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Robustness Checks

Among insolvent banks, does subsidy matter more?

# Among insolvent banks, does subsidy matter more?
m_insolvent_only <- feglm(
  btfp_user ~ run_risk_z + valuation_subsidy_z + log_assets_z,
  data = analysis_omo %>% filter(mtm_insolvent == TRUE),
  family = binomial("logit"), vcov = "hetero"
)

m_solvent_only <- feglm(
  btfp_user ~ run_risk_z + valuation_subsidy_z + log_assets_z,
  data = analysis_omo %>% filter(mtm_insolvent == FALSE),
  family = binomial("logit"), vcov = "hetero"
)

etable(m_insolvent_only, m_solvent_only, headers = c("Insolvent", "Solvent"))

#>                         m_insolvent_only     m_solvent_only
#>                                Insolvent            Solvent
#> Dependent Var.:                btfp_user          btfp_user
#>                                                            
#> Constant              -1.868*** (0.1453) -2.433*** (0.0650)
#> Run Risk (z)             0.1177 (0.1099) 0.3631*** (0.0626)
#> Valuation Subsidy (z)    0.1010 (0.0851)    0.0978 (0.0605)
#> Log Assets (z)        0.4932*** (0.1233) 0.4492*** (0.0539)
#> _____________________ __________________ __________________
#> S.E. type             Heteroskedas.-rob. Heteroskedas.-rob.
#> Observations                         825              3,450
#> Squared Cor.                     0.02822            0.05144
#> Pseudo R2                        0.03796            0.07923
#> BIC                               718.36            2,028.3
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Test if the LEVEL difference is significant
m_pooled_test <- feglm(
  btfp_user ~ mtm_insolvent + run_risk_z + valuation_subsidy_z + log_assets_z,
  data = analysis_omo,
  family = binomial("logit"), vcov = "hetero"
)

etable(m_pooled_test)

#>                            m_pooled_test
#> Dependent Var.:                btfp_user
#>                                         
#> Constant              -2.425*** (0.0635)
#> MTM Insolvent          0.3904** (0.1326)
#> Run Risk (z)          0.2989*** (0.0551)
#> Valuation Subsidy (z)   0.1053* (0.0490)
#> Log Assets (z)        0.4621*** (0.0488)
#> _____________________ __________________
#> S.E. type             Heteroskedas.-rob.
#> Observations                       4,275
#> Squared Cor.                     0.05038
#> Pseudo R2                        0.07561
#> BIC                              2,732.8
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Final Model Comparison Table
etable(m2_btfp, m5_btfp, m_pooled_test, m6_btfp,
       headers = c("Run Risk Only", "Continuous Insolvency", "Binary Insolvency", "Full Controls"),
       title = "Table: Determinants of BTFP Usage")

#>                                  m2_btfp               m5_btfp
#>                            Run Risk Only Continuous Insolvency
#> Dependent Var.:                btfp_user             btfp_user
#>                                                               
#> Constant              -2.330*** (0.0559)    -2.344*** (0.0566)
#> Run Risk (z)          0.3805*** (0.0492)      0.1489. (0.0870)
#> Log Assets (z)        0.4386*** (0.0460)    0.4424*** (0.0479)
#> Insured Coverage (z)                        -0.2642** (0.0879)
#> Valuation Subsidy (z)                        0.1498** (0.0539)
#> Insolvency × Subsidy                          0.0526 (0.0474)
#> MTM Insolvent                                                 
#> FHLB Ratio (z)                                                
#> Wholesale Funding (z)                                         
#> Cash Ratio (z)                                                
#> Book Equity Ratio (z)                                         
#> _____________________ __________________    __________________
#> S.E. type             Heteroskedas.-rob.    Heteroskedas.-rob.
#> Observations                       4,275                 4,275
#> Squared Cor.                     0.04699               0.04991
#> Pseudo R2                        0.06987               0.07551
#> BIC                              2,732.8               2,741.5
#> 
#>                            m_pooled_test             m6_btfp
#>                        Binary Insolvency       Full Controls
#> Dependent Var.:                btfp_user           btfp_user
#>                                                             
#> Constant              -2.425*** (0.0635)  -2.575*** (0.0773)
#> Run Risk (z)          0.2989*** (0.0551)    -0.0189 (0.1122)
#> Log Assets (z)        0.4621*** (0.0488)  0.4141*** (0.0556)
#> Insured Coverage (z)                       -0.2195* (0.1089)
#> Valuation Subsidy (z)   0.1053* (0.0490)    0.1041. (0.0587)
#> Insolvency × Subsidy                        0.0383 (0.0502)
#> MTM Insolvent          0.3904** (0.1326)                    
#> FHLB Ratio (z)                            0.2080*** (0.0539)
#> Wholesale Funding (z)                      0.1441** (0.0462)
#> Cash Ratio (z)                           -0.5452*** (0.1463)
#> Book Equity Ratio (z)                    -0.3349*** (0.0735)
#> _____________________ __________________ ___________________
#> S.E. type             Heteroskedas.-rob. Heteroskedast.-rob.
#> Observations                       4,275               4,270
#> Squared Cor.                     0.05038             0.07621
#> Pseudo R2                        0.07561             0.10978
#> BIC                              2,732.8             2,670.4
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Table: Determinants of BTFP Usage (Logit Models)

Four specifications to build our understanding:

Model 1 (Run Risk Only): Liquidity risk strongly predicts BTFP usage (β = 0.367, p < 0.001)

Model 2 (Continuous Insolvency): Adding insolvency measures reveals both channels matter:

Run Risk remains significant but attenuates (β = 0.202, p < 0.05) Insured Coverage is negative and significant (β = -0.187, p < 0.05), indicating banks closer to insolvency used BTFP more Valuation Subsidy is positive and significant (β = 0.139, p < 0.01), confirming BTFP’s par valuation attracted banks with larger MTM losses on eligible securities

##Model 3 (Binary Insolvency — Preferred Specification): The clearest result emerges with a binary insolvency indicator:

MTM Insolvent banks were 50% more likely to use BTFP (β = 0.406, OR = 1.50, p < 0.01) Run Risk remains highly significant (β = 0.290, p < 0.001) This parsimonious model captures the key economic mechanism

##Model 4 (Full Controls): Adding balance sheet controls improves fit (Pseudo R² = 0.105) but absorbs significance from risk variables:

FHLB Ratio (+): Banks already using wholesale funding also accessed BTFP Book Equity Ratio (−): Better capitalized banks less likely to use BTFP These controls may over-control for the insolvency channel

Alternative Sample: Full Sample (Including Non-OMO Banks)

################################################################################
# ROBUSTNESS: FULL SAMPLE WITH OMO-ELIGIBLE CONTROL
################################################################################

cat("\n=== ROBUSTNESS CHECKS ===\n\n")

#> 
#> === ROBUSTNESS CHECKS ===

# Full sample with has_omo_eligible as control
m_full_btfp <- feglm(
  btfp_user ~ run_risk_z + insured_coverage_z * valuation_subsidy_z +
    has_omo_eligible +
    fhlb_ratio_z + wholesale_ratio_z + cash_ratio_z + book_equity_ratio_z +
    log_assets_z,
  data = analysis_clean, 
  family = binomial("logit"), 
  vcov = "hetero"
)

# Full sample with interaction: Run Risk × OMO-Eligible
m_full_interact <- feglm(
  btfp_user ~ run_risk_z * has_omo_eligible + 
    insured_coverage_z * valuation_subsidy_z +
    fhlb_ratio_z + wholesale_ratio_z + cash_ratio_z + book_equity_ratio_z +
    log_assets_z,
  data = analysis_clean, 
  family = binomial("logit"), 
  vcov = "hetero"
)

cat("✓ Full sample robustness models estimated\n")

#> ✓ Full sample robustness models estimated

etable(m6_btfp, m_full_btfp, m_full_interact,
       title = "Table 11: Robustness - OMO-Eligible Sample vs Full Sample",
       headers = c("OMO Only", "Full + Control", "Full + Interaction"),
       fitstat = ~ n + pr2)

#>                                             m6_btfp         m_full_btfp
#>                                            OMO Only      Full + Control
#> Dependent Var.:                           btfp_user           btfp_user
#>                                                                        
#> Constant                         -2.575*** (0.0773)  -2.990*** (0.2749)
#> Run Risk (z)                       -0.0189 (0.1122)    -0.0311 (0.1093)
#> Insured Coverage (z)              -0.2195* (0.1089)   -0.2338* (0.1062)
#> Valuation Subsidy (z)              0.1041. (0.0587)    0.0985. (0.0586)
#> FHLB Ratio (z)                   0.2080*** (0.0539)  0.1907*** (0.0526)
#> Wholesale Funding (z)             0.1441** (0.0462)  0.1558*** (0.0456)
#> Cash Ratio (z)                  -0.5452*** (0.1463) -0.5641*** (0.1407)
#> Book Equity Ratio (z)           -0.3349*** (0.0735) -0.3511*** (0.0730)
#> Log Assets (z)                   0.4141*** (0.0556)  0.4124*** (0.0550)
#> Insolvency × Subsidy               0.0383 (0.0502)     0.0412 (0.0493)
#> Has OMO-Eligible                                        0.4059 (0.2818)
#> Run Risk (z) x Has OMO-Eligible                                        
#> _______________________________ ___________________ ___________________
#> S.E. type                       Heteroskedast.-rob. Heteroskedast.-rob.
#> Observations                                  4,270               4,651
#> Pseudo R2                                   0.10978             0.11706
#> 
#>                                     m_full_interact
#>                                  Full + Interaction
#> Dependent Var.:                           btfp_user
#>                                                    
#> Constant                         -3.018*** (0.2884)
#> Run Risk (z)                        0.0778 (0.3118)
#> Insured Coverage (z)              -0.2330* (0.1064)
#> Valuation Subsidy (z)              0.0972. (0.0587)
#> FHLB Ratio (z)                   0.1908*** (0.0526)
#> Wholesale Funding (z)            0.1556*** (0.0456)
#> Cash Ratio (z)                  -0.5631*** (0.1409)
#> Book Equity Ratio (z)           -0.3506*** (0.0729)
#> Log Assets (z)                   0.4128*** (0.0550)
#> Insolvency × Subsidy               0.0380 (0.0502)
#> Has OMO-Eligible                    0.4355 (0.2955)
#> Run Risk (z) x Has OMO-Eligible    -0.1122 (0.2960)
#> _______________________________ ___________________
#> S.E. type                       Heteroskedast.-rob.
#> Observations                                  4,651
#> Pseudo R2                                   0.11710
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Alternative Severity Measure (50% Run Scenario)

################################################################################
# ROBUSTNESS: 50% UNINSURED RUN SCENARIO
################################################################################

# Create 50% run coverage variable
analysis_omo <- analysis_omo %>%
  mutate(insured_coverage_50_z = standardize(insured_coverage_50pct))

m_btfp_50pct <- feglm(
  btfp_user ~ run_risk_z + insured_coverage_50_z * valuation_subsidy_z +
    fhlb_ratio_z + wholesale_ratio_z + cash_ratio_z + book_equity_ratio_z +
    log_assets_z,
  data = analysis_omo, 
  family = binomial("logit"), 
  vcov = "hetero"
)

cat("✓ 50% run scenario model estimated\n")

#> ✓ 50% run scenario model estimated

etable(m6_btfp, m_btfp_50pct,
       title = "Table 12: Robustness - 100% vs 50% Uninsured Run",
       headers = c("100% Run", "50% Run"),
       fitstat = ~ n + pr2)

#>                                                           m6_btfp
#>                                                          100% Run
#> Dependent Var.:                                         btfp_user
#>                                                                  
#> Constant                                       -2.575*** (0.0773)
#> Run Risk (z)                                     -0.0189 (0.1122)
#> Insured Coverage (z)                            -0.2195* (0.1089)
#> Valuation Subsidy (z)                            0.1041. (0.0587)
#> FHLB Ratio (z)                                 0.2080*** (0.0539)
#> Wholesale Funding (z)                           0.1441** (0.0462)
#> Cash Ratio (z)                                -0.5452*** (0.1463)
#> Book Equity Ratio (z)                         -0.3349*** (0.0735)
#> Log Assets (z)                                 0.4141*** (0.0556)
#> Insolvency × Subsidy                             0.0383 (0.0502)
#> insured_coverage_50_z                                            
#> insured_coverage_50_z x Valuation Subsidy (z)                    
#> ________________________________________      ___________________
#> S.E. type                                     Heteroskedast.-rob.
#> Observations                                                4,270
#> Pseudo R2                                                 0.10978
#> 
#>                                                      m_btfp_50pct
#>                                                           50% Run
#> Dependent Var.:                                         btfp_user
#>                                                                  
#> Constant                                       -2.573*** (0.0771)
#> Run Risk (z)                                     -0.0175 (0.1080)
#> Insured Coverage (z)                                             
#> Valuation Subsidy (z)                            0.0987. (0.0586)
#> FHLB Ratio (z)                                 0.2063*** (0.0540)
#> Wholesale Funding (z)                           0.1435** (0.0462)
#> Cash Ratio (z)                                -0.5462*** (0.1460)
#> Book Equity Ratio (z)                         -0.3337*** (0.0735)
#> Log Assets (z)                                 0.4166*** (0.0558)
#> Insolvency × Subsidy                                            
#> insured_coverage_50_z                           -0.2172* (0.1029)
#> insured_coverage_50_z x Valuation Subsidy (z)     0.0321 (0.0497)
#> ________________________________________      ___________________
#> S.E. type                                     Heteroskedast.-rob.
#> Observations                                                4,270
#> Pseudo R2                                                 0.10985
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Size Subsamples

################################################################################
# ROBUSTNESS: BY SIZE BIN
################################################################################

if ("size_bin" %in% names(analysis_omo)) {
  # Small banks
  m_small <- tryCatch({
    feglm(
      btfp_user ~ run_risk_z + insured_coverage_z * valuation_subsidy_z +
        fhlb_ratio_z + wholesale_ratio_z + cash_ratio_z + book_equity_ratio_z,
      data = analysis_omo %>% filter(size_bin == "small"),
      family = binomial("logit"), 
      vcov = "hetero"
    )
  }, error = function(e) NULL)
  
  # Medium banks
  m_medium <- tryCatch({
    feglm(
      btfp_user ~ run_risk_z + insured_coverage_z * valuation_subsidy_z +
        fhlb_ratio_z + wholesale_ratio_z + cash_ratio_z + book_equity_ratio_z,
      data = analysis_omo %>% filter(size_bin == "medium"),
      family = binomial("logit"), 
      vcov = "hetero"
    )
  }, error = function(e) NULL)
  
  # Large banks
  m_large <- tryCatch({
    feglm(
      btfp_user ~ run_risk_z + insured_coverage_z * valuation_subsidy_z +
        fhlb_ratio_z + wholesale_ratio_z + cash_ratio_z + book_equity_ratio_z,
      data = analysis_omo %>% filter(size_bin == "large"),
      family = binomial("logit"), 
      vcov = "hetero"
    )
  }, error = function(e) NULL)
  
  models_by_size <- list()
  headers_size <- c()
  
  if (!is.null(m_small)) { models_by_size <- c(models_by_size, list(m_small)); headers_size <- c(headers_size, "small") }
  if (!is.null(m_medium)) { models_by_size <- c(models_by_size, list(m_medium)); headers_size <- c(headers_size, "medium") }
  if (!is.null(m_large)) { models_by_size <- c(models_by_size, list(m_large)); headers_size <- c(headers_size, "large") }
  
  if (length(models_by_size) > 0) {
    cat("✓ Size subsample models estimated\n")
    etable(models_by_size,
           title = "Table 13: Robustness - By Bank Size",
           headers = headers_size,
           fitstat = ~ n + pr2)
  }
}

#> ✓ Size subsample models estimated

#>                                   model 1            model 2
#>                                     small              large
#> Dependent Var.:                 btfp_user          btfp_user
#>                                                             
#> Constant               -2.712*** (0.0903) -1.814*** (0.1626)
#> Run Risk (z)              0.0849 (0.1293)    0.0858 (0.2359)
#> Insured Coverage (z)     -0.2062 (0.1300)   -0.0031 (0.2236)
#> Valuation Subsidy (z)    0.1532* (0.0636)    0.0639 (0.1258)
#> FHLB Ratio (z)         0.2091*** (0.0621)  0.3090** (0.0986)
#> Wholesale Funding (z)   0.1513** (0.0546)   0.1365. (0.0818)
#> Cash Ratio (z)        -0.8561*** (0.1678)   -0.1548 (0.2609)
#> Book Equity Ratio (z)  -0.2479** (0.0780)  -0.3029* (0.1477)
#> Insolvency × Subsidy    0.1125. (0.0633)  -0.1624. (0.0897)
#> _____________________ ___________________ __________________
#> S.E. type             Heteroskedast.-rob. Heteroskedas.-rob.
#> Observations                        3,533                713
#> Pseudo R2                         0.09032            0.05695
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Visualizations

BTFP Usage by Risk Characteristics

################################################################################
# VISUALIZATIONS
################################################################################

cat("\n=== CREATING VISUALIZATIONS ===\n\n")

#> 
#> === CREATING VISUALIZATIONS ===

# Usage rates by insolvency and subsidy quartiles
usage_plot_data <- analysis_omo %>%
  mutate(
    insolvency_quartile = ntile(insured_coverage, 4),
    subsidy_quartile = ntile(valuation_subsidy, 4)
  ) %>%
  group_by(insolvency_quartile, subsidy_quartile) %>%
  summarise(
    btfp_rate = mean(btfp_user) * 100,
    n = n(),
    .groups = "drop"
  ) %>%
  filter(!is.na(insolvency_quartile), !is.na(subsidy_quartile))

p1 <- ggplot(usage_plot_data, aes(x = factor(insolvency_quartile), 
                                   y = factor(subsidy_quartile), 
                                   fill = btfp_rate)) +
  geom_tile(color = "white", linewidth = 0.5) +
  geom_text(aes(label = sprintf("%.1f%%\n(n=%d)", btfp_rate, n)), 
            color = "white", size = 3.5) +
  scale_fill_gradient2(low = "#2166ac", mid = "#f7f7f7", high = "#b2182b",
                       midpoint = mean(usage_plot_data$btfp_rate),
                       name = "BTFP Usage %") +
  labs(
    title = "BTFP Usage Rate by Insolvency Risk and Valuation Subsidy",
    x = "Insured Coverage Quartile (1 = Most Insolvent)",
    y = "Valuation Subsidy Quartile (4 = Highest Subsidy)"
  ) +
  theme_minimal(base_size = 12) +
  theme(legend.position = "right")

print(p1)

# Save figure
ggsave(file.path(FIG_PATH, "btfp_usage_heatmap.png"), p1, 
       width = 10, height = 6, dpi = 300)

BTFP vs DW Usage Comparison

# Usage by insolvency status
usage_by_status <- analysis_omo %>%
  group_by(mtm_insolvent) %>%
  summarise(
    BTFP = mean(btfp_user) * 100,
    DW = mean(dw_user) * 100,
    `Any Fed` = mean(used_any_fed) * 100,
    n = n(),
    .groups = "drop"
  ) %>%
  pivot_longer(cols = c(BTFP, DW, `Any Fed`), 
               names_to = "Facility", values_to = "Usage_Rate") %>%
  mutate(Status = ifelse(mtm_insolvent, "MTM Insolvent", "MTM Solvent"))

p2 <- ggplot(usage_by_status, aes(x = Status, y = Usage_Rate, fill = Facility)) +
  geom_col(position = position_dodge(width = 0.8), width = 0.7) +
  geom_text(aes(label = sprintf("%.1f%%", Usage_Rate)), 
            position = position_dodge(width = 0.8), vjust = -0.5, size = 3.5) +
  scale_fill_manual(values = c("BTFP" = "#2166ac", "DW" = "#b2182b", "Any Fed" = "#7570b3")) +
  labs(
    title = "Facility Usage by MTM Solvency Status",
    subtitle = "Banks with OMO-eligible assets only",
    x = NULL,
    y = "Usage Rate (%)",
    fill = "Facility"
  ) +
  theme_minimal(base_size = 12) +
  theme(legend.position = "bottom") +
  ylim(0, max(usage_by_status$Usage_Rate) * 1.15)

print(p2)

ggsave(file.path(FIG_PATH, "usage_by_insolvency.png"), p2, 
       width = 8, height = 6, dpi = 300)

Coefficient Plot

# Extract coefficients from key models
extract_coefs <- function(model, model_name) {
  tidy(model, conf.int = TRUE) %>%
    mutate(model = model_name)
}

coef_df <- bind_rows(
  extract_coefs(m6_btfp, "BTFP"),
  extract_coefs(m6_dw, "DW")
) %>%
  filter(term %in% c("run_risk_z", "insured_coverage_z", "valuation_subsidy_z",
                     "insured_coverage_z:valuation_subsidy_z")) %>%
  mutate(
    term_label = case_when(
      term == "run_risk_z" ~ "Run Risk",
      term == "insured_coverage_z" ~ "Insured Coverage",
      term == "valuation_subsidy_z" ~ "Valuation Subsidy",
      term == "insured_coverage_z:valuation_subsidy_z" ~ "Coverage × Subsidy",
      TRUE ~ term
    ),
    term_label = factor(term_label, levels = c("Run Risk", "Insured Coverage", 
                                                "Valuation Subsidy", "Coverage × Subsidy"))
  )

p3 <- ggplot(coef_df, aes(x = term_label, y = estimate, color = model)) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray50") +
  geom_point(position = position_dodge(width = 0.5), size = 3) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high),
                position = position_dodge(width = 0.5), width = 0.2) +
  scale_color_manual(values = c("BTFP" = "#2166ac", "DW" = "#b2182b")) +
  labs(
    title = "Coefficient Estimates: BTFP vs DW",
    subtitle = "Logit models with 95% confidence intervals",
    x = NULL,
    y = "Coefficient (log-odds)",
    color = "Facility"
  ) +
  theme_minimal(base_size = 12) +
  theme(legend.position = "bottom",
        axis.text.x = element_text(angle = 15, hjust = 1))

print(p3)

ggsave(file.path(FIG_PATH, "coefficient_comparison.png"), p3, 
       width = 10, height = 6, dpi = 300)

Results Interpretation

Key Findings

################################################################################
# SECTION 9: RESULTS INTERPRETATION
################################################################################

cat("\n")

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

#> ================================================================

cat("KEY FINDINGS SUMMARY\n")

#> KEY FINDINGS SUMMARY

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

#> ================================================================

# Extract key coefficients from main model
main_results <- tidy(m6_btfp, conf.int = TRUE)

cat("1. MAIN MODEL (BTFP Usage - Logit, OMO-Eligible Banks):\n\n")

#> 1. MAIN MODEL (BTFP Usage - Logit, OMO-Eligible Banks):

# Run Risk
run_risk_coef <- main_results %>% filter(term == "run_risk_z")
if (nrow(run_risk_coef) > 0) {
  cat("   Run Risk (z):\n")
  cat("     Coefficient:", round(run_risk_coef$estimate, 3), "\n")
  cat("     95% CI: [", round(run_risk_coef$conf.low, 3), ",", 
      round(run_risk_coef$conf.high, 3), "]\n")
  cat("     p-value:", format.pval(run_risk_coef$p.value, digits = 3), "\n\n")
}

#>    Run Risk (z):
#>      Coefficient: -0.019 
#>      95% CI: [ -0.239 , 0.201 ]
#>      p-value: 0.867

# Insured Coverage
coverage_coef <- main_results %>% filter(term == "insured_coverage_z")
if (nrow(coverage_coef) > 0) {
  cat("   Insured Coverage (z):\n")
  cat("     Coefficient:", round(coverage_coef$estimate, 3), "\n")
  cat("     95% CI: [", round(coverage_coef$conf.low, 3), ",", 
      round(coverage_coef$conf.high, 3), "]\n")
  cat("     p-value:", format.pval(coverage_coef$p.value, digits = 3), "\n\n")
}

#>    Insured Coverage (z):
#>      Coefficient: -0.22 
#>      95% CI: [ -0.433 , -0.006 ]
#>      p-value: 0.0438

# Valuation Subsidy
subsidy_coef <- main_results %>% filter(term == "valuation_subsidy_z")
if (nrow(subsidy_coef) > 0) {
  cat("   Valuation Subsidy (z):\n")
  cat("     Coefficient:", round(subsidy_coef$estimate, 3), "\n")
  cat("     95% CI: [", round(subsidy_coef$conf.low, 3), ",", 
      round(subsidy_coef$conf.high, 3), "]\n")
  cat("     p-value:", format.pval(subsidy_coef$p.value, digits = 3), "\n\n")
}

#>    Valuation Subsidy (z):
#>      Coefficient: 0.104 
#>      95% CI: [ -0.011 , 0.219 ]
#>      p-value: 0.0761

# Interaction
interact_coef <- main_results %>% filter(term == "insured_coverage_z:valuation_subsidy_z")
if (nrow(interact_coef) > 0) {
  cat("   Coverage × Subsidy Interaction:\n")
  cat("     Coefficient:", round(interact_coef$estimate, 3), "\n")
  cat("     95% CI: [", round(interact_coef$conf.low, 3), ",", 
      round(interact_coef$conf.high, 3), "]\n")
  cat("     p-value:", format.pval(interact_coef$p.value, digits = 3), "\n\n")
}

#>    Coverage × Subsidy Interaction:
#>      Coefficient: 0.038 
#>      95% CI: [ -0.06 , 0.137 ]
#>      p-value: 0.445

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

#> ================================================================

Economic Interpretation

Based on the results:

Hypothesis 1: Insolvency × Subsidy Interaction

If the interaction term is negative and significant: Banks with LOWER insured coverage (more insolvent) AND HIGHER valuation subsidy were more likely to use BTFP.
This supports the hypothesis that BTFP’s par valuation attracted distressed banks.

Hypothesis 2: Run Risk

The Run Risk coefficient captures liquidity pressure independent of insolvency.
If significant, it shows that banks facing potential runs (high uninsured deposits × high MTM losses) sought Fed facilities.

Hypothesis 3: BTFP vs DW

Compare coefficients across facilities.
If valuation subsidy matters more for BTFP than DW, it confirms that BTFP’s design feature drove selection.

Save Results

################################################################################
# SAVE ALL OUTPUTS
################################################################################

cat("\n=== SAVING OUTPUTS ===\n\n")

#> 
#> === SAVING OUTPUTS ===

# Save analysis dataset
write_csv(analysis_clean, file.path(TABLE_PATH, "analysis_dataset.csv"))
cat("✓ Analysis dataset saved\n")

#> ✓ Analysis dataset saved

# Save regression tables
etable(m1_btfp, m2_btfp, m3_btfp, m4_btfp, m5_btfp, m6_btfp,
       tex = TRUE,
       file = file.path(TABLE_PATH, "table5_btfp_progressive.tex"),
       replace = TRUE)

etable(m5_btfp, m6_btfp, m5_dw, m6_dw,
       tex = TRUE,
       file = file.path(TABLE_PATH, "table6_btfp_vs_dw.tex"),
       replace = TRUE)

cat("✓ LaTeX tables saved\n")

#> ✓ LaTeX tables saved

# Save variable definitions
var_defs <- tibble::tribble(
  ~Variable, ~Definition, ~Construction,
  "btfp_user", "BTFP usage indicator", "=1 if any BTFP loan during acute crisis",
  "dw_user", "DW usage indicator", "=1 if any DW loan during acute crisis",
  "pct_uninsured_z", "Uninsured deposit share (z-score)", "Uninsured / Total Deposits, standardized",
  "pct_mtm_loss_z", "MTM loss ratio (z-score)", "MTM Loss / Total Assets, standardized",
  "run_risk_z", "Run risk (z-score)", "pct_uninsured × pct_mtm_loss, standardized",
  "insured_coverage_z", "Insured coverage ratio (z-score)", "(MV Assets - Uninsured - Insured) / Insured, standardized",
  "valuation_subsidy_z", "Valuation subsidy (z-score)", "MTM Loss on OMO-Eligible / Assets, standardized",
  "omo_eligible_z", "OMO-eligible ratio (z-score)", "OMO-Eligible / Assets, standardized",
  "mtm_insolvent", "MTM insolvency indicator", "=1 if MV Equity < 0"
)

write_csv(var_defs, file.path(TABLE_PATH, "variable_definitions.csv"))
cat("✓ Variable definitions saved\n")

#> ✓ Variable definitions saved

cat("\n")

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

#> ================================================================

cat("ANALYSIS COMPLETE\n")

#> ANALYSIS COMPLETE

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

#> ================================================================

cat("\nOutput locations:\n")

#> 
#> Output locations:

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

#>   Tables: C:/Users/mferdo2/OneDrive - Louisiana State University/Finance_PhD/DW_Stigma_paper/Liquidity_project_2025/03_documentation/regression_tables/btfp_analysis

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

#>   Figures: C:/Users/mferdo2/OneDrive - Louisiana State University/Finance_PhD/DW_Stigma_paper/Liquidity_project_2025/03_documentation/figures/btfp_analysis

Appendix: Model Diagnostics

################################################################################
# APPENDIX: MODEL DIAGNOSTICS
################################################################################

cat("\n=== MODEL DIAGNOSTICS ===\n\n")

#> 
#> === MODEL DIAGNOSTICS ===

# Main model summary
cat("Main Model (m6_btfp) Summary:\n")

#> Main Model (m6_btfp) Summary:

summary(m6_btfp)

#> GLM estimation, family = binomial, Dep. Var.: btfp_user
#> Observations: 4,270
#> Standard-errors: Heteroskedasticity-robust 
#>                                         Estimate Std. Error    z value
#> (Intercept)                            -2.575430   0.077327 -33.305804
#> run_risk_z                             -0.018858   0.112235  -0.168025
#> insured_coverage_z                     -0.219549   0.108894  -2.016178
#> valuation_subsidy_z                     0.104081   0.058668   1.774052
#> fhlb_ratio_z                            0.208028   0.053946   3.856257
#> wholesale_ratio_z                       0.144126   0.046220   3.118244
#> cash_ratio_z                           -0.545202   0.146319  -3.726109
#> book_equity_ratio_z                    -0.334893   0.073543  -4.553680
#> log_assets_z                            0.414115   0.055581   7.450655
#> insured_coverage_z:valuation_subsidy_z  0.038291   0.050153   0.763481
#>                                                    Pr(>|z|)    
#> (Intercept)                                       < 2.2e-16 ***
#> run_risk_z                             0.866563222134021593    
#> insured_coverage_z                     0.043781318508669557 *  
#> valuation_subsidy_z                    0.076054612875532376 .  
#> fhlb_ratio_z                           0.000115136628832644 ***
#> wholesale_ratio_z                      0.001819323743109357 ** 
#> cash_ratio_z                           0.000194458700604767 ***
#> book_equity_ratio_z                    0.000005271548598430 ***
#> log_assets_z                           0.000000000000092878 ***
#> insured_coverage_z:valuation_subsidy_z 0.445176442546601547    
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> Log-Likelihood: -1,293.4   Adj. Pseudo R2: 0.10358 
#>            BIC:  2,670.4     Squared Cor.: 0.076211

# Pseudo R-squared
cat("\nModel Fit Statistics:\n")

#> 
#> Model Fit Statistics:

cat("  Pseudo R² (McFadden):", round(fitstat(m6_btfp, "pr2")[[1]], 4), "\n")

#>   Pseudo R² (McFadden): 0.1098

# Predicted probabilities - use newdata to handle NAs properly
pred_prob <- predict(m6_btfp, newdata = analysis_omo, type = "response")

cat("\nPredicted Probability Distribution:\n")

#> 
#> Predicted Probability Distribution:

cat("  N predictions:", sum(!is.na(pred_prob)), "\n")

#>   N predictions: 4270

cat("  N missing:", sum(is.na(pred_prob)), "\n")

#>   N missing: 5

cat("  Mean:", round(mean(pred_prob, na.rm = TRUE), 4), "\n")

#>   Mean: 0.107

cat("  SD:", round(sd(pred_prob, na.rm = TRUE), 4), "\n")

#>   SD: 0.0883

cat("  Min:", round(min(pred_prob, na.rm = TRUE), 4), "\n")

#>   Min: 0.0012

cat("  Max:", round(max(pred_prob, na.rm = TRUE), 4), "\n")

#>   Max: 0.6867

# Confusion matrix at 0.5 threshold - only use complete cases
pred_class <- ifelse(pred_prob > 0.5, 1, 0)

# Create confusion matrix excluding NAs
valid_idx <- !is.na(pred_class) & !is.na(analysis_omo$btfp_user)
conf_matrix <- table(
  Predicted = pred_class[valid_idx], 
  Actual = analysis_omo$btfp_user[valid_idx]
)

cat("\nConfusion Matrix (threshold = 0.5):\n")

#> 
#> Confusion Matrix (threshold = 0.5):

print(conf_matrix)

#>          Actual
#> Predicted    0    1
#>         0 3802  450
#>         1   11    7

accuracy <- sum(diag(conf_matrix)) / sum(conf_matrix)
cat("\nAccuracy:", round(accuracy * 100, 2), "%\n")

#> 
#> Accuracy: 89.2 %

cat("Note: Based on", sum(valid_idx), "observations with complete data\n")

#> Note: Based on 4270 observations with complete data

Empirically Estimating the Run Rate (s)

Motivation

The insured deposit coverage ratio depends on the assumed run rate (\(s\)) on uninsured deposits:

\[\text{Insured Coverage Ratio} = \frac{\text{MV Assets} - s \times \text{Uninsured Deposits} - \text{Insured Deposits}}{\text{Insured Deposits}}\]

Rather than assuming arbitrary values for \(s\), we estimate it empirically using actual deposit outflows observed during the 2023 banking crisis.

Approach: Actual Deposit Outflows (2022Q4 to 2023Q1)

We calculate the realized run rate as the percentage decline in uninsured deposits from 2022Q4 (pre-crisis) to 2023Q1 (crisis quarter):

\[s_{i} = \frac{\text{Uninsured Deposits}_{2022Q4} - \text{Uninsured Deposits}_{2023Q1}}{\text{Uninsured Deposits}_{2022Q4}}\]

################################################################################
# EMPIRICALLY ESTIMATE RUN RATE (s) FROM ACTUAL DEPOSIT OUTFLOWS
################################################################################

cat("\n=== ESTIMATING RUN RATE FROM ACTUAL DEPOSIT FLOWS ===\n\n")

#> 
#> === ESTIMATING RUN RATE FROM ACTUAL DEPOSIT FLOWS ===

# Calculate actual deposit changes from 2022Q4 to 2023Q1
deposit_changes <- call_q %>%
  filter(quarter %in% c("2022Q4", "2023Q1")) %>%
  select(idrssd, quarter, uninsured_deposit, insured_deposit, total_deposit) %>%
  pivot_wider(
    id_cols = idrssd,
    names_from = quarter,
    values_from = c(uninsured_deposit, insured_deposit, total_deposit)
  ) %>%
  filter(
    !is.na(uninsured_deposit_2022Q4),
    !is.na(uninsured_deposit_2023Q1),
    uninsured_deposit_2022Q4 > 0  # Need positive baseline to calculate rate
  ) %>%
  mutate(
    # Raw change in uninsured deposits
    uninsured_change = uninsured_deposit_2022Q4 - uninsured_deposit_2023Q1,
    
    # Actual run rate on uninsured deposits
    # Positive = outflow (run), Negative = inflow
    s_actual_raw = uninsured_change / uninsured_deposit_2022Q4,
    
    # Bounded run rate [0, 1] for economic interpretation
    # s = 0 means no run, s = 1 means 100% of uninsured deposits left
    s_actual = pmax(0, pmin(1, s_actual_raw)),
    
    # Also calculate insured deposit change (for comparison)
    insured_change = insured_deposit_2022Q4 - insured_deposit_2023Q1,
    s_insured = pmax(0, pmin(1, insured_change / pmax(insured_deposit_2022Q4, 1))),
    
    # Total deposit change
    total_change = total_deposit_2022Q4 - total_deposit_2023Q1,
    s_total = total_change / pmax(total_deposit_2022Q4, 1)
  )

cat("Sample size for run rate estimation:", nrow(deposit_changes), "banks\n\n")

#> Sample size for run rate estimation: 4640 banks

Distribution of Actual Run Rates

################################################################################
# SUMMARY STATISTICS FOR ACTUAL RUN RATES
################################################################################

# Summary statistics
cat("=== Summary Statistics: Actual Uninsured Deposit Run Rate (s) ===\n\n")

#> === Summary Statistics: Actual Uninsured Deposit Run Rate (s) ===

s_summary <- deposit_changes %>%
  summarise(
    N = n(),
    Mean = mean(s_actual, na.rm = TRUE),
    SD = sd(s_actual, na.rm = TRUE),
    Min = min(s_actual, na.rm = TRUE),
    P10 = quantile(s_actual, 0.10, na.rm = TRUE),
    P25 = quantile(s_actual, 0.25, na.rm = TRUE),
    Median = quantile(s_actual, 0.50, na.rm = TRUE),
    P75 = quantile(s_actual, 0.75, na.rm = TRUE),
    P90 = quantile(s_actual, 0.90, na.rm = TRUE),
    Max = max(s_actual, na.rm = TRUE),
    Pct_Positive_Outflow = 100 * mean(s_actual > 0, na.rm = TRUE),
    Pct_Large_Run = 100 * mean(s_actual > 0.10, na.rm = TRUE)
  )

# Display as table
s_summary_display <- s_summary %>%
  pivot_longer(everything(), names_to = "Statistic", values_to = "Value") %>%
  mutate(Value = round(Value, 4))

kable(s_summary_display, 
      caption = "Table: Distribution of Actual Uninsured Deposit Run Rates (2022Q4 to 2023Q1)",
      col.names = c("Statistic", "Value")) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
                full_width = FALSE)

Table: Distribution of Actual Uninsured Deposit Run Rates (2022Q4 to 2023Q1)
Statistic	Value
N	4640.0000
Mean	0.0815
SD	0.1076
Min	0.0000
P10	0.0000
P25	0.0000
Median	0.0445
P75	0.1257
P90	0.2176
Max	0.9954
Pct_Positive_Outflow	64.3103
Pct_Large_Run	32.3922

################################################################################
# COMPARE RUN RATES: UNINSURED VS INSURED DEPOSITS
################################################################################

cat("\n=== Comparison: Uninsured vs Insured Deposit Outflows ===\n\n")

#> 
#> === Comparison: Uninsured vs Insured Deposit Outflows ===

comparison_summary <- deposit_changes %>%
  summarise(
    `Uninsured Mean Run Rate` = mean(s_actual, na.rm = TRUE),
    `Uninsured Median Run Rate` = median(s_actual, na.rm = TRUE),
    `Insured Mean Run Rate` = mean(s_insured, na.rm = TRUE),
    `Insured Median Run Rate` = median(s_insured, na.rm = TRUE),
    `% Banks with Uninsured Outflow` = 100 * mean(s_actual > 0, na.rm = TRUE),
    `% Banks with Insured Outflow` = 100 * mean(s_insured > 0, na.rm = TRUE)
  ) %>%
  pivot_longer(everything(), names_to = "Metric", values_to = "Value") %>%
  mutate(Value = round(Value, 4))

kable(comparison_summary,
      caption = "Table: Uninsured vs Insured Deposit Outflows",
      col.names = c("Metric", "Value")) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
                full_width = FALSE)

Table: Uninsured vs Insured Deposit Outflows
Metric	Value
Uninsured Mean Run Rate	0.0815
Uninsured Median Run Rate	0.0445
Insured Mean Run Rate	0.0087
Insured Median Run Rate	0.0000
% Banks with Uninsured Outflow	64.3103
% Banks with Insured Outflow	32.2845

Visualization: Distribution of Run Rates

################################################################################
# VISUALIZE RUN RATE DISTRIBUTION
################################################################################

# Histogram of actual run rates
p_hist <- ggplot(deposit_changes, aes(x = s_actual)) +
  geom_histogram(bins = 50, fill = "#2166ac", color = "white", alpha = 0.8) +
  geom_vline(aes(xintercept = mean(s_actual, na.rm = TRUE)), 
             color = "red", linetype = "dashed", linewidth = 1) +
  geom_vline(aes(xintercept = median(s_actual, na.rm = TRUE)), 
             color = "darkgreen", linetype = "dashed", linewidth = 1) +
  annotate("text", x = mean(deposit_changes$s_actual, na.rm = TRUE) + 0.05, y = Inf, 
           label = paste("Mean =", round(mean(deposit_changes$s_actual, na.rm = TRUE), 3)),
           vjust = 2, color = "red", size = 4) +
  annotate("text", x = median(deposit_changes$s_actual, na.rm = TRUE) + 0.05, y = Inf,
           label = paste("Median =", round(median(deposit_changes$s_actual, na.rm = TRUE), 3)),
           vjust = 4, color = "darkgreen", size = 4) +
  labs(
    title = "Distribution of Actual Uninsured Deposit Run Rates",
    subtitle = "2022Q4 to 2023Q1 (Positive = Outflow)",
    x = "Run Rate (s)",
    y = "Number of Banks"
  ) +
  scale_x_continuous(labels = scales::percent_format(), limits = c(0, 1)) +
  theme_minimal(base_size = 12)

print(p_hist)

# Save figure
ggsave(file.path(FIG_PATH, "run_rate_distribution.png"), p_hist, 
       width = 10, height = 5, dpi = 300)

Run Rates by Bank Characteristics

################################################################################
# RUN RATES BY BANK CHARACTERISTICS
################################################################################

# Merge run rates with baseline characteristics
run_rate_analysis <- deposit_changes %>%
  select(idrssd, s_actual, s_actual_raw, uninsured_change) %>%
  left_join(baseline, by = "idrssd") %>%
  filter(!is.na(s_actual), !is.na(mtm_insolvent))

cat("\n=== Run Rates by Bank Characteristics ===\n\n")

#> 
#> === Run Rates by Bank Characteristics ===

# By MTM insolvency status
run_by_insolvency <- run_rate_analysis %>%
  group_by(mtm_insolvent) %>%
  summarise(
    N = n(),
    Mean_Run_Rate = mean(s_actual, na.rm = TRUE),
    Median_Run_Rate = median(s_actual, na.rm = TRUE),
    SD_Run_Rate = sd(s_actual, na.rm = TRUE),
    Pct_Any_Outflow = 100 * mean(s_actual > 0, na.rm = TRUE),
    Pct_Large_Run = 100 * mean(s_actual > 0.10, na.rm = TRUE),
    .groups = "drop"
  ) %>%
  mutate(mtm_insolvent = ifelse(mtm_insolvent, "MTM Insolvent", "MTM Solvent"))

kable(run_by_insolvency,
      caption = "Table: Actual Run Rates by MTM Insolvency Status",
      digits = 4) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
                full_width = FALSE)

Table: Actual Run Rates by MTM Insolvency Status
mtm_insolvent	N	Mean_Run_Rate	Median_Run_Rate	SD_Run_Rate	Pct_Any_Outflow	Pct_Large_Run
MTM Solvent	3764	0.0819	0.0433	0.1097	63.4166	32.3858
MTM Insolvent	861	0.0812	0.0496	0.0983	68.9895	32.8688

# By size bin
if ("size_bin" %in% names(run_rate_analysis)) {
  run_by_size <- run_rate_analysis %>%
    filter(!is.na(size_bin)) %>%
    group_by(size_bin) %>%
    summarise(
      N = n(),
      Mean_Run_Rate = mean(s_actual, na.rm = TRUE),
      Median_Run_Rate = median(s_actual, na.rm = TRUE),
      Pct_Large_Run = 100 * mean(s_actual > 0.10, na.rm = TRUE),
      .groups = "drop"
    )
  
  kable(run_by_size,
        caption = "Table: Actual Run Rates by Bank Size",
        digits = 4) %>%
    kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
                  full_width = FALSE)
}

Table: Actual Run Rates by Bank Size
size_bin	N	Mean_Run_Rate	Median_Run_Rate	Pct_Large_Run
gsib	29	0.0881	0.0424	31.0345
large	724	0.0793	0.0616	32.3204
small	3872	0.0822	0.0407	32.5155

# By BTFP usage
run_rate_with_btfp <- run_rate_analysis %>%
  left_join(btfp_crisis %>% select(idrssd, btfp_user), by = "idrssd") %>%
  mutate(btfp_user = coalesce(btfp_user, 0L))

run_by_btfp <- run_rate_with_btfp %>%
  group_by(btfp_user) %>%
  summarise(
    N = n(),
    Mean_Run_Rate = mean(s_actual, na.rm = TRUE),
    Median_Run_Rate = median(s_actual, na.rm = TRUE),
    SD_Run_Rate = sd(s_actual, na.rm = TRUE),
    Pct_Large_Run = 100 * mean(s_actual > 0.10, na.rm = TRUE),
    .groups = "drop"
  ) %>%
  mutate(btfp_user = ifelse(btfp_user == 1, "BTFP User", "Non-User"))

kable(run_by_btfp,
      caption = "Table: Actual Run Rates by BTFP Usage",
      digits = 4) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
                full_width = FALSE)

Table: Actual Run Rates by BTFP Usage
btfp_user	N	Mean_Run_Rate	Median_Run_Rate	SD_Run_Rate	Pct_Large_Run
Non-User	4152	0.0804	0.0415	0.1077	31.8642
BTFP User	473	0.0933	0.0724	0.1064	37.8436

Key Finding: Empirically Estimated Run Rate

################################################################################
# KEY FINDING: EMPIRICALLY ESTIMATED RUN RATE
################################################################################

# Calculate key statistics
mean_s <- mean(deposit_changes$s_actual, na.rm = TRUE)
median_s <- median(deposit_changes$s_actual, na.rm = TRUE)
p75_s <- quantile(deposit_changes$s_actual, 0.75, na.rm = TRUE)
p90_s <- quantile(deposit_changes$s_actual, 0.90, na.rm = TRUE)

cat("\n")

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

#> ================================================================

cat("KEY FINDING: EMPIRICALLY ESTIMATED RUN RATE (s)\n")

#> KEY FINDING: EMPIRICALLY ESTIMATED RUN RATE (s)

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

#> ================================================================

cat("Based on actual deposit outflows from 2022Q4 to 2023Q1:\n\n")

#> Based on actual deposit outflows from 2022Q4 to 2023Q1:

cat("  Mean run rate:     s =", round(mean_s, 4), "(", round(mean_s * 100, 1), "% )\n")

#>   Mean run rate:     s = 0.0815 ( 8.2 % )

cat("  Median run rate:   s =", round(median_s, 4), "(", round(median_s * 100, 1), "% )\n")

#>   Median run rate:   s = 0.0445 ( 4.4 % )

cat("  75th percentile:   s =", round(p75_s, 4), "(", round(p75_s * 100, 1), "% )\n")

#>   75th percentile:   s = 0.1257 ( 12.6 % )

cat("  90th percentile:   s =", round(p90_s, 4), "(", round(p90_s * 100, 1), "% )\n\n")

#>   90th percentile:   s = 0.2176 ( 21.8 % )

cat("Recommendation: Use s =", round(median_s, 2), "as baseline,\n")

#> Recommendation: Use s = 0.04 as baseline,

cat("                with s =", round(p75_s, 2), "and s =", round(p90_s, 2), "for robustness.\n")

#>                 with s = 0.13 and s = 0.22 for robustness.

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

#> ================================================================

Using Empirically Estimated s in Coverage Ratio

################################################################################
# CALCULATE COVERAGE RATIO WITH EMPIRICALLY ESTIMATED s
################################################################################

# Use median run rate as primary estimate
s_empirical <- median(deposit_changes$s_actual, na.rm = TRUE)

cat("\n=== Calculating Coverage Ratio with Empirical s =", round(s_empirical, 3), "===\n\n")

#> 
#> === Calculating Coverage Ratio with Empirical s = 0.044 ===

# Add to analysis dataset
analysis_omo <- analysis_omo %>%
  mutate(
    # Coverage with empirically estimated s
    insured_coverage_empirical = (mv_asset - s_empirical * uninsured_deposit - insured_deposit) / 
                                  pmax(insured_deposit, 1),
    insured_coverage_empirical_z = standardize(winsorize(insured_coverage_empirical))
  )

# Run model with empirically estimated s
m_empirical_s <- feglm(
  btfp_user ~ insured_coverage_empirical_z + run_risk_z + valuation_subsidy_z + log_assets_z,
  data = analysis_omo,
  family = binomial("logit"),
  vcov = "hetero"
)

cat("Model with empirically estimated s:\n")

#> Model with empirically estimated s:

summary(m_empirical_s)

#> GLM estimation, family = binomial, Dep. Var.: btfp_user
#> Observations: 4,275
#> Standard-errors: Heteroskedasticity-robust 
#>                               Estimate Std. Error   z value  Pr(>|z|)    
#> (Intercept)                  -2.340180   0.055740 -41.98373 < 2.2e-16 ***
#> insured_coverage_empirical_z -0.214565   0.075625  -2.83720  0.004551 ** 
#> run_risk_z                    0.182848   0.080547   2.27007  0.023203 *  
#> valuation_subsidy_z           0.119138   0.048330   2.46508  0.013698 *  
#> log_assets_z                  0.438925   0.048001   9.14404 < 2.2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> Log-Likelihood: -1,346.5   Adj. Pseudo R2: 0.072169
#>            BIC:  2,734.8     Squared Cor.: 0.050105

Robustness: Models Across Different s Values

################################################################################
# ROBUSTNESS: COMPARE MODELS ACROSS DIFFERENT s VALUES
################################################################################

cat("\n=== Robustness: Coverage Ratio Under Different Run Rates ===\n\n")

#> 
#> === Robustness: Coverage Ratio Under Different Run Rates ===

# Calculate coverage for multiple s values
analysis_omo <- analysis_omo %>%
  mutate(
    # s = 0% (no run)
    coverage_s0 = (mv_asset - 0.00 * uninsured_deposit - insured_deposit) / pmax(insured_deposit, 1),
    coverage_s0_z = standardize(winsorize(coverage_s0)),
    
    # s = 10%
    coverage_s10 = (mv_asset - 0.10 * uninsured_deposit - insured_deposit) / pmax(insured_deposit, 1),
    coverage_s10_z = standardize(winsorize(coverage_s10)),
    
    # s = 25%
    coverage_s25 = (mv_asset - 0.25 * uninsured_deposit - insured_deposit) / pmax(insured_deposit, 1),
    coverage_s25_z = standardize(winsorize(coverage_s25)),
    
    # s = 50%
    coverage_s50 = (mv_asset - 0.50 * uninsured_deposit - insured_deposit) / pmax(insured_deposit, 1),
    coverage_s50_z = standardize(winsorize(coverage_s50)),
    
    # s = 75%
    coverage_s75 = (mv_asset - 0.75 * uninsured_deposit - insured_deposit) / pmax(insured_deposit, 1),
    coverage_s75_z = standardize(winsorize(coverage_s75)),
    
    # s = 100% (full run)
    coverage_s100 = (mv_asset - 1.00 * uninsured_deposit - insured_deposit) / pmax(insured_deposit, 1),
    coverage_s100_z = standardize(winsorize(coverage_s100))
  )

# Run models for each s value
m_s0 <- feglm(btfp_user ~ coverage_s0_z + run_risk_z + log_assets_z,
              data = analysis_omo, family = binomial("logit"), vcov = "hetero")
m_s10 <- feglm(btfp_user ~ coverage_s10_z + run_risk_z + log_assets_z,
               data = analysis_omo, family = binomial("logit"), vcov = "hetero")
m_s25 <- feglm(btfp_user ~ coverage_s25_z + run_risk_z + log_assets_z,
               data = analysis_omo, family = binomial("logit"), vcov = "hetero")
m_s50 <- feglm(btfp_user ~ coverage_s50_z + run_risk_z + log_assets_z,
               data = analysis_omo, family = binomial("logit"), vcov = "hetero")
m_s75 <- feglm(btfp_user ~ coverage_s75_z + run_risk_z + log_assets_z,
               data = analysis_omo, family = binomial("logit"), vcov = "hetero")
m_s100 <- feglm(btfp_user ~ coverage_s100_z + run_risk_z + log_assets_z,
                data = analysis_omo, family = binomial("logit"), vcov = "hetero")

# Display comparison table
etable(m_s0, m_s10, m_s25, m_s50, m_s75, m_s100,
       headers = c("s=0%", "s=10%", "s=25%", "s=50%", "s=75%", "s=100%"),
       title = "Table: BTFP Usage Under Different Run Rate Assumptions",
       notes = "Coverage ratio calculated as (MV Assets - s × Uninsured - Insured) / Insured",
       fitstat = ~ n + pr2)

#>                               m_s0              m_s10              m_s25
#>                               s=0%              s=10%              s=25%
#> Dependent Var.:          btfp_user          btfp_user          btfp_user
#>                                                                         
#> Constant        -2.330*** (0.0556) -2.329*** (0.0556) -2.329*** (0.0556)
#> coverage_s0_z   -0.2394** (0.0740)                                      
#> Run Risk (z)      0.1810* (0.0800)   0.1812* (0.0805)   0.1824* (0.0813)
#> Log Assets (z)  0.4563*** (0.0472) 0.4555*** (0.0472) 0.4544*** (0.0471)
#> coverage_s10_z                     -0.2382** (0.0745)                   
#> coverage_s25_z                                        -0.2355** (0.0751)
#> coverage_s50_z                                                          
#> coverage_s75_z                                                          
#> coverage_s100_z                                                         
#> _______________ __________________ __________________ __________________
#> S.E. type       Heteroskedas.-rob. Heteroskedas.-rob. Heteroskedas.-rob.
#> Observations                 4,275              4,275              4,275
#> Pseudo R2                  0.07294            0.07287            0.07275
#> 
#>                              m_s50              m_s75             m_s100
#>                              s=50%              s=75%             s=100%
#> Dependent Var.:          btfp_user          btfp_user          btfp_user
#>                                                                         
#> Constant        -2.328*** (0.0556) -2.328*** (0.0556) -2.327*** (0.0556)
#> coverage_s0_z                                                           
#> Run Risk (z)      0.1857* (0.0824)   0.1906* (0.0835)   0.1963* (0.0845)
#> Log Assets (z)  0.4524*** (0.0470) 0.4503*** (0.0469) 0.4481*** (0.0467)
#> coverage_s10_z                                                          
#> coverage_s25_z                                                          
#> coverage_s50_z  -0.2294** (0.0760)                                      
#> coverage_s75_z                     -0.2220** (0.0768)                   
#> coverage_s100_z                                       -0.2140** (0.0779)
#> _______________ __________________ __________________ __________________
#> S.E. type       Heteroskedas.-rob. Heteroskedas.-rob. Heteroskedas.-rob.
#> Observations                 4,275              4,275              4,275
#> Pseudo R2                  0.07254            0.07232            0.07208
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Visualization: Coverage Coefficient Across s Values

################################################################################
# VISUALIZE SENSITIVITY OF COVERAGE COEFFICIENT TO s
################################################################################

# Extract coefficients
s_sensitivity <- data.frame(
  s = c(0, 0.10, 0.25, 0.50, 0.75, 1.00),
  coef = c(
    coef(m_s0)["coverage_s0_z"],
    coef(m_s10)["coverage_s10_z"],
    coef(m_s25)["coverage_s25_z"],
    coef(m_s50)["coverage_s50_z"],
    coef(m_s75)["coverage_s75_z"],
    coef(m_s100)["coverage_s100_z"]
  ),
  se = c(
    se(m_s0)["coverage_s0_z"],
    se(m_s10)["coverage_s10_z"],
    se(m_s25)["coverage_s25_z"],
    se(m_s50)["coverage_s50_z"],
    se(m_s75)["coverage_s75_z"],
    se(m_s100)["coverage_s100_z"]
  )
) %>%
  mutate(
    ci_low = coef - 1.96 * se,
    ci_high = coef + 1.96 * se
  )

# Plot
p_sensitivity <- ggplot(s_sensitivity, aes(x = s, y = coef)) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray50") +
  geom_ribbon(aes(ymin = ci_low, ymax = ci_high), alpha = 0.2, fill = "#2166ac") +
  geom_line(color = "#2166ac", linewidth = 1) +
  geom_point(color = "#2166ac", size = 3) +
  geom_vline(xintercept = s_empirical, linetype = "dashed", color = "red", linewidth = 0.8) +
  annotate("text", x = s_empirical + 0.05, y = max(s_sensitivity$coef), 
           label = paste("Empirical s =", round(s_empirical, 2)),
           hjust = 0, color = "red", size = 4) +
  labs(
    title = "Sensitivity of Coverage Ratio Coefficient to Assumed Run Rate (s)",
    subtitle = "Coefficient on Insured Coverage Ratio in BTFP Usage Model",
    x = "Assumed Run Rate on Uninsured Deposits (s)",
    y = "Coefficient (with 95% CI)"
  ) +
  scale_x_continuous(labels = scales::percent_format(), breaks = seq(0, 1, 0.25)) +
  theme_minimal(base_size = 12)

print(p_sensitivity)

# Save figure
ggsave(file.path(FIG_PATH, "coverage_sensitivity_to_s.png"), p_sensitivity,
       width = 10, height = 5, dpi = 300)

Summary: Empirically Estimated Run Rate

Based on actual deposit outflows observed during the 2023 banking crisis:

Statistic	Value	Interpretation
Mean s	0.082	Average bank experienced 8.2% uninsured deposit outflow
Median s	0.044	Typical bank experienced 4.4% outflow
75th percentile	0.126	Stressed banks experienced up to 12.6% outflow
90th percentile	0.218	Severely stressed banks saw 21.8% outflow

Key Finding: The empirically estimated run rate is substantially lower than the 100% assumption often used in stress tests, but higher than zero. Using the median empirical value of s = 0.04 provides a data-driven approach to calculating the insured deposit coverage ratio.

Robustness: Results are qualitatively similar across different assumed values of s, with the coverage ratio coefficient remaining negative (indicating that lower coverage/higher insolvency risk increases BTFP usage) regardless of the specific run rate assumption.

Session Info

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] moments_0.14.1    DescTools_0.99.60 kableExtra_1.4.0  knitr_1.50       
#>  [5] patchwork_1.3.0   scales_1.3.0      ggthemes_5.1.0    ggplot2_3.5.1    
#>  [9] broom_1.0.8       lmtest_0.9-40     zoo_1.8-12        sandwich_3.1-1   
#> [13] fixest_0.12.1     readr_2.1.5       stringr_1.5.1     lubridate_1.9.4  
#> [17] tibble_3.2.1      tidyr_1.3.1       dplyr_1.1.4       data.table_1.17.0
#> 
#> loaded via a namespace (and not attached):
#>  [1] tidyselect_1.2.1    Exact_3.3           viridisLite_0.4.2  
#>  [4] rootSolve_1.8.2.4   farver_2.1.1        fastmap_1.1.1      
#>  [7] digest_0.6.33       timechange_0.3.0    lifecycle_1.0.4    
#> [10] dreamerr_1.4.0      lmom_3.2            magrittr_2.0.3     
#> [13] compiler_4.3.1      rlang_1.1.2         sass_0.4.9         
#> [16] tools_4.3.1         utf8_1.2.4          yaml_2.3.7         
#> [19] labeling_0.4.3      bit_4.0.5           xml2_1.3.6         
#> [22] expm_1.0-0          withr_2.5.2         purrr_1.0.2        
#> [25] numDeriv_2016.8-1.1 grid_4.3.1          fansi_1.0.5        
#> [28] e1071_1.7-14        colorspace_2.1-0    MASS_7.3-60        
#> [31] cli_3.6.1           mvtnorm_1.3-3       crayon_1.5.2       
#> [34] rmarkdown_2.29      ragg_1.3.0          generics_0.1.3     
#> [37] rstudioapi_0.16.0   httr_1.4.7          tzdb_0.4.0         
#> [40] readxl_1.4.3        gld_2.6.7           cachem_1.0.8       
#> [43] proxy_0.4-27        parallel_4.3.1      cellranger_1.1.0   
#> [46] stringmagic_1.1.2   vctrs_0.6.4         boot_1.3-28.1      
#> [49] Matrix_1.5-4.1      jsonlite_1.8.7      hms_1.1.3          
#> [52] bit64_4.0.5         Formula_1.2-5       systemfonts_1.0.6  
#> [55] jquerylib_0.1.4     glue_1.6.2          stringi_1.7.12     
#> [58] gtable_0.3.6        munsell_0.5.0       pillar_1.9.0       
#> [61] htmltools_0.5.7     R6_2.5.1            textshaping_0.3.7  
#> [64] vroom_1.6.5         evaluate_0.23       lattice_0.21-8     
#> [67] haven_2.5.4         backports_1.4.1     bslib_0.5.1        
#> [70] class_7.3-22        Rcpp_1.0.12         svglite_2.1.3      
#> [73] nlme_3.1-162        xfun_0.54           fs_1.6.3           
#> [76] forcats_1.0.0       pkgconfig_2.0.3

95% Winzorized

Who Used BTFP and Why? A Causal Analysis

Mohua Ferdousi

2025-12-10