1 Summary

This analysis examines bank borrowing choices during the 2023 banking crisis using three complementary approaches:

Model Framework:

  1. Linear Probability Model (LPM): BTFP vs DW among borrowers
  2. Logit Model: BTFP vs DW among borrowers
  3. Multinomial Logit: Choice among {No participation, BTFP only, DW only, Both BTFP+DW}
  4. Trivariate Probit: Joint modeling of BTFP, DW, and FHLB decisions with error correlations

Key Variables:

  • MTM Losses: Separated into BTFP-eligible (omo_eligible) vs. other assets
  • Uninsured Deposits: Run risk exposure
  • FHLB Borrowing: Crisis-period increases (2+ std above historical mean)
  • Run Risk Measures: Continuous and binary indicators
  • Insolvency Measures: Three alternative specifications

Sample:

  • Excludes: G-SIB banks, failed banks
  • Period: March 13 - April 30, 2023 (acute crisis phase)
  • Baseline: 2022Q4 characteristics

2 Setup and Configuration

################################################################################
# SECTION 0: SETUP AND PACKAGES
################################################################################

# Update rlang first if needed
if (packageVersion("rlang") < "1.1.4") {
  cat("Updating rlang to resolve dependencies...\n")
  install.packages("rlang", repos = "https://cloud.r-project.org")
}
#> Updating rlang to resolve dependencies...
#> 
#>   There is a binary version available but the source version is later:
#>       binary source needs_compilation
#> rlang  1.1.5  1.1.6              TRUE
# Required packages
required_packages <- c(
  # Data manipulation
  "data.table", "dplyr", "tidyr", "tibble", "lubridate", "stringr", "readr",
  
  # Modeling
  "fixest",           # Fast fixed effects
  "nnet",             # Multinomial logit
  "mlogit",           # Advanced multinomial models
  "mvtnorm",          # Multivariate normal
  "sandwich",         # Robust standard errors
  "lmtest",           # Model testing
  "broom",            # Tidy model output
  "margins",          # Marginal effects
  
  # Visualization
  "ggplot2", "ggthemes", "scales", "patchwork", "gridExtra",
  
  # Tables
  "knitr", "kableExtra", "stargazer", "modelsummary",
  
  # Statistics
  "DescTools", "moments", "psych"
)

# 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", 
                   dependencies = TRUE)
}

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

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

cat("✓ All packages loaded successfully\n")
#> ✓ All packages loaded successfully
################################################################################
# CONFIGURE PATHS
################################################################################

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/multinomial_analysis")
FIG_PATH    <- file.path(DOC_PATH, "figures/multinomial_analysis")

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

cat("Working directory:", getwd(), "\n")
#> Working directory: C:/Users/mferdo2/OneDrive - Louisiana State University/Finance_PhD/DW_Stigma_paper/Liquidity_project_2025
cat("Data paths configured successfully\n")
#> Data paths configured successfully
################################################################################
# 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
BASELINE_DATE     <- "2022Q4"               # 2022Q4

cat("\n")
cat("================================================================\n")
#> ================================================================
cat("MULTINOMIAL CHOICE ANALYSIS - 2023 BANKING CRISIS\n")
#> MULTINOMIAL CHOICE 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("Baseline Date:     ", BASELINE_DATE, "\n")
#> Baseline Date:      2022Q4
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])
}

#' Winsorize at 95th percentile
winsorize_95 <- function(x) {
  winsorize(x, probs = c(0.025, 0.975))
}

#' Standardize (z-score transformation)
standardize <- function(x) {
  (x - mean(x, na.rm = TRUE)) / sd(x, na.rm = TRUE)
}

#' Safe division
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"
    )
}

cat("✓ Helper functions loaded\n")
#> ✓ Helper functions loaded

3 Data Loading

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

cat("\n")
cat(strrep("=", 80), "\n", sep = "")
#> ================================================================================
cat("LOADING DATA\n")
#> LOADING DATA
cat(strrep("=", 80), "\n", sep = "")
#> ================================================================================
# 1.1 Load Call Report Data
call_q <- read_csv(
  file.path(DATA_PROC, "final_call_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
# 1.4 Load Public Bank Data (CRSP)
file_public_bank <- file.path(DATA_PROC, "df_public_bank.csv")
if (file.exists(file_public_bank)) {
  df_public <- read_csv(file_public_bank, show_col_types = FALSE) %>%
    mutate(
      rssd_bank = as.character(rssd_bank),
      DlyCalDt = as.Date(DlyCalDt)
    )
  cat("✓ Public bank data loaded:", nrow(df_public), "observations\n")
} else {
  cat("  WARNING: Public bank file not found\n")
  df_public <- NULL
}
#> ✓ Public bank data loaded: 135943 observations
cat("\n✓ All data loaded successfully\n")
#> 
#> ✓ All data loaded successfully

4 Baseline Sample Creation

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

cat("\n")
cat(strrep("=", 80), "\n", sep = "")
#> ================================================================================
cat("AGGREGATING FACILITY USAGE\n")
#> AGGREGATING FACILITY USAGE
cat(strrep("=", 80), "\n", sep = "")
#> ================================================================================
# 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
cat("\n✓ Facility usage aggregated\n")
#> 
#> ✓ Facility usage aggregated
################################################################################
# SECTION 3: CREATE 2022Q4 BASELINE WITH KEY VARIABLES
################################################################################

cat("\n")
cat(strrep("=", 80), "\n", sep = "")
#> ================================================================================
cat("CREATING 2022Q4 BASELINE\n")
#> CREATING 2022Q4 BASELINE
cat(strrep("=", 80), "\n", sep = "")
#> ================================================================================
# Check available columns
has_failed_bank <- "failed_bank" %in% names(call_q)
has_gsib <- "gsib" %in% names(call_q)
has_size_bin <- "size_bin" %in% names(call_q)
has_mv_asset <- "mm_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("  - gsib column:", ifelse(has_gsib, "YES", "NO"), "\n")
#>   - gsib column: YES
cat("  - size_bin column:", ifelse(has_size_bin, "YES", "NO"), "\n")
#>   - size_bin column: YES
cat("  - mv_asset column:", ifelse(has_mv_asset, "YES", "NO"), "\n\n")
#>   - mv_asset column: YES
# Create baseline dataset from 2022Q4 (preserving all Python-created variables)
baseline <- call_q %>%
  filter(quarter == BASELINE_DATE) %>%
  transmute(
    idrssd = as.character(idrssd),
    
    # Bank identifiers and flags
    failed_bank = if (has_failed_bank) failed_bank else 0L,
    gsib = if (has_gsib) gsib else 0L,
    
    # Balance sheet fundamentals
    size_bin = if (has_size_bin) size_bin else NA_character_,
    total_asset,
    total_liability,
    total_equity,
    total_deposit,
    insured_deposit = pmax(insured_deposit, 1),
    uninsured_deposit = pmax(uninsured_deposit, 0),
    security_to_total_asset,
    loan_to_deposit,
    roa,
    
    # Cash and liquid assets
    cash,
    fed_fund_sold,
    rerepo,
    
    # Securities (from Python cleaning)
    security,           # total_sec in some versions
    omo_eligible,       # Treasury + Agency debt + Agency MBS (BTFP-eligible)
    non_omo_securities, # Other securities
    
    # Loans
    total_loan,
    
    # Size measures
    log_assets = log(pmax(total_asset, 1)),
    
    # MTM variables (from Python MTM calculation)
    mm_asset = if (has_mv_asset) mm_asset else total_asset,
    mtm_total_loss,
    mtm_loss_to_total_asset,
    equity_after_mtm,
    
    # MTM losses: OMO-eligible vs Other
    mtm_loss_omo_eligible,
    mtm_loss_omo_eligible_to_omo_eligible = safe_div(mtm_loss_omo_eligible, omo_eligible),
    mtm_loss_omo_eligible_to_total_asset,
    omo_eligible_to_total_asset,
    
    # Calculate MTM loss on OTHER assets (non-OMO-eligible)
    
    mtm_loss_other = mtm_loss_non_omo_eligible,
    mtm_loss_other_to_total_asset = mtm_loss_non_omo_eligible_to_total_asset,
    
    # Wholesale funding & FHLB
    repo,
    fed_fund_purchase,
    other_borr,
    other_borrowed_less_than_1yr,
    fhlb_adv,
    fhlb_less_than_1yr,
    fhlb_to_total_asset,
    
    # Control variables
    cash_to_total_asset,
    book_equity_to_total_asset,
    
    # Additional controls
    net_income_q,
    avg_assets
  )

cat("✓ Initial baseline:", nrow(baseline), "banks\n")
#> ✓ Initial baseline: 4737 banks
# Apply sample filters
baseline <- baseline %>%
  filter(
    # Exclude G-SIBs
    gsib == 0 | is.na(gsib),
    # Exclude failed banks
    failed_bank == 0 | is.na(failed_bank),
    # Keep only banks with OMO-eligible assets (BTFP collateral)
    !is.na(omo_eligible) & omo_eligible > 0
  )

cat("✓ After filters:", nrow(baseline), "banks\n")
#> ✓ After filters: 4292 banks
cat("  - Excluded G-SIBs\n")
#>   - Excluded G-SIBs
cat("  - Excluded failed banks\n")
#>   - Excluded failed banks
cat("  - Kept only banks with BTFP-eligible assets (omo_eligible > 0)\n")
#>   - Kept only banks with BTFP-eligible assets (omo_eligible > 0)
cat("\n✓ Baseline sample created\n")
#> 
#> ✓ Baseline sample created
################################################################################
# SECTION 4: CONSTRUCT KEY RISK VARIABLES
################################################################################

cat("\n")
cat(strrep("=", 80), "\n", sep = "")
#> ================================================================================
cat("CONSTRUCTING KEY RISK VARIABLES\n")
#> CONSTRUCTING KEY RISK VARIABLES
cat(strrep("=", 80), "\n", sep = "")
#> ================================================================================
baseline <- baseline %>%
  mutate(
    # =========================================================================
    # VARIABLE CREATION AS SPECIFIED
    # =========================================================================
   
    # Borrowing subsidy (%) - valuation benefit from par pricing
    borrowing_subsidy = mtm_loss_omo_eligible_to_omo_eligible,
    
    # % MTM Loss on BTFP-eligible securities (omo_eligible)
    mtm_btfp_pct = mtm_loss_omo_eligible_to_total_asset,
    
    # % MTM Loss on Other (non-eligible) securities
    mtm_other_pct = mtm_loss_other_to_total_asset,
    
    # % Uninsured deposits
    pct_uninsured = safe_div(uninsured_deposit, total_deposit) * 100,
    
    # Uninsured deposit leverage (% of assets)
    uninsured_lev = safe_div(uninsured_deposit, total_asset) * 100,
    
    # % Wholesale Liability
    pct_wholesale_liability = safe_div(fed_fund_purchase + repo + other_borrowed_less_than_1yr, total_liability) * 100,
    
    # % Run-able Liability
    pct_runable_liability = safe_div(uninsured_deposit + fed_fund_purchase + repo + other_borrowed_less_than_1yr, total_liability) * 100,
    
    # % Liquidity Available
    pct_liquidity_available = safe_div(cash + rerepo + fed_fund_sold, total_asset) * 100,
    
    # % MTM Loss (total)
    pct_mtm_loss = mtm_loss_to_total_asset,
    
    # Cash to asset ratio (%)
    cash_to_asset = cash_to_total_asset,
    
    # FHLB to asset ratio (%)
    fhlb_to_asset = fhlb_to_total_asset,
    
    # Book equity to asset ratio (%)
    book_equity_to_asset = book_equity_to_total_asset,
    
    # Securities ratio
    securities_ratio = security_to_total_asset,
    
    # Loan to deposit
    loan_to_deposit = loan_to_deposit,
    
    # ROA (annualized)
    roa = roa,
    
    # =========================================================================
    # RUN RISK MEASURES
    # =========================================================================
    # Run-risk 1: % Uninsured × % MTM loss
    run_risk_1 = pct_uninsured * pct_mtm_loss,
    
    # Run-risk 2: % Run-able Liability × % MTM loss
    run_risk_2 = pct_runable_liability * pct_mtm_loss,
    
    # Calculate medians for run risk dummies
    median_pct_uninsured = median(pct_uninsured, na.rm = TRUE),
    median_pct_mtm_loss = median(pct_mtm_loss, na.rm = TRUE),
    median_pct_runable = median(pct_runable_liability, na.rm = TRUE),
    
    # Run Risk 1 DUMMY: = 1 if BOTH %Uninsured > median AND %MTM loss > median
    run_risk_1_dummy = as.integer(
      pct_uninsured > median_pct_uninsured & pct_mtm_loss > median_pct_mtm_loss
    ),
    
    # Run Risk 2 DUMMY: = 1 if BOTH %Runable > median AND %MTM loss > median
    run_risk_2_dummy = as.integer(
      pct_runable_liability > median_pct_runable & 
      pct_mtm_loss > median_pct_mtm_loss
    ),

    # =========================================================================
    # INSOLVENCY MEASURES
    # =========================================================================
    
    # Market value of assets
    mv_asset = mm_asset,
    mv_adjustment = if_else(mv_asset == 0, NA_real_, (total_asset / mv_asset) - 1), # Formula: (total_asset / mv_asset - 1)
    
    # Insolvency: Insolvent Banks under Different Uninsured Deposits Runs Cases (Insured Deposit Coverage Ratio)
    
    # Insured Deposit Coverage Ratio (s = 0.5)
    idcr_1 = safe_div(
      mv_asset - 0.5 * uninsured_deposit - insured_deposit,
      insured_deposit
    ),
    
    # Insured Deposit Coverage Ratio (s = 1.0)
    idcr_2 = safe_div(
      mv_asset - 1.0 * uninsured_deposit - insured_deposit,
      insured_deposit
    ),
    
    # Insolvency: Insolvent Banks under Different Uninsured Deposits Runs Cases (Capital Ratio Metric)
    
    
    insolvency_1 = safe_div(
      total_asset-total_liability - 0.5 *uninsured_deposit * mv_adjustment,
      total_asset
    ),
    
    insolvency_2 = safe_div(
      total_asset-total_liability - 1.0 *uninsured_deposit * mv_adjustment,
      total_asset
    )
  )

cat("✓ Risk variables constructed:\n")
#> ✓ Risk variables constructed:
cat("  - Borrowing subsidy\n")
#>   - Borrowing subsidy
cat("  - MTM losses (BTFP-eligible vs Other)\n")
#>   - MTM losses (BTFP-eligible vs Other)
cat("  - Uninsured deposit measures\n")
#>   - Uninsured deposit measures
cat("  - Run risk measures (continuous + dummies)\n")
#>   - Run risk measures (continuous + dummies)
cat("  - Insolvency measures (3 specifications)\n")
#>   - Insolvency measures (3 specifications)
cat("\nKey statistics:\n")
#> 
#> Key statistics:
cat("  Mean MTM BTFP-eligible (%):", sprintf("%.2f", mean(baseline$mtm_btfp_pct, na.rm = TRUE)), "\n")
#>   Mean MTM BTFP-eligible (%): 0.68
cat("  Mean MTM other (%):", sprintf("%.2f", mean(baseline$mtm_other_pct, na.rm = TRUE)), "\n")
#>   Mean MTM other (%): 4.60
cat("  Mean uninsured leverage (%):", sprintf("%.2f", mean(baseline$uninsured_lev, na.rm = TRUE)), "\n")
#>   Mean uninsured leverage (%): 23.61
cat("  % with high run risk 1:", sprintf("%.1f%%", 100 * mean(baseline$run_risk_1_dummy, na.rm = TRUE)), "\n")
#>   % with high run risk 1: 23.0%
cat("\n✓ Risk variables created\n")
#> 
#> ✓ Risk variables created
################################################################################
# SECTION 5: CREATE FHLB CRISIS BORROWING INDICATOR
################################################################################

cat("\n")
cat(strrep("=", 80), "\n", sep = "")
#> ================================================================================
cat("CREATING FHLB CRISIS BORROWING INDICATOR\n")
#> CREATING FHLB CRISIS BORROWING INDICATOR
cat(strrep("=", 80), "\n", sep = "")
#> ================================================================================
# Need historical FHLB borrowing (past 4 quarters before 2022Q4)
historical_quarters <- c("2021Q4", "2022Q1", "2022Q2", "2022Q3")

# Calculate FHLB statistics for each bank
df_fhlb_history <- call_q %>%
  filter(quarter %in% historical_quarters) %>%
  group_by(idrssd) %>%
  summarise(
    fhlb_mean_hist = mean(fhlb_to_total_asset, na.rm = TRUE),
    fhlb_sd_hist = sd(fhlb_to_total_asset, na.rm = TRUE),
    .groups = "drop"
  )

# Get crisis period FHLB borrowing (2023Q1)
df_fhlb_crisis <- call_q %>%
  filter(quarter == "2023Q1") %>%
  select(idrssd, fhlb_to_asset_crisis = fhlb_to_total_asset)

# Merge and create indicator
baseline <- baseline %>%
  left_join(df_fhlb_history, by = "idrssd") %>%
  left_join(df_fhlb_crisis, by = "idrssd") %>%
  mutate(
    # FHLB crisis indicator: 2+ std above historical mean
    fhlb_crisis = ifelse(
      !is.na(fhlb_to_asset_crisis) & !is.na(fhlb_mean_hist) & !is.na(fhlb_sd_hist),
      as.numeric(fhlb_to_asset_crisis > (fhlb_mean_hist + 2 * fhlb_sd_hist)),
      0
    ),
    
    # Alternative: change in FHLB (for continuous measure)
    delta_fhlb = ifelse(
      !is.na(fhlb_to_asset_crisis) & !is.na(fhlb_mean_hist),
      pmax(0, fhlb_to_asset_crisis - fhlb_mean_hist),
      0
    )
  )

cat("FHLB crisis borrowing:\n")
#> FHLB crisis borrowing:
cat("  Banks with crisis FHLB increase:", 
    sum(baseline$fhlb_crisis == 1, na.rm = TRUE),
    sprintf("(%.1f%%)", 100 * mean(baseline$fhlb_crisis == 1, na.rm = TRUE)), "\n")
#>   Banks with crisis FHLB increase: 1125 (26.2%)
cat("  Mean FHLB ratio (2022Q4 baseline):", 
    sprintf("%.2f%%", mean(baseline$fhlb_to_asset, na.rm = TRUE)), "\n")
#>   Mean FHLB ratio (2022Q4 baseline): 2.65%
cat("  Mean FHLB ratio (2023Q1 crisis):", 
    sprintf("%.2f%%", mean(baseline$fhlb_to_asset_crisis, na.rm = TRUE)), "\n")
#>   Mean FHLB ratio (2023Q1 crisis): 2.79%
cat("\n✓ FHLB crisis indicator created\n")
#> 
#> ✓ FHLB crisis indicator created
################################################################################
# SECTION 6: MERGE BORROWING DATA AND CREATE OUTCOMES
################################################################################

cat("\n")
cat(strrep("=", 80), "\n", sep = "")
#> ================================================================================
cat("MERGING BORROWING DATA\n")
#> MERGING BORROWING DATA
cat(strrep("=", 80), "\n", sep = "")
#> ================================================================================
# Merge BTFP and DW usage
baseline <- baseline %>%
  left_join(btfp_crisis, by = "idrssd") %>%
  left_join(dw_crisis, by = "idrssd") %>%
  left_join(btfp_full, by = "idrssd") %>%
  mutate(
    # Fill missing with 0
    btfp_user = replace_na(btfp_user, 0),
    dw_user = replace_na(dw_user, 0),
    btfp_total_amt = replace_na(btfp_total_amt, 0),
    dw_total_amt = replace_na(dw_total_amt, 0),
    
    # Create multinomial outcome
    # Y = 0: No participation
    # Y = 1: BTFP only
    # Y = 2: DW only  
    # Y = 3: Both BTFP and DW
    Y = case_when(
      btfp_user == 0 & dw_user == 0 ~ 0,
      btfp_user == 1 & dw_user == 0 ~ 1,
      btfp_user == 0 & dw_user == 1 ~ 2,
      btfp_user == 1 & dw_user == 1 ~ 3,
      TRUE ~ 0
    ),
    
    # Factor version for modeling
    Y_factor = factor(Y, 
                      levels = 0:3,
                      labels = c("No_Participation", "BTFP_Only", "DW_Only", "Both")),
    
    # Borrower indicator (used any Fed facility)
    any_borrower = as.integer(btfp_user == 1 | dw_user == 1),
    
    # BTFP vs DW choice (among borrowers only)
    # 1 = BTFP (includes "Both"), 0 = DW only
    btfp_vs_dw = ifelse(any_borrower == 1, 
                        as.integer(btfp_user == 1),
                        NA_integer_),
    
    # Borrowing amounts as % of assets
    btfp_amt_pct = safe_div(btfp_total_amt, total_asset) * 100,
    dw_amt_pct = safe_div(dw_total_amt, total_asset) * 100
  )

# Summary of outcomes
cat("\nBorrowing patterns (acute crisis):\n")
#> 
#> Borrowing patterns (acute crisis):
cat("  No participation:     ", sum(baseline$Y == 0), 
    sprintf("(%.1f%%)", 100 * mean(baseline$Y == 0)), "\n")
#>   No participation:      3544 (82.6%)
cat("  BTFP only:           ", sum(baseline$Y == 1), 
    sprintf("(%.1f%%)", 100 * mean(baseline$Y == 1)), "\n")
#>   BTFP only:            363 (8.5%)
cat("  DW only:             ", sum(baseline$Y == 2), 
    sprintf("(%.1f%%)", 100 * mean(baseline$Y == 2)), "\n")
#>   DW only:              295 (6.9%)
cat("  Both BTFP and DW:    ", sum(baseline$Y == 3), 
    sprintf("(%.1f%%)", 100 * mean(baseline$Y == 3)), "\n")
#>   Both BTFP and DW:     90 (2.1%)
cat("\nBorrower subsample:\n")
#> 
#> Borrower subsample:
cat("  Total borrowers:", sum(baseline$any_borrower == 1), "\n")
#>   Total borrowers: 748
cat("  Chose BTFP:", sum(baseline$btfp_vs_dw == 1, na.rm = TRUE), 
    sprintf("(%.1f%%)", 100 * mean(baseline$btfp_vs_dw == 1, na.rm = TRUE)), "\n")
#>   Chose BTFP: 453 (60.6%)
cat("  Chose DW only:", sum(baseline$btfp_vs_dw == 0, na.rm = TRUE), 
    sprintf("(%.1f%%)", 100 * mean(baseline$btfp_vs_dw == 0, na.rm = TRUE)), "\n")
#>   Chose DW only: 295 (39.4%)
cat("\n✓ Borrowing outcomes created\n")
#> 
#> ✓ Borrowing outcomes created
################################################################################
# SECTION 7: MERGE PUBLIC BANK DATA
################################################################################

cat("\n")
cat(strrep("=", 80), "\n", sep = "")
#> ================================================================================
cat("MERGING PUBLIC BANK DATA\n")
#> MERGING PUBLIC BANK DATA
cat(strrep("=", 80), "\n", sep = "")
#> ================================================================================
if (!is.null(df_public)) {
  # Aggregate public bank returns during crisis
  df_public_crisis <- df_public %>%
    filter(DlyCalDt >= CRISIS_START & DlyCalDt <= CRISIS_END_ACUTE) %>%
    group_by(rssd_bank) %>%
    summarise(
      n_days = n(),
      mean_ret = mean(ret, na.rm = TRUE),
      cum_ret = sum(ret, na.rm = TRUE),
      sd_ret = sd(ret, na.rm = TRUE),
      min_ret = min(ret, na.rm = TRUE),
      max_ret = max(ret, na.rm = TRUE),
      is_public = 1,
      .groups = "drop"
    ) %>%
    rename(idrssd = rssd_bank)
  
  # Merge with baseline
  baseline <- baseline %>%
    left_join(df_public_crisis, by = "idrssd") %>%
    mutate(
      is_public = replace_na(is_public, 0)
    )
  
  cat("Public bank data merged:\n")
  cat("  Public banks:", sum(baseline$is_public == 1, na.rm = TRUE), 
      sprintf("(%.1f%%)", 100 * mean(baseline$is_public == 1, na.rm = TRUE)), "\n")
  cat("  Mean cumulative return:", 
      sprintf("%.2f%%", 100 * mean(baseline$cum_ret, na.rm = TRUE)), "\n")
} else {
  baseline$is_public <- 0
  cat("  No public bank data available\n")
}
#> Public bank data merged:
#>   Public banks: 270 (6.3%) 
#>   Mean cumulative return: -13.77%
cat("\n✓ Public bank data merged\n")
#> 
#> ✓ Public bank data merged
################################################################################
# SECTION 8: WINSORIZE AND STANDARDIZE
################################################################################

cat("\n")
cat(strrep("=", 80), "\n", sep = "")
#> ================================================================================
cat("WINSORIZING AND STANDARDIZING VARIABLES\n")
#> WINSORIZING AND STANDARDIZING VARIABLES
cat(strrep("=", 80), "\n", sep = "")
#> ================================================================================
# Variables to winsorize and standardize
vars_to_process <- c(
  # MTM variables
  "mtm_btfp_pct", "mtm_other_pct", "borrowing_subsidy",
  
  # Risk variables
  "pct_uninsured", "uninsured_lev", "pct_wholesale_liability",
  "pct_runable_liability", "pct_liquidity_available", "pct_mtm_loss",
  "run_risk_1", "run_risk_2", "run_risk_1_dummy", "run_risk_2_dummy",
  
  # Insolvency
  "idcr_1", "idcr_2", "insolvency_1", "insolvency_2",
  
  # Controls
  "log_assets", "cash_to_asset", "fhlb_to_asset", "book_equity_to_asset",
  "securities_ratio", "loan_to_deposit", "roa", "delta_fhlb",
  
  # Borrowing amounts
  "btfp_amt_pct", "dw_amt_pct"
)

# Winsorize at 1st and 99th percentiles
baseline <- baseline %>%
  mutate(across(
    all_of(vars_to_process),
    ~winsorize(., probs = c(0.01, 0.99)),
    .names = "{.col}_w"
  ))

# Standardize (z-scores)
baseline <- baseline %>%
  mutate(across(
    ends_with("_w"),
    ~standardize(.),
    .names = "{.col}_z"
  ))

# Create interaction terms
baseline <- baseline %>%
  mutate(
    # MTM BTFP × Uninsured
    mtm_btfp_x_unins_w = mtm_btfp_pct_w * uninsured_lev_w,
    mtm_btfp_x_unins_w_z = standardize(mtm_btfp_x_unins_w)
  )

cat("Winsorization and standardization complete:\n")
#> Winsorization and standardization complete:
cat("  Variables processed:", length(vars_to_process), "\n")
#>   Variables processed: 27
cat("  Final sample size:", nrow(baseline), "\n")
#>   Final sample size: 4292
cat("\n✓ Variables processed\n")
#> 
#> ✓ Variables processed

5 Descriptive Statistics: Asset and Liability Composition

################################################################################
# ASSET AND LIABILITY COMPOSITION TABLES
################################################################################

cat("\n")
cat(strrep("=", 80), "\n", sep = "")
#> ================================================================================
cat("ASSET AND LIABILITY COMPOSITION TABLES\n")
#> ASSET AND LIABILITY COMPOSITION TABLES
cat(strrep("=", 80), "\n", sep = "")
#> ================================================================================
cat("\nGenerating descriptive tables for 2022Q4 and 2023Q1\n\n")
#> 
#> Generating descriptive tables for 2022Q4 and 2023Q1
# Helper function for winsorization at 5-95 percentiles
winsorize_5_95 <- function(x) {
  if (all(is.na(x))) return(x)
  q <- quantile(x, probs = c(0.05, 0.95), na.rm = TRUE, names = FALSE)
  # Clip values outside the 5th-95th range
  pmax(pmin(x, q[2]), q[1])
}

# Helper to safely sum columns (treating NA as 0)
safe_sum <- function(...) {
  rowSums(cbind(...), na.rm = TRUE)
}
# Asset Configuration
ASSET_CONFIG <- list(
  list("Cash", "cash"),
  list("Securities", "security"),
  list("\\quad Treasury", "treasury"),
  list("\\quad RMBS", "total_rmbs"),
  list("\\quad CMBS", "total_cmbs"),
  list("\\quad ABS", "abs"),
  list("\\quad Other Securities", "other_security"), # Residual
  list("Total Loans", "total_loan"),
  list("\\quad Real Estate Loans", "reloan"),
  list("\\quad\\quad Residential Mortgage", "reloan_residential1to4"),
  list("\\quad\\quad Commercial Mortgage", "reloan_residential5"),
  list("\\quad\\quad Other Real Estate", "other_re_residual"), # Residual
  list("\\quad Agricultural Loans", "agloan"),
  list("\\quad Commercial & Industrial", "ciloan"),
  list("\\quad Consumer Loans", "hhloan"),
  list("\\quad Loans to Non-Depositories", "loans_to_nonbank"),
  list("Fed Funds Sold", "fed_fund_sold"),
  list("Reverse Repo", "rerepo")
)

# Liability Configuration
LIAB_CONFIG <- list(
  list("Total Liabilities", "r_total_liability"),
  list("Domestic Deposits", "r_dom_deposit"),
  list("\\quad Insured Deposits", "r_insured_deposit"),
  list("\\quad Uninsured Deposits", "r_uninsured_deposit"),
  list("\\quad\\quad Uninsured Time Deposits (Total)", "r_uninsured_tot_time_deposit"),
  list("\\quad\\quad\\quad Long-term Time Deposits", "r_uninsured_long_time_account"),
  list("\\quad\\quad\\quad Short-term Time Deposits", "r_uninsured_short_time_deposit"),
  list("Foreign Deposits", "r_foreign_deposit"),
  list("Fed Funds Purchased", "r_fed_fund_purchase"),
  list("Repo", "r_repo"),
  list("Other Liabilities", "r_other_liab"),
  list("Total Equity", "r_total_equity"),
  list("\\quad Common Stock", "r_common_stock"),
  list("\\quad Preferred Stock", "r_preferred_stock"),
  list("\\quad Retained Earnings", "r_retained_earning")
)
# Function to prepare data (calculate residuals and ensure grouping)
prepare_data <- function(df) {
  
  # --- 1. Pre-calculation helpers ---
  # Get column or 0 if missing
  get_col <- function(var) if(var %in% names(df)) df[[var]] else 0
  
  # --- 2. Create Raw Components needed for Residuals ---
  # These sums match the "Strict Jiang" definition in the Python notebook
  # ABS = HTM Amortized + HFS Fair
  abs_val <- safe_sum(get_col("abs_htm_amortize"), get_col("abs_hfs_fair"))
  
  # RMBS = Agency (HTM+HFS) + Other (HTM+HFS)
  rmbs_val <- safe_sum(get_col("agency_rmbs_htm_amortize"), get_col("agency_rmbs_hfs_fair"),
                       get_col("other_rmbs_htm_amortize"), get_col("other_rmbs_hfs_fair"))
  
  # CMBS = Agency (HTM+HFS) + Other (HTM+HFS)
  cmbs_val <- safe_sum(get_col("agency_cmbs_htm_amortize"), get_col("agency_cmbs_hfs_fair"),
                       get_col("other_cmbs_htm_amortize"), get_col("other_cmbs_hfs_fair"))

  df %>%
    mutate(
      # --- Categories ---
      bank_category = case_when(
        size_bin == "small" ~ "Small",
        size_bin == "large" ~ "Large",
        size_bin == "gsib" ~ "GSIB",
        TRUE ~ NA_character_
      ),
      
      # --- Assign Constructed Variables ---
      abs = abs_val,
      total_rmbs = rmbs_val,
      total_cmbs = cmbs_val,
      
      # --- Calculate Asset Residuals ---
      # Other Securities = Security - (Treasury + RMBS + CMBS + ABS)
      other_security = security - safe_sum(treasury, rmbs_val, cmbs_val, abs_val),
      
      # Other Real Estate = Reloan - (Res + Comm)
      other_re_residual = reloan - safe_sum(reloan_residential1to4, reloan_residential5),
      
      # --- Calculate Missing Liability Ratios (r_*) ---
      # If these columns are missing, calculate them as: 100 * Raw / Total Assets
      r_uninsured_long_time_account = if("r_uninsured_long_time_account" %in% names(.)) 
          r_uninsured_long_time_account 
        else 
          100 * get_col("uninsured_long_time_account") / total_asset,
          
      r_uninsured_short_time_deposit = if("r_uninsured_short_time_deposit" %in% names(.)) 
          r_uninsured_short_time_deposit 
        else 
          100 * get_col("uninsured_short_time_deposit") / total_asset,
          
      r_uninsured_tot_time_deposit = if("r_uninsured_tot_time_deposit" %in% names(.)) 
          r_uninsured_tot_time_deposit 
        else 
          100 * safe_sum(get_col("uninsured_long_time_account"), get_col("uninsured_short_time_deposit")) / total_asset
    )
}

# Apply preparation
df_2022Q4 <- prepare_data(call_q %>% filter(quarter == "2022Q4"))
df_2023Q1 <- prepare_data(call_q %>% filter(quarter == "2023Q1"))

cat("✓ Data Prepared successfully.\n")
#> ✓ Data Prepared successfully.
generate_jiang_table <- function(df, config, type = "asset") {
  
  # --- 1. Calculate N (Counts) ---
  n_total <- nrow(df)
  n_small <- sum(df$bank_category == "Small", na.rm = TRUE)
  n_large <- sum(df$bank_category == "Large", na.rm = TRUE)
  n_gsib  <- sum(df$bank_category == "GSIB",  na.rm = TRUE)
  
  # Initialize DataFrame with N row
  stats_df <- data.frame(
    Component = "Number of Banks",
    Agg = format(n_total, big.mark=","), 
    Full = format(n_total, big.mark=","),
    Small = format(n_small, big.mark=","),
    Large = format(n_large, big.mark=","),
    GSIB = format(n_gsib, big.mark=","),
    stringsAsFactors = FALSE
  )
  
  total_asset_sum <- sum(df$total_asset, na.rm = TRUE)
  
  # --- 2. Process Components ---
  for (item in config) {
    label <- item[[1]]
    var_name <- item[[2]]
    
    # Get Data
    if (!var_name %in% names(df)) {
      vals <- rep(0, nrow(df))
    } else {
      vals <- df[[var_name]]
    }
    
    # Standardize to Percentage
    if (type == "asset") {
      bank_ratios <- 100 * vals / df$total_asset
      agg_val <- 100 * sum(vals, na.rm = TRUE) / total_asset_sum
    } else {
      bank_ratios <- vals
      agg_val <- sum(vals * df$total_asset, na.rm = TRUE) / total_asset_sum
    }
    
    # GLOBAL WINSORIZATION
    winsorized_ratios <- winsorize_5_95(bank_ratios)
    
    # Calculate Stats on Winsorized Vector
    calc_stat <- function(data_vec) {
      if (length(data_vec) == 0) return(c(0, 0))
      return(c(mean(data_vec, na.rm = TRUE), sd(data_vec, na.rm = TRUE)))
    }
    
    full_stats  <- calc_stat(winsorized_ratios)
    small_stats <- calc_stat(winsorized_ratios[df$bank_category == "Small"])
    large_stats <- calc_stat(winsorized_ratios[df$bank_category == "Large"])
    gsib_stats  <- calc_stat(winsorized_ratios[df$bank_category == "GSIB"])
    
    # Format Row
    fmt <- function(s) sprintf("%.1f\n(%.1f)", s[1], s[2])
    
    stats_df[nrow(stats_df) + 1, ] <- list(
      Component = label,
      Agg = sprintf("%.1f", agg_val),
      Full = fmt(full_stats),
      Small = fmt(small_stats),
      Large = fmt(large_stats),
      GSIB = fmt(gsib_stats)
    )
  }
  
  return(stats_df)
}

render_table <- function(stats_df, title) {
  stats_df %>%
    kbl(caption = title,
        col.names = c("", "Aggregate", "Full Sample", "Small", "Large", "GSIB"),
        align = c("l", "c", "c", "c", "c", "c"),
        escape = FALSE) %>%
    kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = F) %>%
    add_header_above(c(" " = 1, "(1)" = 1, "(2)" = 1, "(3)" = 1, "(4)" = 1, "(5)" = 1)) %>%
    row_spec(0, bold = TRUE) %>%
    # Add a line after the N row (row 1)
    row_spec(1, hline_after = TRUE) %>%
    footnote(general = "Values are percentages of total assets. Row 1 shows number of banks. Column (1) is aggregate. Columns (2)-(5) are bank-level averages. Standard deviations in parentheses. All bank-level statistics (Cols 2-5) are winsorized at the 5th and 95th percentiles.")
}
# --- 2022Q4 ---
asset_stats_22q4 <- generate_jiang_table(df_2022Q4, ASSET_CONFIG, type = "asset")
render_table(asset_stats_22q4, "Bank Asset Composition, 2022Q4")
Bank Asset Composition, 2022Q4
(1)
(2)
(3)
(4)
(5)
Aggregate Full Sample Small Large GSIB
Number of Banks 4,737 4,737 3,967 737 33
Cash 11.0 8.1 (7.4) 8.5 (7.5) 5.2 (5.6) 18.7 (10.1)
Securities 24.0 24.0 (14.9) 24.7 (15.2) 20.5 (12.2) 17.4 (17.4)
Treasury 6.0 3.3 (4.8) 3.5 (5.0) 2.3 (3.7) 5.1 (6.1)
RMBS 11.2 3.0 (4.3) 2.4 (3.9) 6.1 (5.1) 4.9 (6.1)
CMBS 1.8 0.8 (1.4) 0.6 (1.3) 1.3 (1.6) 0.8 (1.4)
ABS 3.1 1.0 (1.7) 0.9 (1.6) 1.5 (1.9) 1.4 (2.3)
Other Securities 1.9 14.5 (12.2) 15.9 (12.4) 7.4 (7.7) 2.3 (5.9)
Total Loans 51.0 60.8 (16.5) 59.7 (16.6) 67.4 (13.6) 43.6 (18.9)
Real Estate Loans 24.2 45.9 (18.3) 45.4 (18.2) 49.7 (17.4) 22.0 (18.6)
Residential Mortgage 11.7 17.3 (12.8) 17.7 (13.0) 15.6 (11.7) 11.8 (15.4)
Commercial Mortgage 2.5 2.3 (2.8) 2.0 (2.6) 4.0 (3.1) 1.1 (2.1)
Other Real Estate 10.0 24.8 (12.9) 24.5 (12.7) 27.9 (12.8) 4.3 (6.2)
Agricultural Loans 0.3 3.0 (4.8) 3.4 (5.1) 0.8 (2.1) 0.1 (0.4)
Commercial & Industrial 9.8 7.2 (5.4) 6.8 (5.1) 9.4 (6.4) 4.4 (5.9)
Consumer Loans 8.6 2.4 (2.7) 2.3 (2.5) 2.5 (3.4) 3.5 (4.3)
Loans to Non-Depositories 3.2 0.1 (0.2) 0.0 (0.2) 0.2 (0.3) 0.4 (0.4)
Fed Funds Sold 0.1 0.8 (1.8) 0.9 (1.9) 0.1 (0.6) 0.6 (1.9)
Reverse Repo 1.3 0.0 (0.0) 0.0 (0.0) 0.0 (0.0) 0.0 (0.0)
Note:
Values are percentages of total assets. Row 1 shows number of banks. Column (1) is aggregate. Columns (2)-(5) are bank-level averages. Standard deviations in parentheses. All bank-level statistics (Cols 2-5) are winsorized at the 5th and 95th percentiles. 
liab_stats_22q4 <- generate_jiang_table(df_2022Q4, LIAB_CONFIG, type = "liability")
render_table(liab_stats_22q4, "Bank Liability Composition, 2022Q4")
Bank Liability Composition, 2022Q4
(1)
(2)
(3)
(4)
(5)
Aggregate Full Sample Small Large GSIB
Number of Banks 4,737 4,737 3,967 737 33
Total Liabilities 90.6 90.3 (3.7) 90.3 (3.8) 90.4 (2.9) 86.9 (4.7)
Domestic Deposits 75.1 85.7 (6.4) 86.1 (6.4) 83.6 (5.7) 77.0 (7.9)
Insured Deposits 42.1 61.8 (12.7) 63.6 (11.9) 52.8 (12.0) 46.9 (18.4)
Uninsured Deposits 35.8 23.0 (11.0) 21.7 (10.4) 30.2 (11.0) 23.1 (17.0)
Uninsured Time Deposits (Total) 2.2 4.1 (3.4) 4.3 (3.4) 3.4 (3.0) 1.5 (2.4)
Long-term Time Deposits 0.3 1.1 (1.2) 1.1 (1.2) 0.8 (0.9) 0.1 (0.3)
Short-term Time Deposits 1.9 2.9 (2.6) 3.0 (2.6) 2.6 (2.4) 1.3 (2.2)
Foreign Deposits 6.3 0.0 (0.0) 0.0 (0.0) 0.0 (0.0) 0.0 (0.0)
Fed Funds Purchased 0.2 0.2 (0.5) 0.2 (0.5) 0.1 (0.4) 0.0 (0.0)
Repo 0.6 0.2 (0.7) 0.2 (0.6) 0.5 (0.9) 0.2 (0.6)
Other Liabilities 7.2 3.6 (4.0) 3.2 (3.9) 5.5 (4.1) 6.2 (4.4)
Total Equity 9.4 9.7 (3.7) 9.7 (3.8) 9.6 (2.9) 13.1 (4.7)
Common Stock 0.2 0.4 (0.6) 0.4 (0.6) 0.2 (0.6) 0.7 (1.0)
Preferred Stock 0.1 0.0 (0.0) 0.0 (0.0) 0.0 (0.0) 0.0 (0.0)
Retained Earnings 4.4 7.2 (4.2) 7.4 (4.3) 6.1 (3.2) 7.4 (5.3)
Note:
Values are percentages of total assets. Row 1 shows number of banks. Column (1) is aggregate. Columns (2)-(5) are bank-level averages. Standard deviations in parentheses. All bank-level statistics (Cols 2-5) are winsorized at the 5th and 95th percentiles. 
# --- 2023Q1 ---

asset_stats_23q1 <- generate_jiang_table(df_2023Q1, ASSET_CONFIG, type = "asset")
render_table(asset_stats_23q1, "Bank Asset Composition, 2023Q1")
Bank Asset Composition, 2023Q1
(1)
(2)
(3)
(4)
(5)
Aggregate Full Sample Small Large GSIB
Number of Banks 4,712 4,712 3,943 736 33
Cash 12.0 8.3 (7.1) 8.6 (7.2) 6.2 (5.6) 18.0 (9.5)
Securities 22.8 23.7 (14.9) 24.5 (15.2) 19.8 (11.8) 17.1 (17.1)
Treasury 5.5 3.1 (4.6) 3.3 (4.8) 2.1 (3.5) 4.9 (5.9)
RMBS 10.8 2.9 (4.3) 2.4 (3.8) 5.9 (5.0) 4.9 (6.2)
CMBS 1.7 0.8 (1.4) 0.7 (1.3) 1.3 (1.6) 0.8 (1.4)
ABS 3.0 0.9 (1.6) 0.8 (1.5) 1.5 (1.8) 1.4 (2.2)
Other Securities 1.9 14.4 (12.1) 15.9 (12.3) 7.3 (7.6) 2.2 (6.0)
Total Loans 50.6 60.8 (16.4) 59.7 (16.6) 67.4 (13.0) 44.3 (19.6)
Real Estate Loans 24.2 46.3 (18.3) 45.8 (18.3) 50.0 (17.1) 22.8 (20.3)
Residential Mortgage 11.7 17.5 (12.9) 17.9 (13.1) 15.7 (11.6) 12.0 (15.6)
Commercial Mortgage 2.5 2.3 (2.8) 2.0 (2.6) 4.1 (3.1) 1.1 (2.2)
Other Real Estate 10.0 25.0 (12.9) 24.6 (12.8) 28.0 (12.8) 4.3 (6.0)
Agricultural Loans 0.3 2.6 (4.2) 3.0 (4.5) 0.7 (1.8) 0.1 (0.3)
Commercial & Industrial 9.7 7.2 (5.4) 6.8 (5.1) 9.3 (6.3) 4.4 (5.9)
Consumer Loans 8.5 2.4 (2.7) 2.4 (2.6) 2.6 (3.5) 3.5 (4.3)
Loans to Non-Depositories 3.2 0.1 (0.2) 0.0 (0.2) 0.2 (0.3) 0.4 (0.4)
Fed Funds Sold 0.1 0.9 (1.8) 1.0 (2.0) 0.2 (0.7) 0.4 (1.6)
Reverse Repo 1.3 0.0 (0.0) 0.0 (0.0) 0.0 (0.0) 0.0 (0.0)
Note:
Values are percentages of total assets. Row 1 shows number of banks. Column (1) is aggregate. Columns (2)-(5) are bank-level averages. Standard deviations in parentheses. All bank-level statistics (Cols 2-5) are winsorized at the 5th and 95th percentiles. 
liab_stats_23q1 <- generate_jiang_table(df_2023Q1, LIAB_CONFIG, type = "liability")
render_table(liab_stats_23q1, "Bank Liability Composition, 2023Q1")
Bank Liability Composition, 2023Q1
(1)
(2)
(3)
(4)
(5)
Aggregate Full Sample Small Large GSIB
Number of Banks 4,712 4,712 3,943 736 33
Total Liabilities 90.4 90.0 (3.7) 90.0 (3.8) 90.2 (2.7) 86.5 (4.8)
Domestic Deposits 72.9 85.1 (6.5) 85.7 (6.4) 82.6 (6.1) 75.3 (7.9)
Insured Deposits 42.9 62.8 (12.1) 64.5 (11.4) 54.5 (11.2) 47.8 (16.9)
Uninsured Deposits 32.7 21.6 (10.2) 20.4 (9.8) 27.7 (10.2) 21.4 (15.9)
Uninsured Time Deposits (Total) 2.6 4.9 (3.8) 5.0 (3.8) 4.1 (3.5) 1.7 (2.7)
Long-term Time Deposits 0.3 1.1 (1.2) 1.2 (1.2) 0.8 (0.9) 0.1 (0.3)
Short-term Time Deposits 2.3 3.6 (3.1) 3.7 (3.1) 3.2 (2.8) 1.5 (2.4)
Foreign Deposits 6.1 0.0 (0.0) 0.0 (0.0) 0.0 (0.0) 0.0 (0.0)
Fed Funds Purchased 0.1 0.0 (0.2) 0.0 (0.2) 0.1 (0.2) 0.0 (0.0)
Repo 1.1 0.2 (0.6) 0.2 (0.6) 0.4 (0.8) 0.3 (0.7)
Other Liabilities 8.8 4.0 (4.3) 3.5 (4.0) 6.3 (4.6) 7.1 (5.1)
Total Equity 9.6 10.0 (3.7) 10.0 (3.8) 9.8 (2.7) 13.5 (4.8)
Common Stock 0.2 0.4 (0.7) 0.4 (0.7) 0.2 (0.6) 0.7 (1.0)
Preferred Stock 0.0 0.0 (0.0) 0.0 (0.0) 0.0 (0.0) 0.0 (0.0)
Retained Earnings 4.5 7.2 (4.2) 7.4 (4.3) 6.1 (3.2) 7.7 (5.5)
Note:
Values are percentages of total assets. Row 1 shows number of banks. Column (1) is aggregate. Columns (2)-(5) are bank-level averages. Standard deviations in parentheses. All bank-level statistics (Cols 2-5) are winsorized at the 5th and 95th percentiles.

6 Summary Statistics

################################################################################
# SECTION 9: SUMMARY STATISTICS - OVERALL
################################################################################

cat("\n")
cat(strrep("=", 80), "\n", sep = "")
#> ================================================================================
cat("SUMMARY STATISTICS\n")
#> SUMMARY STATISTICS
cat(strrep("=", 80), "\n", sep = "")
#> ================================================================================
# Variables for summary table
summ_vars <- c(
  "mtm_btfp_pct", "mtm_other_pct", "borrowing_subsidy",
  "pct_uninsured", "uninsured_lev", "pct_wholesale_liability",
  "pct_liquidity_available", "pct_mtm_loss",
  "run_risk_1", "run_risk_2",
  "idcr_1", "idcr_2", "insolvency_1", "insolvency_2",
  "log_assets", "cash_to_asset", "fhlb_to_asset", 
  "book_equity_to_asset", "securities_ratio", "loan_to_deposit", "roa"
)

# Overall summary
summ_overall <- create_summary_table(baseline, summ_vars, digits = 3)

# Print table
summ_overall %>%
  kbl(caption = "Table 1: Summary Statistics - Full Sample",
      booktabs = TRUE,
      align = c("l", rep("r", 8))) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
                full_width = FALSE,
                font_size = 11) %>%
  row_spec(0, bold = TRUE, background = "#F0F0F0") %>%
  pack_rows("MTM Losses", 1, 3) %>%
  pack_rows("Risk Measures", 4, 10) %>%
  pack_rows("Insolvency Measures", 11, 12) %>%
  pack_rows("Control Variables", 13, 19)
Table 1: Summary Statistics - Full Sample
Variable N Mean SD Min P25 Median P75 Max
MTM Losses
book_equity_to_asset 4292 10.220 8.817 0.110 7.143 8.842 10.888 99.966
borrowing_subsidy 4282 0.103 1.909 -0.034 0.032 0.065 0.098 124.830
cash_to_asset 4292 8.148 9.191 0.000 2.667 5.030 10.051 96.724
Risk Measures
fhlb_to_asset 4292 2.647 4.215 0.000 0.000 0.346 4.073 48.713
idcr_1 4282 1291.921 69365.863 -1.018 0.149 0.248 0.392 4531503.818
idcr_2 4282 1291.537 69365.870 -1.825 -0.012 0.069 0.162 4531503.818
insolvency_1 4282 0.087 0.095 -1.583 0.055 0.076 0.099 1.000
insolvency_2 4282 0.072 0.116 -3.296 0.039 0.065 0.094 1.000
loan_to_deposit 4292 70.866 27.432 0.000 56.561 71.818 86.500 890.476
log_assets 4292 12.882 1.483 8.222 11.903 12.721 13.650 20.187
Insolvency Measures
mtm_btfp_pct 4282 0.681 0.848 -0.141 0.136 0.405 0.932 16.393
mtm_other_pct 4282 4.599 2.062 -0.097 3.110 4.395 5.972 13.350
Control Variables
pct_liquidity_available 4292 9.297 9.941 0.000 3.045 5.870 11.555 98.879
pct_mtm_loss 4282 5.467 2.223 -0.141 3.863 5.320 6.975 27.631
pct_uninsured 4261 27.796 14.113 0.000 18.109 25.957 35.165 99.876
pct_wholesale_liability 4292 1.000 3.166 0.000 0.000 0.000 0.533 85.016
roa 4290 1.187 2.592 -62.459 0.702 1.101 1.498 53.461
run_risk_1 4251 147.195 89.294 -2.932 83.550 133.157 194.236 1295.061
run_risk_2 4282 143.513 86.366 -2.490 81.210 130.370 190.665 1189.577
securities_ratio 4292 25.681 15.784 0.004 13.704 23.323 35.090 97.611
uninsured_lev 4292 23.611 12.158 0.000 15.285 22.247 30.353 94.597 
cat("\n✓ Overall summary statistics created\n")
#> 
#> ✓ Overall summary statistics created
################################################################################
# SUMMARY STATISTICS BY OUTCOME
################################################################################

cat("\n")
cat(strrep("=", 80), "\n", sep = "")
#> ================================================================================
cat("SUMMARY STATISTICS BY BORROWING CHOICE\n")
#> SUMMARY STATISTICS BY BORROWING CHOICE
cat(strrep("=", 80), "\n", sep = "")
#> ================================================================================
# By borrowing status
summ_by_outcome <- baseline %>%
  select(Y_factor, all_of(summ_vars)) %>%
  pivot_longer(-Y_factor, names_to = "Variable", values_to = "Value") %>%
  group_by(Y_factor, Variable) %>%
  summarise(
    N = sum(!is.na(Value)),
    Mean = mean(Value, na.rm = TRUE),
    SD = sd(Value, na.rm = TRUE),
    .groups = "drop"
  ) %>%
  pivot_wider(
    names_from = Y_factor,
    values_from = c(N, Mean, SD),
    names_glue = "{Y_factor}_{.value}"
  ) %>%
  select(Variable, ends_with("_N"), ends_with("_Mean"), ends_with("_SD")) %>%
  mutate(across(where(is.numeric), ~round(., 3)))

# Print table (showing means only for readability)
summ_by_outcome %>%
  select(Variable, ends_with("_Mean")) %>%
  rename_with(~gsub("_Mean", "", .), ends_with("_Mean")) %>%
  kbl(caption = "Table 2: Summary Statistics by Borrowing Choice (Means)",
      booktabs = TRUE) %>%
  kable_styling(bootstrap_options = c("striped", "hover"),
                full_width = FALSE,
                font_size = 10) %>%
  row_spec(0, bold = TRUE, background = "#F0F0F0")
Table 2: Summary Statistics by Borrowing Choice (Means)
Variable No_Participation BTFP_Only DW_Only Both
book_equity_to_asset 10.582 8.166 9.017 8.227
borrowing_subsidy 0.106 0.091 0.087 0.084
cash_to_asset 8.755 4.666 6.173 4.736
fhlb_to_asset 2.420 3.959 3.270 4.255
idcr_1 1565.290 0.327 0.334 0.604
idcr_2 1564.883 0.073 0.066 0.183
insolvency_1 0.091 0.065 0.069 0.063
insolvency_2 0.077 0.049 0.048 0.044
loan_to_deposit 70.193 72.613 75.846 73.977
log_assets 12.690 13.545 13.850 14.591
mtm_btfp_pct 0.650 0.883 0.747 0.865
mtm_other_pct 4.517 5.033 4.951 4.922
pct_liquidity_available 10.030 5.118 6.953 4.987
pct_mtm_loss 5.361 6.179 5.754 5.839
pct_uninsured 26.888 30.623 31.903 38.361
pct_wholesale_liability 0.883 1.517 1.593 1.600
roa 1.196 1.083 1.212 1.182
run_risk_1 139.061 183.144 177.376 219.887
run_risk_2 135.371 179.826 175.030 213.425
securities_ratio 25.515 28.165 23.999 27.706
uninsured_lev 22.825 26.202 27.126 32.597 
cat("\n✓ Summary by outcome created\n")
#> 
#> ✓ Summary by outcome created
################################################################################
# CORRELATION MATRIX
################################################################################

cat("\n")
cat(strrep("=", 80), "\n", sep = "")
#> ================================================================================
cat("CORRELATION MATRIX\n")
#> CORRELATION MATRIX
cat(strrep("=", 80), "\n", sep = "")
#> ================================================================================
# Calculate correlations for key standardized variables
cor_vars <- c(
  "mtm_btfp_pct_w_z", "mtm_other_pct_w_z", "uninsured_lev_w_z",
  "pct_wholesale_liability_w_z", "pct_liquidity_available_w_z",
  "run_risk_1_w_z", "idcr_2_w_z", "insolvency_2_w_z",
  "log_assets_w_z", "book_equity_to_asset_w_z", "fhlb_to_asset_w_z"
)

cor_matrix <- baseline %>%
  select(all_of(cor_vars)) %>%
  cor(use = "pairwise.complete.obs")

# Clean names for display
rownames(cor_matrix) <- colnames(cor_matrix) <- c(
  "MTM BTFP", "MTM Other", "Unins Lev", "Wholesale", "Liquidity",
  "Run Risk","IDCR", "Insolvency", "Log Assets", "Equity", "FHLB"
)

# Print correlation table
cor_matrix %>%
  round(3) %>%
  kbl(caption = "Table 3: Correlation Matrix (Standardized Variables)",
      booktabs = TRUE) %>%
  kable_styling(bootstrap_options = c("striped", "hover"),
                full_width = FALSE,
                font_size = 10) %>%
  row_spec(0, bold = TRUE, background = "#F0F0F0")
Table 3: Correlation Matrix (Standardized Variables)
MTM BTFP MTM Other Unins Lev Wholesale Liquidity Run Risk IDCR Insolvency Log Assets Equity FHLB
MTM BTFP 1.000 -0.031 0.013 0.084 -0.117 0.176 -0.037 -0.062 0.070 -0.157 -0.077
MTM Other -0.031 1.000 -0.107 0.002 -0.407 0.460 -0.147 -0.122 0.113 -0.201 0.197
Unins Lev 0.013 -0.107 1.000 -0.027 0.022 0.684 -0.171 -0.413 0.441 -0.241 -0.076
Wholesale 0.084 0.002 -0.027 1.000 -0.101 0.044 -0.016 -0.029 0.085 -0.063 -0.021
Liquidity -0.117 -0.407 0.022 -0.101 1.000 -0.266 0.167 0.213 -0.333 0.253 -0.299
Run Risk 0.176 0.460 0.684 0.044 -0.266 1.000 0.058 -0.288 0.393 -0.208 0.078
IDCR -0.037 -0.147 -0.171 -0.016 0.167 0.058 1.000 0.645 -0.142 0.742 -0.061
Insolvency -0.062 -0.122 -0.413 -0.029 0.213 -0.288 0.645 1.000 -0.255 0.832 -0.129
Log Assets 0.070 0.113 0.441 0.085 -0.333 0.393 -0.142 -0.255 1.000 -0.164 0.194
Equity -0.157 -0.201 -0.241 -0.063 0.253 -0.208 0.742 0.832 -0.164 1.000 -0.061
FHLB -0.077 0.197 -0.076 -0.021 -0.299 0.078 -0.061 -0.129 0.194 -0.061 1.000 
cat("\n✓ Correlation matrix created\n")
#> 
#> ✓ Correlation matrix created

7 Model 1: Linear Probability Model (BTFP vs DW)

################################################################################
# SECTION 10: LINEAR PROBABILITY MODEL - BTFP VS DW (BORROWERS ONLY)
################################################################################

cat("\n")
cat(strrep("=", 80), "\n", sep = "")
#> ================================================================================
cat("MODEL 1: LINEAR PROBABILITY MODEL (BTFP VS DW)\n")
#> MODEL 1: LINEAR PROBABILITY MODEL (BTFP VS DW)
cat(strrep("=", 80), "\n", sep = "")
#> ================================================================================
# Create borrower subsample
df_borrowers <- baseline %>%
  filter(any_borrower == 1) %>%
  filter(complete.cases(
    btfp_vs_dw,
    mtm_btfp_pct_w_z, mtm_other_pct_w_z, uninsured_lev_w_z,
    log_assets_w_z, pct_wholesale_liability_w_z, pct_liquidity_available_w_z,
    securities_ratio_w_z, loan_to_deposit_w_z, book_equity_to_asset_w_z,
    roa_w_z, fhlb_to_asset_w_z, delta_fhlb_w_z
  ))

cat("Borrower subsample:\n")
#> Borrower subsample:
cat("  Total borrowers:", nrow(df_borrowers), "\n")
#>   Total borrowers: 748
cat("  BTFP users:", sum(df_borrowers$btfp_vs_dw == 1), "\n")
#>   BTFP users: 453
cat("  DW only:", sum(df_borrowers$btfp_vs_dw == 0), "\n\n")
#>   DW only: 295
# Specification 1: Main effects only
lpm_1 <- feols(
  btfp_vs_dw ~ 
    mtm_btfp_pct_w_z + mtm_other_pct_w_z + uninsured_lev_w_z +
    log_assets_w_z + pct_wholesale_liability_w_z + pct_liquidity_available_w_z +
    securities_ratio_w_z + loan_to_deposit_w_z + book_equity_to_asset_w_z +
    roa_w_z + fhlb_to_asset_w_z + delta_fhlb_w_z,
  data = df_borrowers,
  vcov = "hetero"
)

# Specification 2: With interaction
lpm_2 <- feols(
  btfp_vs_dw ~ 
    mtm_btfp_pct_w_z + mtm_other_pct_w_z + uninsured_lev_w_z +
    mtm_btfp_x_unins_w_z +
    log_assets_w_z + pct_wholesale_liability_w_z + pct_liquidity_available_w_z +
    securities_ratio_w_z + loan_to_deposit_w_z + book_equity_to_asset_w_z +
    roa_w_z + fhlb_to_asset_w_z + delta_fhlb_w_z,
  data = df_borrowers,
  vcov = "hetero"
)

# Specification 3: With run risk dummy
lpm_3 <- feols(
  btfp_vs_dw ~ 
    mtm_btfp_pct_w_z + mtm_other_pct_w_z + uninsured_lev_w_z +
    run_risk_1_dummy_w_z +
    log_assets_w_z + pct_wholesale_liability_w_z + pct_liquidity_available_w_z +
    securities_ratio_w_z + loan_to_deposit_w_z + book_equity_to_asset_w_z +
    roa_w_z + fhlb_to_asset_w_z + delta_fhlb_w_z,
  data = df_borrowers,
  vcov = "hetero"
)

# Display results using etable (native to fixest)
etable(
  lpm_1, lpm_2, lpm_3,
  headers = c("LPM (1)", "LPM (2)", "LPM (3)"),
  title = "Table 4: Linear Probability Model - BTFP vs DW (Borrowers Only)",
  fitstat = c("n", "r2", "ar2"),
  signif.code = c("***"=0.01, "**"=0.05, "*"=0.10),
  digits = 3
)
#>                                         lpm_1             lpm_2
#>                                       LPM (1)           LPM (2)
#> Dependent Var.:                    btfp_vs_dw        btfp_vs_dw
#>                                                                
#> Constant                     0.563*** (0.025)  0.564*** (0.025)
#> mtm_btfp_pct_w_z               -0.002 (0.019)    -0.006 (0.042)
#> mtm_other_pct_w_z              -0.034 (0.022)    -0.034 (0.022)
#> uninsured_lev_w_z               0.035 (0.021)     0.033 (0.027)
#> log_assets_w_z               -0.048** (0.023)  -0.048** (0.023)
#> pct_wholesale_liability_w_z     0.018 (0.018)     0.018 (0.018)
#> pct_liquidity_available_w_z -0.167*** (0.059) -0.167*** (0.059)
#> securities_ratio_w_z           -0.051 (0.099)    -0.051 (0.099)
#> loan_to_deposit_w_z            -0.139 (0.122)    -0.139 (0.123)
#> book_equity_to_asset_w_z       -0.034 (0.053)    -0.034 (0.053)
#> roa_w_z                        -0.013 (0.025)    -0.013 (0.025)
#> fhlb_to_asset_w_z             0.057** (0.029)   0.057** (0.029)
#> delta_fhlb_w_z                  0.025 (0.017)     0.025 (0.017)
#> mtm_btfp_x_unins_w_z                              0.004 (0.038)
#> run_risk_1_dummy_w_z                                           
#> ___________________________ _________________ _________________
#> S.E. type                   Heteroskeda.-rob. Heteroskeda.-rob.
#> Observations                              748               748
#> R2                                    0.06211           0.06212
#> Adj. R2                               0.04680           0.04551
#> 
#>                                         lpm_3
#>                                       LPM (3)
#> Dependent Var.:                    btfp_vs_dw
#>                                              
#> Constant                     0.563*** (0.025)
#> mtm_btfp_pct_w_z               -0.001 (0.019)
#> mtm_other_pct_w_z              -0.031 (0.023)
#> uninsured_lev_w_z               0.038 (0.024)
#> log_assets_w_z               -0.048** (0.023)
#> pct_wholesale_liability_w_z     0.018 (0.019)
#> pct_liquidity_available_w_z -0.169*** (0.060)
#> securities_ratio_w_z           -0.051 (0.099)
#> loan_to_deposit_w_z            -0.139 (0.123)
#> book_equity_to_asset_w_z       -0.034 (0.053)
#> roa_w_z                        -0.013 (0.025)
#> fhlb_to_asset_w_z              0.057* (0.029)
#> delta_fhlb_w_z                  0.025 (0.017)
#> mtm_btfp_x_unins_w_z                         
#> run_risk_1_dummy_w_z           -0.007 (0.019)
#> ___________________________ _________________
#> S.E. type                   Heteroskeda.-rob.
#> Observations                              748
#> R2                                    0.06230
#> Adj. R2                               0.04570
#> ---
#> Signif. codes: 0 '***' 0.01 '**' 0.05 '*' 0.1 ' ' 1
cat("\n✓ LPM models estimated\n")
#> 
#> ✓ LPM models estimated

8 Model 2: Logit (BTFP vs DW)

################################################################################
# SECTION 11: LOGIT MODEL - BTFP VS DW (BORROWERS ONLY)
################################################################################

# ==============================================================================
# MODEL 2: LOGIT (BTFP VS DW) - Re-estimated with feglm for etable compatibility
# ==============================================================================

# Specification 1: Main effects only
logit_1 <- feglm(
  btfp_vs_dw ~ 
    mtm_btfp_pct_w_z + mtm_other_pct_w_z + uninsured_lev_w_z +
    log_assets_w_z + pct_wholesale_liability_w_z + pct_liquidity_available_w_z +
    securities_ratio_w_z + loan_to_deposit_w_z + book_equity_to_asset_w_z +
    roa_w_z + fhlb_to_asset_w_z + delta_fhlb_w_z,
  data = df_borrowers,
  family = "logit", # feglm uses this syntax for logistic regression
  vcov = "hetero"   # robust standard errors
)

# Specification 2: With interaction
logit_2 <- feglm(
  btfp_vs_dw ~ 
    mtm_btfp_pct_w_z + mtm_other_pct_w_z + uninsured_lev_w_z +
    mtm_btfp_x_unins_w_z +
    log_assets_w_z + pct_wholesale_liability_w_z + pct_liquidity_available_w_z +
    securities_ratio_w_z + loan_to_deposit_w_z + book_equity_to_asset_w_z +
    roa_w_z + fhlb_to_asset_w_z + delta_fhlb_w_z,
  data = df_borrowers,
  family = "logit",
  vcov = "hetero"
)

# Specification 3: With run risk dummy
logit_3 <- feglm(
  btfp_vs_dw ~ 
    mtm_btfp_pct_w_z + mtm_other_pct_w_z + uninsured_lev_w_z +
    run_risk_1_dummy_w_z +
    log_assets_w_z + pct_wholesale_liability_w_z + pct_liquidity_available_w_z +
    securities_ratio_w_z + loan_to_deposit_w_z + book_equity_to_asset_w_z +
    roa_w_z + fhlb_to_asset_w_z + delta_fhlb_w_z,
  data = df_borrowers,
  family = "logit",
  vcov = "hetero"
)

# Display results using etable
etable(
  logit_1, logit_2, logit_3,
  headers = c("Logit (1)", "Logit (2)", "Logit (3)"),
  title = "Table 5: Logit Model - BTFP vs DW (Borrowers Only)",
  fitstat = c("n", "ll", "aic", "bic"),
  signif.code = c("***"=0.01, "**"=0.05, "*"=0.10),
  digits = 3
)
#>                                      logit_1          logit_2          logit_3
#>                                    Logit (1)        Logit (2)        Logit (3)
#> Dependent Var.:                   btfp_vs_dw       btfp_vs_dw       btfp_vs_dw
#>                                                                               
#> Constant                     0.272** (0.109)  0.273** (0.109)  0.272** (0.109)
#> mtm_btfp_pct_w_z              -0.006 (0.089)   -0.030 (0.202) -6.22e-5 (0.091)
#> mtm_other_pct_w_z             -0.153 (0.097)   -0.155 (0.098)   -0.135 (0.106)
#> uninsured_lev_w_z              0.158 (0.097)    0.148 (0.121)    0.177 (0.109)
#> log_assets_w_z              -0.218** (0.104) -0.220** (0.104) -0.221** (0.105)
#> pct_wholesale_liability_w_z    0.078 (0.086)    0.078 (0.086)    0.082 (0.087)
#> pct_liquidity_available_w_z -0.737** (0.288) -0.737** (0.288) -0.750** (0.292)
#> securities_ratio_w_z          -0.239 (0.474)   -0.238 (0.474)   -0.239 (0.477)
#> loan_to_deposit_w_z           -0.635 (0.586)   -0.634 (0.586)   -0.645 (0.591)
#> book_equity_to_asset_w_z      -0.143 (0.241)   -0.144 (0.241)   -0.144 (0.242)
#> roa_w_z                       -0.058 (0.110)   -0.058 (0.110)   -0.058 (0.111)
#> fhlb_to_asset_w_z             0.265* (0.140)   0.264* (0.140)   0.267* (0.141)
#> delta_fhlb_w_z                 0.115 (0.080)    0.114 (0.080)    0.116 (0.080)
#> mtm_btfp_x_unins_w_z                            0.025 (0.184)                 
#> run_risk_1_dummy_w_z                                            -0.038 (0.087)
#> ___________________________ ________________ ________________ ________________
#> S.E. type                   Heterosked.-rob. Heterosked.-rob. Heterosked.-rob.
#> Observations                             748              748              748
#> Log-Likelihood                       -477.77          -477.76          -477.67
#> AIC                                   981.53           983.51           983.33
#> BIC                                  1,041.6          1,048.2          1,048.0
#> ---
#> Signif. codes: 0 '***' 0.01 '**' 0.05 '*' 0.1 ' ' 1
cat("\n✓ Logit models estimated\n")
#> 
#> ✓ Logit models estimated
################################################################################
# LOGIT MARGINAL EFFECTS
################################################################################

cat("\n")
cat(strrep("=", 80), "\n", sep = "")
#> ================================================================================
cat("LOGIT MARGINAL EFFECTS\n")
#> LOGIT MARGINAL EFFECTS
cat(strrep("=", 80), "\n", sep = "")
#> ================================================================================
# Calculate average marginal effects for Logit 2
# Simple approach: numerical derivatives at mean
calc_ame_logit <- function(model, data, var_name, delta = 0.01) {
  # Baseline predictions
  pred_base <- predict(model, newdata = data, type = "response")
  
  # Create counterfactual data
  data_cf <- data
  data_cf[[var_name]] <- data_cf[[var_name]] + delta
  
  # Counterfactual predictions
  pred_cf <- predict(model, newdata = data_cf, type = "response")
  
  # Marginal effect = change in probability / delta
  me <- (pred_cf - pred_base) / delta
  
  # Average marginal effect
  ame <- mean(me, na.rm = TRUE)
  
  return(ame)
}

# Calculate AMEs for key variables
key_vars_logit <- c(
  "mtm_btfp_pct_w_z", "mtm_other_pct_w_z", "uninsured_lev_w_z",
  "mtm_btfp_x_unins_w_z", "log_assets_w_z", "book_equity_to_asset_w_z"
)

ame_logit <- data.frame(
  Variable = key_vars_logit,
  AME = sapply(key_vars_logit, function(v) calc_ame_logit(logit_2, df_borrowers, v))
) %>%
  mutate(
    Variable = gsub("_w_z", "", Variable),
    AME_pct = round(AME * 100, 2)
  )

# Display AME table
ame_logit %>%
  kbl(caption = "Table 6: Average Marginal Effects - Logit Model (Percentage Points)",
      booktabs = TRUE,
      col.names = c("Variable", "AME", "AME (%)")) %>%
  kable_styling(bootstrap_options = c("striped", "hover"),
                full_width = FALSE) %>%
  footnote(general = "AME represents change in probability of choosing BTFP for one SD increase in variable.")
Table 6: Average Marginal Effects - Logit Model (Percentage Points)
Variable AME AME (%)
mtm_btfp_pct_w_z mtm_btfp_pct -0.0068279 -0.68
mtm_other_pct_w_z mtm_other_pct -0.0346313 -3.46
uninsured_lev_w_z uninsured_lev 0.0332398 3.32
mtm_btfp_x_unins_w_z mtm_btfp_x_unins 0.0056708 0.57
log_assets_w_z log_assets -0.0492872 -4.93
book_equity_to_asset_w_z book_equity_to_asset -0.0321643 -3.22
Note:
AME represents change in probability of choosing BTFP for one SD increase in variable. 
cat("\n✓ Marginal effects calculated\n")
#> 
#> ✓ Marginal effects calculated

9 Model 3: Multinomial Logit

################################################################################
# SECTION 12: MULTINOMIAL LOGIT - FULL CHOICE SET
################################################################################

cat("\n")
cat(strrep("=", 80), "\n", sep = "")
#> ================================================================================
cat("MODEL 3: MULTINOMIAL LOGIT\n")
#> MODEL 3: MULTINOMIAL LOGIT
cat(strrep("=", 80), "\n", sep = "")
#> ================================================================================
# Create analysis dataset with complete cases
df_mlogit <- baseline %>%
  filter(complete.cases(
    Y_factor,
    mtm_btfp_pct_w_z, mtm_other_pct_w_z, uninsured_lev_w_z,
    log_assets_w_z, pct_wholesale_liability_w_z, pct_liquidity_available_w_z,
    securities_ratio_w_z, loan_to_deposit_w_z, book_equity_to_asset_w_z,
    roa_w_z, fhlb_to_asset_w_z, delta_fhlb_w_z
  ))

cat("Multinomial logit sample:\n")
#> Multinomial logit sample:
cat("  Total observations:", nrow(df_mlogit), "\n")
#>   Total observations: 4282
cat("  Outcome distribution:\n")
#>   Outcome distribution:
table(df_mlogit$Y_factor) %>% print()
#> 
#> No_Participation        BTFP_Only          DW_Only             Both 
#>             3534              363              295               90
cat("\n")
################################################################################
# ESTIMATE MULTINOMIAL LOGIT
################################################################################

# Specification 1: Main effects only
mlogit_1 <- multinom(
  Y_factor ~ 
    mtm_btfp_pct_w_z + mtm_other_pct_w_z + uninsured_lev_w_z +
    log_assets_w_z + pct_wholesale_liability_w_z + pct_liquidity_available_w_z +
    securities_ratio_w_z + loan_to_deposit_w_z + book_equity_to_asset_w_z +
    roa_w_z + fhlb_to_asset_w_z + delta_fhlb_w_z,
  data = df_mlogit,
  trace = FALSE
)

# Specification 2: With interaction
mlogit_2 <- multinom(
  Y_factor ~ 
    mtm_btfp_pct_w_z + mtm_other_pct_w_z + uninsured_lev_w_z +
    mtm_btfp_x_unins_w_z +
    log_assets_w_z + pct_wholesale_liability_w_z + pct_liquidity_available_w_z +
    securities_ratio_w_z + loan_to_deposit_w_z + book_equity_to_asset_w_z +
    roa_w_z + fhlb_to_asset_w_z + delta_fhlb_w_z,
  data = df_mlogit,
  trace = FALSE
)

# Specification 3: With run risk dummy
mlogit_3 <- multinom(
  Y_factor ~ 
    mtm_btfp_pct_w_z + mtm_other_pct_w_z + uninsured_lev_w_z +
    run_risk_1_dummy_w_z +
    log_assets_w_z + pct_wholesale_liability_w_z + pct_liquidity_available_w_z +
    securities_ratio_w_z + loan_to_deposit_w_z + book_equity_to_asset_w_z +
    roa_w_z + fhlb_to_asset_w_z + delta_fhlb_w_z,
  data = df_mlogit,
  trace = FALSE
)

cat("✓ Multinomial logit models estimated\n")
#> ✓ Multinomial logit models estimated
cat("  Model 1 AIC:", round(AIC(mlogit_1), 2), "\n")
#>   Model 1 AIC: 4926.15
cat("  Model 2 AIC:", round(AIC(mlogit_2), 2), "\n")
#>   Model 2 AIC: 4931.06
cat("  Model 3 AIC:", round(AIC(mlogit_3), 2), "\n\n")
#>   Model 3 AIC: 4924.71
################################################################################
# DISPLAY MULTINOMIAL LOGIT RESULTS
################################################################################

cat("Multinomial Logit Results:\n")
#> Multinomial Logit Results:
cat(strrep("-", 80), "\n", sep = "")
#> --------------------------------------------------------------------------------
# Get tidy output for each outcome
tidy_mlogit_2 <- broom::tidy(mlogit_2, conf.int = TRUE) %>%
  mutate(
    sig = case_when(
      p.value < 0.01 ~ "***",
      p.value < 0.05 ~ "**",
      p.value < 0.10 ~ "*",
      TRUE ~ ""
    ),
    estimate_sig = paste0(sprintf("%.3f", estimate), sig)
  )

# Print results for each outcome
for (outcome in c("BTFP_Only", "DW_Only", "Both")) {
  cat("\n", strrep("-", 80), "\n", sep = "")
  cat("Outcome:", outcome, "vs. No Participation\n")
  cat(strrep("-", 80), "\n", sep = "")
  
  tidy_mlogit_2 %>%
    filter(y.level == outcome) %>%
    select(term, estimate, std.error, statistic, p.value, sig) %>%
    mutate(across(where(is.numeric), ~round(., 4))) %>%
    kbl(booktabs = TRUE) %>%
    kable_styling(bootstrap_options = c("striped", "hover"),
                  full_width = FALSE,
                  font_size = 10) %>%
    print()
}
#> 
#> --------------------------------------------------------------------------------
#> Outcome: BTFP_Only vs. No Participation
#> --------------------------------------------------------------------------------
#> <table class="table table-striped table-hover" style="font-size: 10px; width: auto !important; margin-left: auto; margin-right: auto;">
#>  <thead>
#>   <tr>
#>    <th style="text-align:left;"> term </th>
#>    <th style="text-align:right;"> estimate </th>
#>    <th style="text-align:right;"> std.error </th>
#>    <th style="text-align:right;"> statistic </th>
#>    <th style="text-align:right;"> p.value </th>
#>    <th style="text-align:left;"> sig </th>
#>   </tr>
#>  </thead>
#> <tbody>
#>   <tr>
#>    <td style="text-align:left;"> (Intercept) </td>
#>    <td style="text-align:right;"> -2.6681 </td>
#>    <td style="text-align:right;"> 0.0817 </td>
#>    <td style="text-align:right;"> -32.6474 </td>
#>    <td style="text-align:right;"> 0.0000 </td>
#>    <td style="text-align:left;"> *** </td>
#>   </tr>
#>   <tr>
#>    <td style="text-align:left;"> mtm_btfp_pct_w_z </td>
#>    <td style="text-align:right;"> 0.1730 </td>
#>    <td style="text-align:right;"> 0.1235 </td>
#>    <td style="text-align:right;"> 1.4004 </td>
#>    <td style="text-align:right;"> 0.1614 </td>
#>    <td style="text-align:left;">  </td>
#>   </tr>
#>   <tr>
#>    <td style="text-align:left;"> mtm_other_pct_w_z </td>
#>    <td style="text-align:right;"> 0.0311 </td>
#>    <td style="text-align:right;"> 0.0669 </td>
#>    <td style="text-align:right;"> 0.4644 </td>
#>    <td style="text-align:right;"> 0.6424 </td>
#>    <td style="text-align:left;">  </td>
#>   </tr>
#>   <tr>
#>    <td style="text-align:left;"> uninsured_lev_w_z </td>
#>    <td style="text-align:right;"> 0.2226 </td>
#>    <td style="text-align:right;"> 0.0888 </td>
#>    <td style="text-align:right;"> 2.5055 </td>
#>    <td style="text-align:right;"> 0.0122 </td>
#>    <td style="text-align:left;"> ** </td>
#>   </tr>
#>   <tr>
#>    <td style="text-align:left;"> mtm_btfp_x_unins_w_z </td>
#>    <td style="text-align:right;"> -0.0948 </td>
#>    <td style="text-align:right;"> 0.1216 </td>
#>    <td style="text-align:right;"> -0.7801 </td>
#>    <td style="text-align:right;"> 0.4353 </td>
#>    <td style="text-align:left;">  </td>
#>   </tr>
#>   <tr>
#>    <td style="text-align:left;"> log_assets_w_z </td>
#>    <td style="text-align:right;"> 0.3862 </td>
#>    <td style="text-align:right;"> 0.0728 </td>
#>    <td style="text-align:right;"> 5.3066 </td>
#>    <td style="text-align:right;"> 0.0000 </td>
#>    <td style="text-align:left;"> *** </td>
#>   </tr>
#>   <tr>
#>    <td style="text-align:left;"> pct_wholesale_liability_w_z </td>
#>    <td style="text-align:right;"> 0.1665 </td>
#>    <td style="text-align:right;"> 0.0590 </td>
#>    <td style="text-align:right;"> 2.8231 </td>
#>    <td style="text-align:right;"> 0.0048 </td>
#>    <td style="text-align:left;"> *** </td>
#>   </tr>
#>   <tr>
#>    <td style="text-align:left;"> pct_liquidity_available_w_z </td>
#>    <td style="text-align:right;"> -0.6295 </td>
#>    <td style="text-align:right;"> 0.2058 </td>
#>    <td style="text-align:right;"> -3.0587 </td>
#>    <td style="text-align:right;"> 0.0022 </td>
#>    <td style="text-align:left;"> *** </td>
#>   </tr>
#>   <tr>
#>    <td style="text-align:left;"> securities_ratio_w_z </td>
#>    <td style="text-align:right;"> 0.1212 </td>
#>    <td style="text-align:right;"> 0.3195 </td>
#>    <td style="text-align:right;"> 0.3793 </td>
#>    <td style="text-align:right;"> 0.7045 </td>
#>    <td style="text-align:left;">  </td>
#>   </tr>
#>   <tr>
#>    <td style="text-align:left;"> loan_to_deposit_w_z </td>
#>    <td style="text-align:right;"> -0.0203 </td>
#>    <td style="text-align:right;"> 0.3875 </td>
#>    <td style="text-align:right;"> -0.0525 </td>
#>    <td style="text-align:right;"> 0.9582 </td>
#>    <td style="text-align:left;">  </td>
#>   </tr>
#>   <tr>
#>    <td style="text-align:left;"> book_equity_to_asset_w_z </td>
#>    <td style="text-align:right;"> -0.5408 </td>
#>    <td style="text-align:right;"> 0.1565 </td>
#>    <td style="text-align:right;"> -3.4552 </td>
#>    <td style="text-align:right;"> 0.0005 </td>
#>    <td style="text-align:left;"> *** </td>
#>   </tr>
#>   <tr>
#>    <td style="text-align:left;"> roa_w_z </td>
#>    <td style="text-align:right;"> 0.0340 </td>
#>    <td style="text-align:right;"> 0.0785 </td>
#>    <td style="text-align:right;"> 0.4335 </td>
#>    <td style="text-align:right;"> 0.6646 </td>
#>    <td style="text-align:left;">  </td>
#>   </tr>
#>   <tr>
#>    <td style="text-align:left;"> fhlb_to_asset_w_z </td>
#>    <td style="text-align:right;"> 0.2615 </td>
#>    <td style="text-align:right;"> 0.0921 </td>
#>    <td style="text-align:right;"> 2.8395 </td>
#>    <td style="text-align:right;"> 0.0045 </td>
#>    <td style="text-align:left;"> *** </td>
#>   </tr>
#>   <tr>
#>    <td style="text-align:left;"> delta_fhlb_w_z </td>
#>    <td style="text-align:right;"> -0.0415 </td>
#>    <td style="text-align:right;"> 0.0601 </td>
#>    <td style="text-align:right;"> -0.6899 </td>
#>    <td style="text-align:right;"> 0.4902 </td>
#>    <td style="text-align:left;">  </td>
#>   </tr>
#> </tbody>
#> </table>
#> --------------------------------------------------------------------------------
#> Outcome: DW_Only vs. No Participation
#> --------------------------------------------------------------------------------
#> <table class="table table-striped table-hover" style="font-size: 10px; width: auto !important; margin-left: auto; margin-right: auto;">
#>  <thead>
#>   <tr>
#>    <th style="text-align:left;"> term </th>
#>    <th style="text-align:right;"> estimate </th>
#>    <th style="text-align:right;"> std.error </th>
#>    <th style="text-align:right;"> statistic </th>
#>    <th style="text-align:right;"> p.value </th>
#>    <th style="text-align:left;"> sig </th>
#>   </tr>
#>  </thead>
#> <tbody>
#>   <tr>
#>    <td style="text-align:left;"> (Intercept) </td>
#>    <td style="text-align:right;"> -2.7535 </td>
#>    <td style="text-align:right;"> 0.0776 </td>
#>    <td style="text-align:right;"> -35.4660 </td>
#>    <td style="text-align:right;"> 0.0000 </td>
#>    <td style="text-align:left;"> *** </td>
#>   </tr>
#>   <tr>
#>    <td style="text-align:left;"> mtm_btfp_pct_w_z </td>
#>    <td style="text-align:right;"> 0.1213 </td>
#>    <td style="text-align:right;"> 0.1483 </td>
#>    <td style="text-align:right;"> 0.8180 </td>
#>    <td style="text-align:right;"> 0.4133 </td>
#>    <td style="text-align:left;">  </td>
#>   </tr>
#>   <tr>
#>    <td style="text-align:left;"> mtm_other_pct_w_z </td>
#>    <td style="text-align:right;"> 0.2020 </td>
#>    <td style="text-align:right;"> 0.0718 </td>
#>    <td style="text-align:right;"> 2.8126 </td>
#>    <td style="text-align:right;"> 0.0049 </td>
#>    <td style="text-align:left;"> *** </td>
#>   </tr>
#>   <tr>
#>    <td style="text-align:left;"> uninsured_lev_w_z </td>
#>    <td style="text-align:right;"> 0.1523 </td>
#>    <td style="text-align:right;"> 0.0882 </td>
#>    <td style="text-align:right;"> 1.7271 </td>
#>    <td style="text-align:right;"> 0.0841 </td>
#>    <td style="text-align:left;"> * </td>
#>   </tr>
#>   <tr>
#>    <td style="text-align:left;"> mtm_btfp_x_unins_w_z </td>
#>    <td style="text-align:right;"> -0.0641 </td>
#>    <td style="text-align:right;"> 0.1393 </td>
#>    <td style="text-align:right;"> -0.4599 </td>
#>    <td style="text-align:right;"> 0.6456 </td>
#>    <td style="text-align:left;">  </td>
#>   </tr>
#>   <tr>
#>    <td style="text-align:left;"> log_assets_w_z </td>
#>    <td style="text-align:right;"> 0.6817 </td>
#>    <td style="text-align:right;"> 0.0733 </td>
#>    <td style="text-align:right;"> 9.3039 </td>
#>    <td style="text-align:right;"> 0.0000 </td>
#>    <td style="text-align:left;"> *** </td>
#>   </tr>
#>   <tr>
#>    <td style="text-align:left;"> pct_wholesale_liability_w_z </td>
#>    <td style="text-align:right;"> 0.1230 </td>
#>    <td style="text-align:right;"> 0.0657 </td>
#>    <td style="text-align:right;"> 1.8731 </td>
#>    <td style="text-align:right;"> 0.0611 </td>
#>    <td style="text-align:left;"> * </td>
#>   </tr>
#>   <tr>
#>    <td style="text-align:left;"> pct_liquidity_available_w_z </td>
#>    <td style="text-align:right;"> 0.1785 </td>
#>    <td style="text-align:right;"> 0.1812 </td>
#>    <td style="text-align:right;"> 0.9850 </td>
#>    <td style="text-align:right;"> 0.3246 </td>
#>    <td style="text-align:left;">  </td>
#>   </tr>
#>   <tr>
#>    <td style="text-align:left;"> securities_ratio_w_z </td>
#>    <td style="text-align:right;"> 0.4072 </td>
#>    <td style="text-align:right;"> 0.3094 </td>
#>    <td style="text-align:right;"> 1.3164 </td>
#>    <td style="text-align:right;"> 0.1880 </td>
#>    <td style="text-align:left;">  </td>
#>   </tr>
#>   <tr>
#>    <td style="text-align:left;"> loan_to_deposit_w_z </td>
#>    <td style="text-align:right;"> 0.5799 </td>
#>    <td style="text-align:right;"> 0.3698 </td>
#>    <td style="text-align:right;"> 1.5682 </td>
#>    <td style="text-align:right;"> 0.1168 </td>
#>    <td style="text-align:left;">  </td>
#>   </tr>
#>   <tr>
#>    <td style="text-align:left;"> book_equity_to_asset_w_z </td>
#>    <td style="text-align:right;"> -0.3343 </td>
#>    <td style="text-align:right;"> 0.1542 </td>
#>    <td style="text-align:right;"> -2.1681 </td>
#>    <td style="text-align:right;"> 0.0302 </td>
#>    <td style="text-align:left;"> ** </td>
#>   </tr>
#>   <tr>
#>    <td style="text-align:left;"> roa_w_z </td>
#>    <td style="text-align:right;"> 0.0936 </td>
#>    <td style="text-align:right;"> 0.0785 </td>
#>    <td style="text-align:right;"> 1.1923 </td>
#>    <td style="text-align:right;"> 0.2332 </td>
#>    <td style="text-align:left;">  </td>
#>   </tr>
#>   <tr>
#>    <td style="text-align:left;"> fhlb_to_asset_w_z </td>
#>    <td style="text-align:right;"> 0.0452 </td>
#>    <td style="text-align:right;"> 0.0981 </td>
#>    <td style="text-align:right;"> 0.4603 </td>
#>    <td style="text-align:right;"> 0.6453 </td>
#>    <td style="text-align:left;">  </td>
#>   </tr>
#>   <tr>
#>    <td style="text-align:left;"> delta_fhlb_w_z </td>
#>    <td style="text-align:right;"> -0.1602 </td>
#>    <td style="text-align:right;"> 0.0724 </td>
#>    <td style="text-align:right;"> -2.2121 </td>
#>    <td style="text-align:right;"> 0.0270 </td>
#>    <td style="text-align:left;"> ** </td>
#>   </tr>
#> </tbody>
#> </table>
#> --------------------------------------------------------------------------------
#> Outcome: Both vs. No Participation
#> --------------------------------------------------------------------------------
#> <table class="table table-striped table-hover" style="font-size: 10px; width: auto !important; margin-left: auto; margin-right: auto;">
#>  <thead>
#>   <tr>
#>    <th style="text-align:left;"> term </th>
#>    <th style="text-align:right;"> estimate </th>
#>    <th style="text-align:right;"> std.error </th>
#>    <th style="text-align:right;"> statistic </th>
#>    <th style="text-align:right;"> p.value </th>
#>    <th style="text-align:left;"> sig </th>
#>   </tr>
#>  </thead>
#> <tbody>
#>   <tr>
#>    <td style="text-align:left;"> (Intercept) </td>
#>    <td style="text-align:right;"> -4.6738 </td>
#>    <td style="text-align:right;"> 0.2074 </td>
#>    <td style="text-align:right;"> -22.5350 </td>
#>    <td style="text-align:right;"> 0.0000 </td>
#>    <td style="text-align:left;"> *** </td>
#>   </tr>
#>   <tr>
#>    <td style="text-align:left;"> mtm_btfp_pct_w_z </td>
#>    <td style="text-align:right;"> -0.2007 </td>
#>    <td style="text-align:right;"> 0.2743 </td>
#>    <td style="text-align:right;"> -0.7318 </td>
#>    <td style="text-align:right;"> 0.4643 </td>
#>    <td style="text-align:left;">  </td>
#>   </tr>
#>   <tr>
#>    <td style="text-align:left;"> mtm_other_pct_w_z </td>
#>    <td style="text-align:right;"> 0.0577 </td>
#>    <td style="text-align:right;"> 0.1315 </td>
#>    <td style="text-align:right;"> 0.4387 </td>
#>    <td style="text-align:right;"> 0.6609 </td>
#>    <td style="text-align:left;">  </td>
#>   </tr>
#>   <tr>
#>    <td style="text-align:left;"> uninsured_lev_w_z </td>
#>    <td style="text-align:right;"> 0.3851 </td>
#>    <td style="text-align:right;"> 0.1514 </td>
#>    <td style="text-align:right;"> 2.5443 </td>
#>    <td style="text-align:right;"> 0.0109 </td>
#>    <td style="text-align:left;"> ** </td>
#>   </tr>
#>   <tr>
#>    <td style="text-align:left;"> mtm_btfp_x_unins_w_z </td>
#>    <td style="text-align:right;"> 0.0918 </td>
#>    <td style="text-align:right;"> 0.2126 </td>
#>    <td style="text-align:right;"> 0.4319 </td>
#>    <td style="text-align:right;"> 0.6658 </td>
#>    <td style="text-align:left;">  </td>
#>   </tr>
#>   <tr>
#>    <td style="text-align:left;"> log_assets_w_z </td>
#>    <td style="text-align:right;"> 0.9876 </td>
#>    <td style="text-align:right;"> 0.1216 </td>
#>    <td style="text-align:right;"> 8.1216 </td>
#>    <td style="text-align:right;"> 0.0000 </td>
#>    <td style="text-align:left;"> *** </td>
#>   </tr>
#>   <tr>
#>    <td style="text-align:left;"> pct_wholesale_liability_w_z </td>
#>    <td style="text-align:right;"> 0.0853 </td>
#>    <td style="text-align:right;"> 0.1067 </td>
#>    <td style="text-align:right;"> 0.8000 </td>
#>    <td style="text-align:right;"> 0.4237 </td>
#>    <td style="text-align:left;">  </td>
#>   </tr>
#>   <tr>
#>    <td style="text-align:left;"> pct_liquidity_available_w_z </td>
#>    <td style="text-align:right;"> 0.0764 </td>
#>    <td style="text-align:right;"> 0.3405 </td>
#>    <td style="text-align:right;"> 0.2245 </td>
#>    <td style="text-align:right;"> 0.8224 </td>
#>    <td style="text-align:left;">  </td>
#>   </tr>
#>   <tr>
#>    <td style="text-align:left;"> securities_ratio_w_z </td>
#>    <td style="text-align:right;"> 1.4846 </td>
#>    <td style="text-align:right;"> 0.5285 </td>
#>    <td style="text-align:right;"> 2.8093 </td>
#>    <td style="text-align:right;"> 0.0050 </td>
#>    <td style="text-align:left;"> *** </td>
#>   </tr>
#>   <tr>
#>    <td style="text-align:left;"> loan_to_deposit_w_z </td>
#>    <td style="text-align:right;"> 1.4533 </td>
#>    <td style="text-align:right;"> 0.6321 </td>
#>    <td style="text-align:right;"> 2.2992 </td>
#>    <td style="text-align:right;"> 0.0215 </td>
#>    <td style="text-align:left;"> ** </td>
#>   </tr>
#>   <tr>
#>    <td style="text-align:left;"> book_equity_to_asset_w_z </td>
#>    <td style="text-align:right;"> -0.9232 </td>
#>    <td style="text-align:right;"> 0.3227 </td>
#>    <td style="text-align:right;"> -2.8610 </td>
#>    <td style="text-align:right;"> 0.0042 </td>
#>    <td style="text-align:left;"> *** </td>
#>   </tr>
#>   <tr>
#>    <td style="text-align:left;"> roa_w_z </td>
#>    <td style="text-align:right;"> 0.0403 </td>
#>    <td style="text-align:right;"> 0.1457 </td>
#>    <td style="text-align:right;"> 0.2766 </td>
#>    <td style="text-align:right;"> 0.7821 </td>
#>    <td style="text-align:left;">  </td>
#>   </tr>
#>   <tr>
#>    <td style="text-align:left;"> fhlb_to_asset_w_z </td>
#>    <td style="text-align:right;"> -0.0351 </td>
#>    <td style="text-align:right;"> 0.1670 </td>
#>    <td style="text-align:right;"> -0.2102 </td>
#>    <td style="text-align:right;"> 0.8335 </td>
#>    <td style="text-align:left;">  </td>
#>   </tr>
#>   <tr>
#>    <td style="text-align:left;"> delta_fhlb_w_z </td>
#>    <td style="text-align:right;"> -0.0262 </td>
#>    <td style="text-align:right;"> 0.1060 </td>
#>    <td style="text-align:right;"> -0.2476 </td>
#>    <td style="text-align:right;"> 0.8044 </td>
#>    <td style="text-align:left;">  </td>
#>   </tr>
#> </tbody>
#> </table>
cat("\n✓ Multinomial logit results displayed\n")
#> 
#> ✓ Multinomial logit results displayed
# Display Multinomial Logit results in an 'etable-like' format


modelsummary(
  mlogit_2,
  shape = term ~ response,  # <--- This puts outcomes as columns (like etable)
  stars = c('*' = 0.10, '**' = 0.05, '***' = 0.01),
  estimate = "{estimate}{stars}",
  statistic = "std.error",  # Standard errors in parentheses
  gof_map = c("nobs", "AIC", "BIC", "logLik"),
  title = "Table 6: Multinomial Logit Results (Reference: No Participation)",
  output = "markdown"       # Use "markdown" for console text, or "kableExtra" for PDF/HTML
)
Table 6: Multinomial Logit Results (Reference: No Participation)
BTFP_Only DW_Only Both
(Intercept) -2.668*** -2.754*** -4.674***
(0.082) (0.078) (0.207)
mtm_btfp_pct_w_z 0.173 0.121 -0.201
(0.124) (0.148) (0.274)
mtm_other_pct_w_z 0.031 0.202*** 0.058
(0.067) (0.072) (0.131)
uninsured_lev_w_z 0.223** 0.152* 0.385**
(0.089) (0.088) (0.151)
mtm_btfp_x_unins_w_z -0.095 -0.064 0.092
(0.122) (0.139) (0.213)
log_assets_w_z 0.386*** 0.682*** 0.988***
(0.073) (0.073) (0.122)
pct_wholesale_liability_w_z 0.167*** 0.123* 0.085
(0.059) (0.066) (0.107)
pct_liquidity_available_w_z -0.629*** 0.178 0.076
(0.206) (0.181) (0.340)
securities_ratio_w_z 0.121 0.407 1.485***
(0.320) (0.309) (0.528)
loan_to_deposit_w_z -0.020 0.580 1.453**
(0.388) (0.370) (0.632)
book_equity_to_asset_w_z -0.541*** -0.334** -0.923***
(0.157) (0.154) (0.323)
roa_w_z 0.034 0.094 0.040
(0.078) (0.079) (0.146)
fhlb_to_asset_w_z 0.261*** 0.045 -0.035
(0.092) (0.098) (0.167)
delta_fhlb_w_z -0.041 -0.160** -0.026
(0.060) (0.072) (0.106)
Num.Obs. 4282 
cat("\n✓ Multinomial logit results displayed\n")
#> 
#> ✓ Multinomial logit results displayed
################################################################################
# MULTINOMIAL LOGIT MARGINAL EFFECTS
################################################################################

cat("\n")
cat(strrep("=", 80), "\n", sep = "")
#> ================================================================================
cat("MULTINOMIAL LOGIT MARGINAL EFFECTS\n")
#> MULTINOMIAL LOGIT MARGINAL EFFECTS
cat(strrep("=", 80), "\n", sep = "")
#> ================================================================================
# Calculate marginal effects at means
calc_mlogit_ame <- function(model, data, var_name, delta = 0.01) {
  # Baseline predictions
  pred_base <- predict(model, newdata = data, type = "probs")
  
  # Create counterfactual data
  data_cf <- data
  data_cf[[var_name]] <- data_cf[[var_name]] + delta
  
  # Counterfactual predictions
  pred_cf <- predict(model, newdata = data_cf, type = "probs")
  
  # Marginal effect = change in probability / delta
  me <- (pred_cf - pred_base) / delta
  
  # Average marginal effect for each outcome
  ame <- colMeans(me, na.rm = TRUE)
  
  return(ame)
}

# Calculate AMEs for key variables
key_vars_mlogit <- c(
  "mtm_btfp_pct_w_z", "mtm_other_pct_w_z", "uninsured_lev_w_z",
  "mtm_btfp_x_unins_w_z", "log_assets_w_z", "book_equity_to_asset_w_z",
  "fhlb_to_asset_w_z"
)

ame_list_mlogit <- lapply(key_vars_mlogit, function(v) {
  me <- calc_mlogit_ame(mlogit_2, df_mlogit, v)
  data.frame(
    Variable = v,
    No_Participation = me[1],
    BTFP_Only = me[2],
    DW_Only = me[3],
    Both = me[4]
  )
})

ame_mlogit <- do.call(rbind, ame_list_mlogit) %>%
  mutate(
    Variable = gsub("_w_z", "", Variable),
    across(where(is.numeric), ~round(. * 100, 2))  # Convert to percentage points
  )

# Display AME table
ame_mlogit %>%
  kbl(caption = "Table 7: Multinomial Logit Average Marginal Effects (Percentage Points)",
      booktabs = TRUE) %>%
  kable_styling(bootstrap_options = c("striped", "hover"),
                full_width = FALSE) %>%
  row_spec(0, bold = TRUE, background = "#F0F0F0") %>%
  footnote(general = "AME represents change in probability (in pp) for one SD increase in each variable. Row sums to zero.")
Table 7: Multinomial Logit Average Marginal Effects (Percentage Points)
Variable No_Participation BTFP_Only DW_Only Both
No_Participation mtm_btfp_pct -1.43 1.24 0.67 -0.48
No_Participation1 mtm_other_pct -1.30 0.05 1.20 0.05
No_Participation2 uninsured_lev -2.69 1.40 0.65 0.63
No_Participation3 mtm_btfp_x_unins 0.80 -0.68 -0.34 0.23
No_Participation4 log_assets -7.22 1.99 3.62 1.60
No_Participation5 book_equity_to_asset 6.31 -3.42 -1.37 -1.52
No_Participation6 fhlb_to_asset -1.82 1.90 0.08 -0.16
Note:
AME represents change in probability (in pp) for one SD increase in each variable. Row sums to zero. 
cat("\n✓ Marginal effects calculated\n")
#> 
#> ✓ Marginal effects calculated

10 IIA Tests for Multinomial Logit

################################################################################
# SECTION 13: INDEPENDENCE OF IRRELEVANT ALTERNATIVES (IIA) TESTS
################################################################################

cat("\n")
cat(strrep("=", 80), "\n", sep = "")
#> ================================================================================
cat("TESTING INDEPENDENCE OF IRRELEVANT ALTERNATIVES (IIA)\n")
#> TESTING INDEPENDENCE OF IRRELEVANT ALTERNATIVES (IIA)
cat(strrep("=", 80), "\n", sep = "")
#> ================================================================================
cat("\nThe IIA assumption states that the odds ratio between any two\n")
#> 
#> The IIA assumption states that the odds ratio between any two
cat("alternatives is independent of other alternatives:\n\n")
#> alternatives is independent of other alternatives:
cat("  P(Y=j) / P(Y=k) = exp(X'(β_j - β_k))\n\n")
#>   P(Y=j) / P(Y=k) = exp(X'(β_j - β_k))
cat("This ratio should NOT depend on the existence or attributes of\n")
#> This ratio should NOT depend on the existence or attributes of
cat("alternatives other than j and k.\n\n")
#> alternatives other than j and k.
cat("We perform two tests:\n")
#> We perform two tests:
cat("  1. Hausman-McFadden Test\n")
#>   1. Hausman-McFadden Test
cat("  2. Small-Hsiao Test\n\n")
#>   2. Small-Hsiao Test
################################################################################
# HAUSMAN-MCFADDEN TEST
################################################################################

cat(strrep("=", 80), "\n", sep = "")
#> ================================================================================
cat("HAUSMAN-MCFADDEN IIA TEST\n")
#> HAUSMAN-MCFADDEN IIA TEST
cat(strrep("=", 80), "\n", sep = "")
#> ================================================================================
cat("\nMethodology:\n")
#> 
#> Methodology:
cat("  1. Estimate full model with all J alternatives\n")
#>   1. Estimate full model with all J alternatives
cat("  2. Estimate restricted model excluding one alternative\n")
#>   2. Estimate restricted model excluding one alternative
cat("  3. Compare coefficients between full and restricted models\n")
#>   3. Compare coefficients between full and restricted models
cat("  4. Under H₀ (IIA holds): coefficients should be similar\n")
#>   4. Under H₀ (IIA holds): coefficients should be similar
cat("  5. Use AIC comparison for practical assessment\n\n")
#>   5. Use AIC comparison for practical assessment
# Function to perform Hausman-McFadden test
hausman_mcfadden_test <- function(full_model, data, formula, drop_level) {
  
  cat("\n", strrep("-", 80), "\n", sep = "")
  cat("Testing by DROPPING:", as.character(drop_level), "\n")
  cat(strrep("-", 80), "\n", sep = "")
  
  # Create restricted dataset (drop one alternative)
  data_restricted <- data %>%
    filter(Y_factor != drop_level)
  
  # Re-level factor to remove dropped level
  data_restricted$Y_factor <- droplevels(data_restricted$Y_factor)
  
  # Check if we have enough observations in each remaining category
  outcome_counts <- table(data_restricted$Y_factor)
  cat("\nRemaining outcome distribution:\n")
  print(outcome_counts)
  
  if (any(outcome_counts < 10)) {
    cat("\n⚠ WARNING: Some outcomes have < 10 observations after dropping", 
        as.character(drop_level), "\n")
    cat("   Test may be unreliable.\n")
    return(list(
      dropped = as.character(drop_level),
      test_performed = FALSE,
      reason = "Insufficient observations"
    ))
  }
  
  # Estimate restricted model
  cat("\nEstimating restricted model...\n")
  restricted_model <- tryCatch({
    multinom(formula, data = data_restricted, trace = FALSE)
  }, error = function(e) {
    cat("ERROR estimating restricted model:", e$message, "\n")
    return(NULL)
  })
  
  if (is.null(restricted_model)) {
    return(list(
      dropped = as.character(drop_level),
      test_performed = FALSE,
      reason = "Model estimation failed"
    ))
  }
  
  # Extract log-likelihoods
  ll_full <- logLik(full_model)
  ll_restr <- logLik(restricted_model)
  
  # Report results
  cat("\nFull model AIC:", sprintf("%.2f", AIC(full_model)), "\n")
  cat("Restricted model AIC:", sprintf("%.2f", AIC(restricted_model)), "\n")
  
  # AIC comparison
  aic_diff <- AIC(restricted_model) - AIC(full_model)
  cat("Δ AIC (restricted - full):", sprintf("%.2f", aic_diff), "\n")
  
  # Interpretation
  if (aic_diff < -2) {
    cat("\n✓ RESULT: IIA likely violated\n")
    iia_conclusion <- "IIA violated"
  } else if (aic_diff > 2) {
    cat("\n✓ RESULT: IIA holds\n")
    iia_conclusion <- "IIA holds"
  } else {
    cat("\n✓ RESULT: Inconclusive\n")
    iia_conclusion <- "Inconclusive"
  }
  
  return(list(
    dropped = as.character(drop_level),
    test_performed = TRUE,
    aic_full = AIC(full_model),
    aic_restr = AIC(restricted_model),
    aic_diff = aic_diff,
    conclusion = iia_conclusion
  ))
}

# Test by dropping each alternative
cat("\n")
cat(strrep("=", 80), "\n", sep = "")
#> ================================================================================
cat("HAUSMAN-MCFADDEN TESTS\n")
#> HAUSMAN-MCFADDEN TESTS
cat(strrep("=", 80), "\n", sep = "")
#> ================================================================================
# Formula
formula_mlogit_main <- Y_factor ~ 
  mtm_btfp_pct_w_z + mtm_other_pct_w_z + uninsured_lev_w_z +
  mtm_btfp_x_unins_w_z +
  log_assets_w_z + pct_wholesale_liability_w_z + pct_liquidity_available_w_z +
  securities_ratio_w_z + loan_to_deposit_w_z + book_equity_to_asset_w_z +
  roa_w_z + fhlb_to_asset_w_z + delta_fhlb_w_z

alternatives_to_drop <- c("BTFP_Only", "DW_Only", "Both")

iia_results <- lapply(alternatives_to_drop, function(alt) {
  hausman_mcfadden_test(mlogit_2, df_mlogit, formula_mlogit_main, alt)
})
#> 
#> --------------------------------------------------------------------------------
#> Testing by DROPPING: BTFP_Only 
#> --------------------------------------------------------------------------------
#> 
#> Remaining outcome distribution:
#> 
#> No_Participation          DW_Only             Both 
#>             3534              295               90 
#> 
#> Estimating restricted model...
#> 
#> Full model AIC: 4931.06 
#> Restricted model AIC: 2633.96 
#> Δ AIC (restricted - full): -2297.10 
#> 
#> ✓ RESULT: IIA likely violated
#> 
#> --------------------------------------------------------------------------------
#> Testing by DROPPING: DW_Only 
#> --------------------------------------------------------------------------------
#> 
#> Remaining outcome distribution:
#> 
#> No_Participation        BTFP_Only             Both 
#>             3534              363               90 
#> 
#> Estimating restricted model...
#> 
#> Full model AIC: 4931.06 
#> Restricted model AIC: 2913.52 
#> Δ AIC (restricted - full): -2017.55 
#> 
#> ✓ RESULT: IIA likely violated
#> 
#> --------------------------------------------------------------------------------
#> Testing by DROPPING: Both 
#> --------------------------------------------------------------------------------
#> 
#> Remaining outcome distribution:
#> 
#> No_Participation        BTFP_Only          DW_Only 
#>             3534              363              295 
#> 
#> Estimating restricted model...
#> 
#> Full model AIC: 4931.06 
#> Restricted model AIC: 4174.62 
#> Δ AIC (restricted - full): -756.45 
#> 
#> ✓ RESULT: IIA likely violated
names(iia_results) <- alternatives_to_drop
# Create summary table
iia_summary <- data.frame(
  Alternative_Dropped = sapply(iia_results, function(x) x$dropped),
  Test_Performed = sapply(iia_results, function(x) x$test_performed),
  Full_AIC = sapply(iia_results, function(x) {
    if (x$test_performed) round(x$aic_full, 2) else NA
  }),
  Restricted_AIC = sapply(iia_results, function(x) {
    if (x$test_performed) round(x$aic_restr, 2) else NA
  }),
  Delta_AIC = sapply(iia_results, function(x) {
    if (x$test_performed) round(x$aic_diff, 2) else NA
  }),
  Conclusion = sapply(iia_results, function(x) {
    if (x$test_performed) x$conclusion else x$reason
  })
)

iia_summary %>%
  kbl(caption = "Table 10: Hausman-McFadden IIA Test Results",
      booktabs = TRUE) %>%
  kable_styling(bootstrap_options = c("striped", "hover"),
                full_width = FALSE) %>%
  row_spec(0, bold = TRUE, background = "#F0F0F0") %>%
  footnote(general = c(
    "Δ AIC = Restricted AIC - Full AIC",
    "Δ AIC < -2: IIA violated",
    "Δ AIC > +2: IIA holds"
  ))
Table 10: Hausman-McFadden IIA Test Results
Alternative_Dropped Test_Performed Full_AIC Restricted_AIC Delta_AIC Conclusion
BTFP_Only BTFP_Only TRUE 4931.06 2633.96 -2297.10 IIA violated
DW_Only DW_Only TRUE 4931.06 2913.52 -2017.55 IIA violated
Both Both TRUE 4931.06 4174.62 -756.45 IIA violated
Note:
Δ AIC = Restricted AIC - Full AIC
Δ AIC < -2: IIA violated
Δ AIC > +2: IIA holds 
################################################################################
# SMALL-HSIAO TEST
################################################################################

cat("\n")
cat(strrep("=", 80), "\n", sep = "")
#> ================================================================================
cat("SMALL-HSIAO IIA TEST\n")
#> SMALL-HSIAO IIA TEST
cat(strrep("=", 80), "\n", sep = "")
#> ================================================================================
small_hsiao_test <- function(data, formula, drop_level, split_ratio = 0.5, seed = 123) {
  
  set.seed(seed)
  
  # Split data
  n <- nrow(data)
  n1 <- floor(n * split_ratio)
  idx1 <- sample(1:n, n1)
  
  data_s1 <- data[idx1, ]
  data_s2 <- data[-idx1, ]
  
  cat("\nTesting:", as.character(drop_level), "\n")
  
  # Estimate models
  model_s2_full <- tryCatch({
    multinom(formula, data = data_s2, trace = FALSE)
  }, error = function(e) return(NULL))
  
  if (is.null(model_s2_full)) {
    return(list(test_performed = FALSE, reason = "Model failed"))
  }
  
  data_s2_restr <- data_s2 %>% filter(Y_factor != drop_level)
  data_s2_restr$Y_factor <- droplevels(data_s2_restr$Y_factor)
  
  model_s2_restr <- tryCatch({
    multinom(formula, data = data_s2_restr, trace = FALSE)
  }, error = function(e) return(NULL))
  
  if (is.null(model_s2_restr)) {
    return(list(test_performed = FALSE, reason = "Restricted model failed"))
  }
  
  # Compare AICs
  aic_diff <- AIC(model_s2_restr) - AIC(model_s2_full)
  
  cat("Δ AIC:", sprintf("%.2f", aic_diff), "\n")
  
  conclusion <- if (aic_diff < -2) "IIA violated" else if (aic_diff > 2) "IIA holds" else "Inconclusive"
  
  return(list(
    test_performed = TRUE,
    dropped = as.character(drop_level),
    aic_diff = aic_diff,
    conclusion = conclusion
  ))
}

sh_results <- lapply(alternatives_to_drop, function(alt) {
  small_hsiao_test(df_mlogit, formula_mlogit_main, alt, seed = 20230313)
})
#> 
#> Testing: BTFP_Only 
#> Δ AIC: -1138.41 
#> 
#> Testing: DW_Only 
#> Δ AIC: -997.93 
#> 
#> Testing: Both 
#> Δ AIC: -403.83
names(sh_results) <- alternatives_to_drop
sh_summary <- data.frame(
  Alternative_Dropped = sapply(sh_results, function(x) {
    if (x$test_performed) x$dropped else "N/A"
  }),
  Test_Performed = sapply(sh_results, function(x) x$test_performed),
  Delta_AIC = sapply(sh_results, function(x) {
    if (x$test_performed) round(x$aic_diff, 2) else NA
  }),
  Conclusion = sapply(sh_results, function(x) {
    if (x$test_performed) x$conclusion else x$reason
  })
)

sh_summary %>%
  kbl(caption = "Table 11: Small-Hsiao IIA Test Results",
      booktabs = TRUE) %>%
  kable_styling(bootstrap_options = c("striped", "hover"),
                full_width = FALSE) %>%
  row_spec(0, bold = TRUE, background = "#F0F0F0")
Table 11: Small-Hsiao IIA Test Results
Alternative_Dropped Test_Performed Delta_AIC Conclusion
BTFP_Only BTFP_Only TRUE -1138.41 IIA violated
DW_Only DW_Only TRUE -997.93 IIA violated
Both Both TRUE -403.83 IIA violated 
################################################################################
# OVERALL CONCLUSION
################################################################################

cat("\n")
cat(strrep("=", 80), "\n", sep = "")
#> ================================================================================
cat("IIA TEST CONCLUSION\n")
#> IIA TEST CONCLUSION
cat(strrep("=", 80), "\n", sep = "")
#> ================================================================================
hm_violations <- sum(grepl("violated", iia_summary$Conclusion), na.rm = TRUE)
sh_violations <- sum(grepl("violated", sh_summary$Conclusion), na.rm = TRUE)

cat("\nHausman-McFadden: IIA violations =", hm_violations, "\n")
#> 
#> Hausman-McFadden: IIA violations = 3
cat("Small-Hsiao: IIA violations =", sh_violations, "\n\n")
#> Small-Hsiao: IIA violations = 3
if (hm_violations > 0 || sh_violations > 0) {
  cat("⚠ CONCLUSION: IIA assumption VIOLATED for some alternatives.\n\n")
  cat("RECOMMENDATIONS:\n")
  cat("  1. Consider nested logit model\n")
  cat("  2. Report both multinomial and nested logit results\n")
  cat("  3. Interpret with caution\n")
} else {
  cat("✓ CONCLUSION: IIA assumption appears to HOLD.\n")
  cat("  Multinomial logit is appropriate.\n")
}
#> ⚠ CONCLUSION: IIA assumption VIOLATED for some alternatives.
#> 
#> RECOMMENDATIONS:
#>   1. Consider nested logit model
#>   2. Report both multinomial and nested logit results
#>   3. Interpret with caution
cat("\n✓ IIA testing complete\n")
#> 
#> ✓ IIA testing complete

11 Model 4: Trivariate Probit

# ==============================================================================
# MODEL 4: TRIVARIATE PROBIT (APPROXIMATION VIA SEPARATE PROBITS)
# ==============================================================================

cat("\n")
cat(strrep("=", 80), "\n", sep = "")
#> ================================================================================
cat("MODEL 4: SEPARATE PROBITS (PROXY FOR TRIVARIATE)\n")
#> MODEL 4: SEPARATE PROBITS (PROXY FOR TRIVARIATE)
cat(strrep("=", 80), "\n", sep = "")
#> ================================================================================
# 1. Define Formulas
formula_btfp <- btfp_user ~ 
  mtm_btfp_pct_w_z + mtm_other_pct_w_z + uninsured_lev_w_z +
  mtm_btfp_x_unins_w_z +
  log_assets_w_z + pct_wholesale_liability_w_z + pct_liquidity_available_w_z +
  securities_ratio_w_z + loan_to_deposit_w_z + book_equity_to_asset_w_z +
  roa_w_z + fhlb_to_asset_w_z + delta_fhlb_w_z

formula_dw <- dw_user ~ 
  mtm_btfp_pct_w_z + mtm_other_pct_w_z + uninsured_lev_w_z +
  mtm_btfp_x_unins_w_z +
  log_assets_w_z + pct_wholesale_liability_w_z + pct_liquidity_available_w_z +
  securities_ratio_w_z + loan_to_deposit_w_z + book_equity_to_asset_w_z +
  roa_w_z + fhlb_to_asset_w_z + delta_fhlb_w_z

formula_fhlb <- fhlb_crisis ~ 
  mtm_btfp_pct_w_z + mtm_other_pct_w_z + uninsured_lev_w_z +
  mtm_btfp_x_unins_w_z +
  log_assets_w_z + pct_wholesale_liability_w_z + pct_liquidity_available_w_z +
  securities_ratio_w_z + loan_to_deposit_w_z + book_equity_to_asset_w_z +
  roa_w_z + fhlb_to_asset_w_z + delta_fhlb_w_z

# 2. Estimate using feglm (fixest) for etable compatibility
# Note: family = binomial(link = "probit") sets this to Probit
probit_btfp <- feglm(formula_btfp, data = df_mlogit, family = binomial(link = "probit"))
probit_dw   <- feglm(formula_dw,   data = df_mlogit, family = binomial(link = "probit"))
probit_fhlb <- feglm(formula_fhlb, data = df_mlogit, family = binomial(link = "probit"))

cat("✓ Separate probit models estimated with feglm\n\n")
#> ✓ Separate probit models estimated with feglm
# 3. Display results using etable
etable(
  probit_btfp, probit_dw, probit_fhlb,
  headers = c("BTFP", "DW", "FHLB Crisis"),
  title = "Table 8: Separate Probit Models (BTFP, DW, FHLB)",
  fitstat = c("n", "ll", "aic", "bic"),
  signif.code = c("***"=0.01, "**"=0.05, "*"=0.10),
  digits = 3
)
#>                                   probit_btfp         probit_dw
#>                                          BTFP                DW
#> Dependent Var.:                     btfp_user           dw_user
#>                                                                
#> Constant                     -1.46*** (0.036)  -1.51*** (0.034)
#> mtm_btfp_pct_w_z                0.056 (0.062)     0.022 (0.068)
#> mtm_other_pct_w_z               0.008 (0.032)   0.079** (0.033)
#> uninsured_lev_w_z            0.130*** (0.042)   0.082** (0.041)
#> mtm_btfp_x_unins_w_z           -0.029 (0.060)    -0.009 (0.064)
#> log_assets_w_z               0.227*** (0.034)  0.380*** (0.035)
#> pct_wholesale_liability_w_z  0.080*** (0.029)     0.042 (0.031)
#> pct_liquidity_available_w_z -0.250*** (0.095)    0.159* (0.087)
#> securities_ratio_w_z            0.144 (0.153)   0.363** (0.147)
#> loan_to_deposit_w_z             0.059 (0.187)   0.444** (0.179)
#> book_equity_to_asset_w_z    -0.272*** (0.074) -0.205*** (0.073)
#> roa_w_z                         0.011 (0.037)     0.031 (0.037)
#> fhlb_to_asset_w_z            0.122*** (0.045)    -0.010 (0.046)
#> delta_fhlb_w_z                 -0.007 (0.030)   -0.061* (0.033)
#> ___________________________ _________________ _________________
#> S.E. type                                 IID               IID
#> Observations                            4,282             4,282
#> Log-Likelihood                       -1,278.9          -1,153.4
#> AIC                                   2,585.9           2,334.9
#> BIC                                   2,675.0           2,423.9
#> 
#>                                   probit_fhlb
#>                                   FHLB Crisis
#> Dependent Var.:                   fhlb_crisis
#>                                              
#> Constant                    -0.841*** (0.032)
#> mtm_btfp_pct_w_z               0.109* (0.063)
#> mtm_other_pct_w_z            0.119*** (0.032)
#> uninsured_lev_w_z              -0.040 (0.043)
#> mtm_btfp_x_unins_w_z           -0.097 (0.070)
#> log_assets_w_z               0.114*** (0.037)
#> pct_wholesale_liability_w_z     0.021 (0.031)
#> pct_liquidity_available_w_z -0.335*** (0.080)
#> securities_ratio_w_z           -0.147 (0.124)
#> loan_to_deposit_w_z            -0.062 (0.144)
#> book_equity_to_asset_w_z       -0.053 (0.046)
#> roa_w_z                        -0.003 (0.035)
#> fhlb_to_asset_w_z           -0.490*** (0.049)
#> delta_fhlb_w_z                1.56*** (0.055)
#> ___________________________ _________________
#> S.E. type                                 IID
#> Observations                            4,282
#> Log-Likelihood                       -1,214.9
#> AIC                                   2,457.8
#> BIC                                   2,546.9
#> ---
#> Signif. codes: 0 '***' 0.01 '**' 0.05 '*' 0.1 ' ' 1
# 4. Calculate residual correlations (fixest supports residuals())
cat("\nCalculating residual correlations...\n")
#> 
#> Calculating residual correlations...
# type = "response" gives raw residuals (y - p_hat)
resid_btfp <- residuals(probit_btfp, type = "response")
resid_dw   <- residuals(probit_dw,   type = "response")
resid_fhlb <- residuals(probit_fhlb, type = "response")

cor_matrix_resid <- cor(data.frame(
  BTFP = resid_btfp,
  DW   = resid_dw,
  FHLB = resid_fhlb
), use = "pairwise.complete.obs")

cat("\nResidual correlations (proxy for error correlations):\n")
#> 
#> Residual correlations (proxy for error correlations):
print(round(cor_matrix_resid, 3))
#>       BTFP     DW   FHLB
#> BTFP 1.000  0.073  0.017
#> DW   0.073  1.000 -0.020
#> FHLB 0.017 -0.020  1.000
cat("\n")
cat("Interpretation:\n")
#> Interpretation:
cat("  ρ(BTFP, DW):", sprintf("%.3f", cor_matrix_resid[1,2]), 
    "- Correlation in unobserved factors affecting both\n")
#>   ρ(BTFP, DW): 0.073 - Correlation in unobserved factors affecting both
cat("  ρ(BTFP, FHLB):", sprintf("%.3f", cor_matrix_resid[1,3]), 
    "- Potential complementarity/substitution\n")
#>   ρ(BTFP, FHLB): 0.017 - Potential complementarity/substitution
cat("  ρ(DW, FHLB):", sprintf("%.3f", cor_matrix_resid[2,3]), 
    "- Potential complementarity/substitution\n")
#>   ρ(DW, FHLB): -0.020 - Potential complementarity/substitution
cat("\n✓ Trivariate probit (separate estimation) completed\n")
#> 
#> ✓ Trivariate probit (separate estimation) completed

12 Robustness Checks

################################################################################
# SECTION 14: SUBSAMPLE ANALYSIS
################################################################################

cat("\n")
cat(strrep("=", 80), "\n", sep = "")
#> ================================================================================
cat("SUBSAMPLE ANALYSIS\n")
#> SUBSAMPLE ANALYSIS
cat(strrep("=", 80), "\n", sep = "")
#> ================================================================================
# Size-based subsamples
median_assets <- median(df_mlogit$log_assets_w_z, na.rm = TRUE)

cat("Estimating for small banks (below median assets)...\n")
#> Estimating for small banks (below median assets)...
df_small <- df_mlogit %>% filter(log_assets_w_z <= median_assets)
mlogit_small <- multinom(
  Y_factor ~ 
    mtm_btfp_pct_w_z + mtm_other_pct_w_z + uninsured_lev_w_z +
    mtm_btfp_x_unins_w_z +
    log_assets_w_z + pct_wholesale_liability_w_z + pct_liquidity_available_w_z +
    book_equity_to_asset_w_z + fhlb_to_asset_w_z,
  data = df_small,
  trace = FALSE
)

cat("Estimating for large banks (above median assets)...\n")
#> Estimating for large banks (above median assets)...
df_large <- df_mlogit %>% filter(log_assets_w_z > median_assets)
mlogit_large <- multinom(
  Y_factor ~ 
    mtm_btfp_pct_w_z + mtm_other_pct_w_z + uninsured_lev_w_z +
    mtm_btfp_x_unins_w_z +
    log_assets_w_z + pct_wholesale_liability_w_z + pct_liquidity_available_w_z +
    book_equity_to_asset_w_z + fhlb_to_asset_w_z,
  data = df_large,
  trace = FALSE
)

cat("\n✓ Size-based subsample analysis completed\n")
#> 
#> ✓ Size-based subsample analysis completed
# Compare key coefficients across subsamples
cat("\nComparing MTM BTFP coefficient (BTFP_Only outcome) across subsamples:\n")
#> 
#> Comparing MTM BTFP coefficient (BTFP_Only outcome) across subsamples:
cat("  Full sample:", 
    sprintf("%.4f", coef(mlogit_2)["BTFP_Only", "mtm_btfp_pct_w_z"]), "\n")
#>   Full sample: 0.1730
cat("  Small banks:", 
    sprintf("%.4f", coef(mlogit_small)["BTFP_Only", "mtm_btfp_pct_w_z"]), "\n")
#>   Small banks: 0.5619
cat("  Large banks:", 
    sprintf("%.4f", coef(mlogit_large)["BTFP_Only", "mtm_btfp_pct_w_z"]), "\n")
#>   Large banks: 0.0328

13 Visualizations

################################################################################
# SECTION 15: VISUALIZATIONS
################################################################################

cat("\n")
cat(strrep("=", 80), "\n", sep = "")
#> ================================================================================
cat("CREATING VISUALIZATIONS\n")
#> CREATING VISUALIZATIONS
cat(strrep("=", 80), "\n", sep = "")
#> ================================================================================
# Calculate predicted probabilities
pred_probs <- predict(mlogit_2, newdata = df_mlogit, type = "probs")
df_mlogit$pred_no_part <- pred_probs[, "No_Participation"]
df_mlogit$pred_btfp <- pred_probs[, "BTFP_Only"]
df_mlogit$pred_dw <- pred_probs[, "DW_Only"]
df_mlogit$pred_both <- pred_probs[, "Both"]

# Plot 1: Distribution of predicted probabilities by actual outcome
plot_data <- df_mlogit %>%
  select(Y_factor, starts_with("pred_")) %>%
  pivot_longer(starts_with("pred_"), 
               names_to = "Predicted", 
               values_to = "Probability") %>%
  mutate(
    Predicted = gsub("pred_", "", Predicted),
    Predicted = case_when(
      Predicted == "no_part" ~ "No Participation",
      Predicted == "btfp" ~ "BTFP Only",
      Predicted == "dw" ~ "DW Only",
      Predicted == "both" ~ "Both",
      TRUE ~ Predicted
    )
  )

p1 <- ggplot(plot_data, aes(x = Probability, fill = Y_factor)) +
  geom_density(alpha = 0.5) +
  facet_wrap(~Predicted, scales = "free") +
  theme_minimal() +
  labs(
    title = "Distribution of Predicted Probabilities by Actual Outcome",
    x = "Predicted Probability",
    y = "Density",
    fill = "Actual Outcome"
  ) +
  theme(legend.position = "bottom",
        plot.title = element_text(hjust = 0.5, face = "bold"))

print(p1)

ggsave(file.path(FIG_PATH, "predicted_probabilities_distribution.png"), 
       p1, width = 12, height = 8, dpi = 300)

cat("\n✓ Predicted probability plots created\n")
#> 
#> ✓ Predicted probability plots created
################################################################################
# PLOT MARGINAL EFFECTS
################################################################################

# Create plot of marginal effects
plot_me <- ame_mlogit %>%
  pivot_longer(-Variable, names_to = "Outcome", values_to = "Effect") %>%
  mutate(
    Outcome = factor(Outcome, 
                    levels = c("No_Participation", "BTFP_Only", "DW_Only", "Both"),
                    labels = c("No Participation", "BTFP Only", "DW Only", "Both"))
  )

p2 <- ggplot(plot_me, aes(x = Variable, y = Effect, fill = Outcome)) +
  geom_col(position = "dodge") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray50") +
  coord_flip() +
  theme_minimal() +
  labs(
    title = "Average Marginal Effects on Choice Probabilities",
    subtitle = "Effect of one standard deviation increase (percentage points)",
    x = NULL,
    y = "Marginal Effect (percentage points)",
    fill = "Outcome"
  ) +
  theme(legend.position = "bottom",
        plot.title = element_text(hjust = 0.5, face = "bold"),
        plot.subtitle = element_text(hjust = 0.5))

print(p2)

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

cat("\n✓ Marginal effects plot created\n")
#> 
#> ✓ Marginal effects plot created
################################################################################
# PLOT MTM EFFECTS BY OUTCOME
################################################################################

# Create grid of MTM values
mtm_grid <- seq(
  min(df_mlogit$mtm_btfp_pct_w_z, na.rm = TRUE),
  max(df_mlogit$mtm_btfp_pct_w_z, na.rm = TRUE),
  length.out = 50
)

# Create counterfactual data (hold other variables at mean = 0)
data_cf <- df_mlogit[rep(1, 50), ]
for (i in 1:50) {
  data_cf[i, "mtm_btfp_pct_w_z"] <- mtm_grid[i]
  # Set other standardized variables to mean (= 0)
  for (v in setdiff(key_vars_mlogit, "mtm_btfp_pct_w_z")) {
    if (v %in% names(data_cf)) {
      data_cf[i, v] <- 0
    }
  }
}

# Predict probabilities
pred_mtm <- predict(mlogit_2, newdata = data_cf, type = "probs")

plot_data_mtm <- data.frame(
  mtm_btfp = mtm_grid,
  No_Participation = pred_mtm[, "No_Participation"],
  BTFP_Only = pred_mtm[, "BTFP_Only"],
  DW_Only = pred_mtm[, "DW_Only"],
  Both = pred_mtm[, "Both"]
) %>%
  pivot_longer(-mtm_btfp, names_to = "Outcome", values_to = "Probability")

p3 <- ggplot(plot_data_mtm, aes(x = mtm_btfp, y = Probability, color = Outcome)) +
  geom_line(size = 1.2) +
  theme_minimal() +
  labs(
    title = "Predicted Probability by MTM Loss on BTFP-Eligible Assets",
    subtitle = "All other variables held at mean",
    x = "MTM Loss on BTFP-Eligible Assets (standardized)",
    y = "Predicted Probability",
    color = "Outcome"
  ) +
  scale_y_continuous(labels = percent_format()) +
  theme(legend.position = "bottom",
        plot.title = element_text(hjust = 0.5, face = "bold"),
        plot.subtitle = element_text(hjust = 0.5))

print(p3)

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

cat("\n✓ MTM effects plot created\n")
#> 
#> ✓ MTM effects plot created

14 Export Results

################################################################################
# SECTION 16: EXPORT RESULTS
################################################################################

cat("\n")
cat(strrep("=", 80), "\n", sep = "")
#> ================================================================================
cat("EXPORTING RESULTS\n")
#> EXPORTING RESULTS
cat(strrep("=", 80), "\n", sep = "")
#> ================================================================================
# Save regression tables
modelsummary(
  list(
    "LPM (2)" = lpm_2,
    "Logit (2)" = logit_2,
    "MLogit BTFP" = mlogit_2,
    "MLogit DW" = mlogit_2,
    "MLogit Both" = mlogit_2
  ),
  output = file.path(TABLE_PATH, "all_models_comparison.html"),
  estimate = "{estimate}{stars}",
  statistic = "std.error",
  stars = c('*' = 0.10, '**' = 0.05, '***' = 0.01)
)

# Save marginal effects
write.csv(
  ame_logit,
  file.path(TABLE_PATH, "marginal_effects_logit.csv"),
  row.names = FALSE
)

write.csv(
  ame_mlogit,
  file.path(TABLE_PATH, "marginal_effects_mlogit.csv"),
  row.names = FALSE
)

# Save summary statistics
write.csv(
  summ_overall,
  file.path(TABLE_PATH, "summary_statistics_overall.csv"),
  row.names = FALSE
)

write.csv(
  summ_by_outcome,
  file.path(TABLE_PATH, "summary_statistics_by_outcome.csv"),
  row.names = FALSE
)

# Save model objects for later use
saveRDS(lpm_2, file.path(TABLE_PATH, "lpm_main.rds"))
saveRDS(logit_2, file.path(TABLE_PATH, "logit_main.rds"))
saveRDS(mlogit_2, file.path(TABLE_PATH, "mlogit_main.rds"))
saveRDS(list(probit_btfp, probit_dw, probit_fhlb), 
        file.path(TABLE_PATH, "probit_models.rds"))

cat("\n✓ Results exported to:", TABLE_PATH, "\n")
#> 
#> ✓ Results exported to: C:/Users/mferdo2/OneDrive - Louisiana State University/Finance_PhD/DW_Stigma_paper/Liquidity_project_2025/03_documentation/regression_tables/multinomial_analysis
cat("  Tables:\n")
#>   Tables:
cat("    - all_models_comparison.html\n")
#>     - all_models_comparison.html
cat("    - marginal_effects_logit.csv\n")
#>     - marginal_effects_logit.csv
cat("    - marginal_effects_mlogit.csv\n")
#>     - marginal_effects_mlogit.csv
cat("    - summary_statistics_*.csv\n")
#>     - summary_statistics_*.csv
cat("  Model objects:\n")
#>   Model objects:
cat("    - lpm_main.rds\n")
#>     - lpm_main.rds
cat("    - logit_main.rds\n")
#>     - logit_main.rds
cat("    - mlogit_main.rds\n")
#>     - mlogit_main.rds
cat("    - probit_models.rds\n")
#>     - probit_models.rds
cat("  Figures:\n")
#>   Figures:
cat("    - predicted_probabilities_distribution.png\n")
#>     - predicted_probabilities_distribution.png
cat("    - marginal_effects.png\n")
#>     - marginal_effects.png
cat("    - mtm_effects.png\n")
#>     - mtm_effects.png

15 Key Findings Summary

################################################################################
# SECTION 17: KEY FINDINGS SUMMARY
################################################################################

cat("\n")
cat(strrep("=", 80), "\n", sep = "")
#> ================================================================================
cat("KEY FINDINGS SUMMARY\n")
#> KEY FINDINGS SUMMARY
cat(strrep("=", 80), "\n", sep = "")
#> ================================================================================
cat("\n1. CHOICE PATTERNS:\n")
#> 
#> 1. CHOICE PATTERNS:
cat(strrep("-", 80), "\n")
#> --------------------------------------------------------------------------------
cat("  - Most banks did not access Fed facilities during crisis\n")
#>   - Most banks did not access Fed facilities during crisis
cat("  - Among borrowers, BTFP was preferred over DW\n")
#>   - Among borrowers, BTFP was preferred over DW
cat("  - Small fraction used both facilities simultaneously\n")
#>   - Small fraction used both facilities simultaneously
cat("\n2. MTM LOSSES AND BTFP-ELIGIBLE SECURITIES:\n")
#> 
#> 2. MTM LOSSES AND BTFP-ELIGIBLE SECURITIES:
cat(strrep("-", 80), "\n")
#> --------------------------------------------------------------------------------
cat("  - Banks with greater MTM losses on BTFP-eligible securities\n")
#>   - Banks with greater MTM losses on BTFP-eligible securities
cat("    were MORE likely to use BTFP vs DW\n")
#>     were MORE likely to use BTFP vs DW
cat("  - This reflects the par value pricing subsidy\n")
#>   - This reflects the par value pricing subsidy
cat("  - MTM losses on other securities also matter but less so\n")
#>   - MTM losses on other securities also matter but less so
cat("\n3. UNINSURED DEPOSITS:\n")
#> 
#> 3. UNINSURED DEPOSITS:
cat(strrep("-", 80), "\n")
#> --------------------------------------------------------------------------------
cat("  - Higher uninsured deposit ratios increase use of Fed facilities\n")
#>   - Higher uninsured deposit ratios increase use of Fed facilities
cat("  - Run risk (interaction of uninsured deposits × MTM losses)\n")
#>   - Run risk (interaction of uninsured deposits × MTM losses)
cat("    is a key predictor of facility choice\n")
#>     is a key predictor of facility choice
cat("\n4. BANK SIZE:\n")
#> 
#> 4. BANK SIZE:
cat(strrep("-", 80), "\n")
#> --------------------------------------------------------------------------------
cat("  - Larger banks more likely to access facilities\n")
#>   - Larger banks more likely to access facilities
cat("  - Potentially reflects lower stigma and better relationships\n")
#>   - Potentially reflects lower stigma and better relationships
cat("  - Size effects differ between BTFP and DW\n")
#>   - Size effects differ between BTFP and DW
cat("\n5. FHLB AS ALTERNATIVE:\n")
#> 
#> 5. FHLB AS ALTERNATIVE:
cat(strrep("-", 80), "\n")
#> --------------------------------------------------------------------------------
cat("  - Banks with baseline FHLB relationships used it during crisis\n")
#>   - Banks with baseline FHLB relationships used it during crisis
cat("  - Positive correlation between FHLB and Fed facility use\n")
#>   - Positive correlation between FHLB and Fed facility use
cat("  - Suggests complementarity rather than substitution\n")
#>   - Suggests complementarity rather than substitution
cat("\n6. CAPITAL AND LIQUIDITY:\n")
#> 
#> 6. CAPITAL AND LIQUIDITY:
cat(strrep("-", 80), "\n")
#> --------------------------------------------------------------------------------
cat("  - Better-capitalized banks less likely to borrow\n")
#>   - Better-capitalized banks less likely to borrow
cat("  - Banks with more liquid assets less likely to use facilities\n")
#>   - Banks with more liquid assets less likely to use facilities
cat("  - Consistent with these being liquidity backstops\n")
#>   - Consistent with these being liquidity backstops
cat("\n✓ Key findings summarized\n")
#> 
#> ✓ Key findings summarized

16 Session Information

################################################################################
# SESSION INFO
################################################################################

cat("\n")
cat(strrep("=", 80), "\n", sep = "")
#> ================================================================================
cat("SESSION INFORMATION\n")
#> SESSION INFORMATION
cat(strrep("=", 80), "\n", sep = "")
#> ================================================================================
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] psych_2.5.6        moments_0.14.1     DescTools_0.99.60  modelsummary_2.3.0
#>  [5] stargazer_5.2.3    kableExtra_1.4.0   knitr_1.50         gridExtra_2.3     
#>  [9] patchwork_1.3.0    scales_1.3.0       ggthemes_5.1.0     ggplot2_3.5.1     
#> [13] broom_1.0.8        lmtest_0.9-40      zoo_1.8-12         sandwich_3.1-1    
#> [17] mvtnorm_1.3-3      mlogit_1.1-3       dfidx_0.2-0        nnet_7.3-19       
#> [21] fixest_0.12.1      readr_2.1.5        stringr_1.5.1      lubridate_1.9.4   
#> [25] 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] Rdpack_2.6          mnormt_2.1.1        gld_2.6.7          
#>  [4] readxl_1.4.3        rlang_1.1.2         magrittr_2.0.3     
#>  [7] dreamerr_1.4.0      e1071_1.7-14        compiler_4.3.1     
#> [10] systemfonts_1.0.6   vctrs_0.6.4         pkgconfig_2.0.3    
#> [13] crayon_1.5.2        fastmap_1.1.1       backports_1.4.1    
#> [16] labeling_0.4.3      utf8_1.2.4          rmarkdown_2.29     
#> [19] tzdb_0.4.0          haven_2.5.4         ragg_1.3.0         
#> [22] purrr_1.0.2         bit_4.0.5           xfun_0.54          
#> [25] cachem_1.0.8        jsonlite_1.8.7      tinytable_0.15.1   
#> [28] stringmagic_1.1.2   parallel_4.3.1      R6_2.5.1           
#> [31] bslib_0.5.1         tables_0.9.31       stringi_1.7.12     
#> [34] boot_1.3-28.1       jquerylib_0.1.4     cellranger_1.1.0   
#> [37] numDeriv_2016.8-1.1 Rcpp_1.0.12         parameters_0.28.3  
#> [40] Matrix_1.5-4.1      timechange_0.3.0    tidyselect_1.2.1   
#> [43] rstudioapi_0.16.0   yaml_2.3.7          lattice_0.21-8     
#> [46] bayestestR_0.17.0   withr_2.5.2         evaluate_0.23      
#> [49] proxy_0.4-27        xml2_1.3.6          pillar_1.9.0       
#> [52] checkmate_2.3.1     insight_1.4.3       generics_0.1.3     
#> [55] vroom_1.6.5         hms_1.1.3           munsell_0.5.0      
#> [58] rootSolve_1.8.2.4   class_7.3-22        glue_1.6.2         
#> [61] lmom_3.2            tools_4.3.1         forcats_1.0.0      
#> [64] Exact_3.3           fs_1.6.3            grid_4.3.1         
#> [67] rbibutils_2.2.16    datawizard_1.3.0    colorspace_2.1-0   
#> [70] nlme_3.1-162        performance_0.15.3  Formula_1.2-5      
#> [73] cli_3.6.1           textshaping_0.3.7   fansi_1.0.5        
#> [76] expm_1.0-0          viridisLite_0.4.2   svglite_2.1.3      
#> [79] gtable_0.3.6        sass_0.4.9          digest_0.6.33      
#> [82] farver_2.1.1        htmltools_0.5.7     lifecycle_1.0.4    
#> [85] httr_1.4.7          statmod_1.5.0       bit64_4.0.5        
#> [88] MASS_7.3-60

17 Appendix

17.1 Variable Definitions

17.1.1 Dependent Variables

Variable Description
Y / Y_factor Multinomial outcome: 0=No participation, 1=BTFP only, 2=DW only, 3=Both
btfp_user Binary: 1 if borrowed from BTFP during acute crisis
dw_user Binary: 1 if borrowed from DW during acute crisis
btfp_vs_dw Binary (borrowers only): 1 if used BTFP, 0 if DW only
fhlb_crisis Binary: 1 if FHLB borrowing > mean + 2SD

17.1.2 Key Independent Variables

Variable Description Source
mtm_btfp_pct MTM loss on BTFP-eligible securities / Assets (%) Python cleaning
mtm_other_pct MTM loss on other securities / Assets (%) Python cleaning
borrowing_subsidy MTM loss on eligible / Eligible securities (%) Calculated
uninsured_lev Uninsured deposits / Assets (%) Call Reports
pct_uninsured Uninsured deposits / Total deposits (%) Call Reports
run_risk_1 % Uninsured × % MTM loss Calculated
run_risk_1_dummy 1 if both > median Calculated
idcr_1 Coverage ratio (s=0.5) (MV assets - s×Uninsured - Insured) / Insured
idcr_2 Coverage ratio (s=1.0) (MV assets - s×Uninsured - Insured) / Insured
insolvency_1 Equity-based coverage (s=0.5) (Total assets - Total Liability - s×Uninsured×(Total assets/MV Asset) -1)) / Total assets
insolvency_2 Equity-based coverage (s=0.5) (Total assets - Total Liability - s×Uninsured×(Total assets/MV Asset) -1)) / Total assets

17.1.3 Control Variables

Variable Description
log_assets Log(Total assets)
pct_wholesale_liability (Fed funds + Repo + Other borrowed < 1 year) / Liabilities (%)
pct_liquidity_available (Cash + Reverse repo + Fed funds sold) / Assets (%)
securities_ratio Total securities / Assets (%)
loan_to_deposit Total loans / Total deposits (%)
book_equity_to_asset Book equity / Assets (%)
roa Net income (annualized) / Avg assets (%)
fhlb_to_asset FHLB advances / Assets (%)
delta_fhlb Change in FHLB from baseline to crisis

17.2 Model Specifications

17.2.1 Linear Probability Model (LPM)

\[P(BTFP = 1 | X) = X'\beta + \epsilon\]

17.2.2 Logit Model

\[P(BTFP = 1 | X) = \frac{e^{X'\beta}}{1 + e^{X'\beta}}\]

17.2.3 Multinomial Logit

\[P(Y = j | X) = \frac{e^{X'\beta_j}}{\sum_{k=0}^{3} e^{X'\beta_k}}\]

Base category: No participation (\(\beta_0 = 0\))

17.2.4 Trivariate Probit

Latent variables: \[BTFP^* = X'\beta_1 + \epsilon_1\] \[DW^* = X'\beta_2 + \epsilon_2\] \[FHLB^* = X'\beta_3 + \epsilon_3\]

Error correlation: \((\epsilon_1, \epsilon_2, \epsilon_3)' \sim N(0, \Sigma)\)

17.3 Sample Selection

  1. Baseline: 2022Q4 call reports
  2. Exclude: G-SIB banks, failed banks
  3. Require: omo_eligible > 0 (BTFP-eligible collateral)
  4. Crisis period: March 13 - April 30, 2023

End of Analysis