###
###
### PROGRAM NOTES
###
### The following R program will execute the Minneapolis Fed version of the CLASS model.  The default
### structure of the program (as it exists in the zip file) will generate projections using the 
### "extreme U-shaped recession" that is shown on the "COVID-19 stress test tool" web page listed below:
###
### https://www.minneapolisfed.org/banking/financial-studies-and-community-banking/run-your-own-stress-test-on-large-banks
###
###
### A more detailed description of the Minneapolis version of the model can be found in the article on stress
### testing for large banks in the COVID-19 pandemic available here:
###
### https://www.minneapolisfed.org/article/2020/capital-stress-testing-for-large-banks-in-light-of-the-covid-19-pandemic
###
###
### The model generates invidual projections for the 21 domestic banks described in the above article.  The firm-level
### projections are combined to form a representative projection for the industry.  Importantly, the individual
### firm balance sheets are held constant over the projection horizon.  As such, total assets, liabilities, 
### and risk-weighted assets remain fixed at their starting levels over the projection period.  In addition,
### the model assumes that dividends are held constant over the projection period at the 2019q4 levels and all
### share buybacks are eliminated in the projection period.  
###
### The default structure of the program and underlying data is based on a starting point of 2019q4.  Data from 
### this period are actual observations and all future observations are projections.  Note that the losses 
### used for operational risk, counterparty default, and the global market shock are based on the official 
### results from the DFAST2020 exercise.  The zip file contains five separate hypothetical scenarios that can be
### used to generate projections of capital under stress.  Alternatively, users can construct their own 
### scenarios and run those through the model if they choose.
###
###
### The program generates four CSV output files that are described below:
###
###   estimated_model_coefficients.csv
###     This file contains all of the necessary model coefficients as estimated by 
###     the Minneapolis Fed version of the CLASS model (the coefficient value, the standard error,
###     and the associated t-stat and p-value are included in the file).
###
###
###   aggregate_projections.csv
###     This file contains the combined firm-level projections that represent the industry aggregate
###     shown in the chart on the "COVID-19 stress test tool" web page.  The aggregate totals of 
###     common equity tier 1 capital (cet1), risk-weighted assets (rwa), average total assets (avg_assets),
###     the CET1 capital ratio (cap_ratio), the simple leverage ratio (lev_ratio) and the change in total
###     CET1 from the q0 starting point (cet1_drop) are included.  Values listed for the ratios are 
###     percentages (i.e., 7.35 = 7.35%) and the remaining items are all reported in billions of dollars.
###
###
###   firm_projections.csv
###     This file contains firm-level data for the same variables that are reported in the aggregate_projections.csv file.  
###     In addition, the file contains a new variable (breach) which takes on the value of 1 for quarters
###     in which a firm's CET1 capital level below the "capital threshold" (defined as the statutory minimum of 4.5% of RWA plus the
###     applicable GSIB surcharge).
###
###   firm_breaches.csv
###       This is a summary file which simply reports whether or not a firm had their projected CET1 level fall
###       below the firm-specific "capital threshold" level.  A value of 1 indicates that such a breach occured 
###       at some point during the projection horizon and a value of 0 indicates that the capital level stayed above
###       the appropriate threshold.
###
###
###
###   Instructions for running the Minneapolis Fed version of the CLASS model
###   1. The model is programmed in R and can be run using either R or Rstudio.
###
###   2. Extract the full set of material contained in the zip file into a 
###      a folder (for example, C:\mpls_stress_test).
###
###   3. The folder should contain the following items once the extraction is complete:
###       
###       C:\mpls_stress_test\macro_data_hist.csv
###           The historical macroeconomic data used to estimate the model.
###
###       C:\mpls_stress_test\macro_data_proj.csv
###           The hypothetical macroeconomic data used to generate the projections.  This
###           file contains 5 separate scenarios:  the DFAST2020 baseline scenario (fed_baseline),
###           the DFAST2020 severely adverse scenario (fed_severe), and the three U-shaped COVID
###           recessions (u_normal, u_severe, and u_extreme).
###
###       C:\mpls_stress_test\y9_bhc_data.csv
###           The historical firm-level balance sheet and income statement data (derived from the Y-9C report).
###
###       C:\mpls_stress_test\mpls_stress_test_model.R 
###           The R program that estimates the model and generates the projections associated with the 
###           hypothetical scenario data.
###
###       C:\mpls_stress_test\read_me.TXT
###           A text file version of the "program notes".
###
###       C:\mpls_stress_test\output
###           The folder that contains the 4 CSV output files generated by the model (and described above).
###
###   4. Open the R program once the zip file has been successfully extracted.  Adjust the following variables
###      in Part 2 (specify starting parameters) of the program:
###       
###       wksp      - the name of the folder where the zip file was extracted (i.e., "c:\\mpls_stress_test\\)
### 
###       scen_nm   - the name of the hypothetical scenario to use for generating the projections.  It
###                   needs to be set equal to one of the following values:  fed_baseline, fed_severe,
###                   u_normal, u_severe, or u_extreme.
###
###   5. Run the program from within R or RStudio
###
###






#############################################################################
#############################################################################
#############################################################################

### Part 1 -- remove all objects from the current R environment and load necessary packages

rm(list=ls())
  

require(dplyr)
require(haven)
require(plm)
require(xlsx)
require(zoo)
      
      
      
#############################################################################
#############################################################################
#############################################################################
      
### Part 2 -- specify starting parameters

wksp              <- "c:\\mpls_stress_test\\"       

scen_nm           <- "u_extreme"



fn_class_data     <- paste(wksp,"y9c_bhc_data.csv",sep="")
fn_macro_hist     <- paste(wksp,"macro_data_hist.csv",sep="")
fn_macro_proj     <- paste(wksp,"macro_data_proj.csv",sep="")
      
  

### adjusting these parameters will generate results for alternative 
### experiments (they will not match the ones shown in the Mpls Fed stress test tool)

ops_year    <- 2020
gms_year    <- 2020

do_special  <- 1
    
exclude_div  <- 0

date_q0     <- 20191231
      
true_horiz  <- 8
    
tax_rate    <- 0.21
      
ops_2019    <- 0.917714957 * 123 * 1e6
ops_2020    <- 0.874301107 * 144 * 1e6

      
gsib <- data.frame(
  entity=c(
  1039502,
  1951350,
  1073757,
  1120754,
  2380443,
  1111435,
  2162966,
  3587146
  ),
  charge=c(
  2.5,
  2,
  1.5,
  1.5,
  1.5,
  1,
  1,
  1
  ))      



gms_2019 <- data.frame(
  entity=c(
  1039502,
  1951350,
  2380443,
  2162966,
  1073757,
  1120754,
  1111435,
  3587146),
  loss_gms=c(
  21200000,
  15700000,
  13600000,
  9800000,
  8900000,
  8500000,
  1200000,
  1100000
  ))
      
      
gms_2020 <- data.frame(
  entity=c(
  1039502,
  1951350,
  2380443,
  2162966,
  1073757,
  1120754,
  1111435,
  3587146),
  loss_gms=c(
  21800000,
  6000000,
  18400000,
  10100000,
  10500000,
  900000,
  600000,
  800000
  ))
      
      
if (ops_year==2019) {
  official_ops <- ops_2019
} else if (ops_year==2020) {  
  official_ops <- ops_2020
}
      
      
if (gms_year==2019) {
  official_gms <- gms_2019
} else if (gms_year==2020) {  
  official_gms <- gms_2020
}
      
      
      
      
banks <- data.frame(
  entity=c(
  1037003, 
  1039502, 
  1068025, 
  1068191, 
  1069778, 
  1070345, 
  1073757, 
  1074156, 
  1111435, 
  1119794, 
  1120754, 
  1199611, 
  1275216, 
  1562859, 
  1951350, 
  2162966, 
  2277860, 
  2380443, 
  3242838, 
  3587146, 
  3846375),
  tkr=c(
  "MT",
  "JPM",
  "KEY",
  "HBAN",
  "PNC",
  "FITB",
  "BAC",
  "TFC",
  "STT",
  "USB",
  "WFC",
  "NTRS",
  "AXP",
  "ALLY",
  "C",
  "MS",
  "COF",
  "GS",
  "RF",
  "BK",
  "DFS"
  ))
      
      
      
#############################################################################
#############################################################################
#############################################################################
      
### Part 3 -- read the necessary input datasets (bank level data, historical macro data, and projected macro data)
      
      
bb      <- read.csv(fn_class_data)
bb$date <- bb$dt
bb$dt   <- NULL
      
      
      
mm                <- read.csv(fn_macro_hist)
mm                <- mm[order(mm$date),]
mm                <- mm %>% arrange(date) %>% mutate(d10y = dplyr::lag(cpr_t10y,1)) 
mm$d10y           <- mm$cpr_t10y - mm$d10y
mm$dbspread_pos   <- ifelse(mm$dbspread > 0, mm$dbspread, 0)
      
      
      
pp              <- read.csv(fn_macro_proj,stringsAsFactors = FALSE)
pp              <- pp[which(pp$type==scen_nm),]
pp              <- pp[order(pp$date),]
pp$type         <- NULL
pp$dt           <- NULL
pp$id           <- 1
pp              <- pp[order(pp$id,pp$date),]
      
      
pp              <- pdata.frame(pp,index=c("id","date"))
pp$spread       <- pp$cpr_t10y - pp$cpr_t3m
pp$d10y         <- pp$cpr_t10y - plm::lag(pp$cpr_t10y,1,shift="row")
pp$stockgrowth  <- log(pp$djtsm / plm::lag(pp$djtsm,1,shift="row"))*100
pp$unchange     <- (pp$cpr_ur - plm::lag(pp$cpr_ur,1,shift="row"))*4
      
pp$yyhpi        <- log(pp$hpi / plm::lag(pp$hpi,4,shift="row"))*100
pp$yyhpi_zero   <- ifelse(pp$yyhpi <= 0, pp$yyhpi, 0)
pp$yycppi       <- log(pp$cppi / plm::lag(pp$cppi,4,shift="row"))*100
pp$yycppi_zero  <- ifelse(pp$yycppi <= 0, pp$yycppi, 0)
pp$bspread      <- pp$bbb_corp - pp$cpr_t10y
pp$dbspread     <- pp$bspread - plm::lag(pp$bspread,1,shift="row")
pp$dbspread_pos <- ifelse(pp$dbspread >= 0, pp$dbspread, 0)
      
      
      
pp      <- pp[,names(mm)]
pp$ddt  <- as.numeric(substring(pp$date,1,4)) + as.numeric(substring(pp$date,5,6))/12
pp      <- pp[which(pp$ddt > as.numeric(substring(date_q0,1,4))+as.numeric(substring(date_q0,5,6))/12),]


      
#############################################################################
#############################################################################
#############################################################################
      
### Part 4 -- assemble the data for the regression estimations
      
      
      
###
###   merge the scenario data onto the BHC data frame and order it by [entity,date]
###
      
      
bb      <- merge(bb,mm,by="date",all.x=T)
bb      <- bb[order(bb$entity,bb$date),]
      
      
###
### date variable stuff
###
      
bb$dt   <- paste(substring(bb$date,5,6),"/",substring(bb$date,7,8),"/",substring(bb$date,1,4),sep="")
bb$dt   <- as.Date(bb$dt,"%m/%d/%Y")
      
bb$ddt  <- as.numeric(substring(bb$date,1,4)) + as.numeric(substring(bb$date,5,6))/12
      
      
      
bb_ind  <- bb[which(bb$entity == 9999999),]
bb      <- bb[which(bb$entity != 9999999),]
      
      
### 
###   create the necessary share variables (used as weights in regressions)
###
      
for (jj in 1:2) {
      
  if (jj==1) {
    ff      <- bb  
  } else {
    ff      <- bb_ind
  }
        
  ff      <- ff %>% group_by(date) %>% mutate(total_assets = sum(cl_assets,na.rm=T)) 
  ff      <- ff %>% group_by(date) %>% mutate(total_intearn_assets = sum(cl_intearn_ass,na.rm=T))
  ff      <- ff %>% group_by(date) %>%  mutate(total_afssec = sum(cl_afssec,na.rm=T))
        
  ff      <- ff %>% mutate(shr_assets = 100*(cl_assets / total_assets))
  ff      <- ff %>% mutate(shr_intearn_assets = 100*(cl_intearn_ass / total_intearn_assets))
  ff      <- ff %>% mutate(shr_afssec = 100*(cl_afssec / total_afssec))
      
  if (jj==1) {
    bb      <- ff  
  } else {
    bb_ind  <- ff
  }
        
}
      
bb        <- rbind(bb,bb_ind)
      
      
      
### 
###   create other required variables
###
      
bb$ddt    <- floor(bb$date/1e4) + as.numeric(substring(bb$date,5,6))/12
bb$time   <- bb$ddt - min(bb$ddt)
      
      
bb$unsafe_dbspread      <- bb$unsafe_ratio * bb$dbspread
bb$unsafe_dbspread_pos  <- ifelse(bb$dbspread > 0, bb$unsafe_dbspread, 0)
      
      
        
### 
###   create the necessary lagged variables
###
      
bb      <- bb[order(bb$entity,bb$date),]
      
pbb     <- pdata.frame(bb,index=c("entity","dt"))
      
      
pbb$cl_nim_lag            <- plm::lag(pbb$cl_nim,1)
pbb$cl_nintrat_lag        <- plm::lag(pbb$cl_nintrat,1)
pbb$cl_nie_comp_lag       <- plm::lag(pbb$cl_nie_comp,1)
pbb$cl_nie_fass_lag       <- plm::lag(pbb$cl_nie_fass,1)
pbb$cl_nie_allother_lag   <- plm::lag(pbb$cl_nie_allother,1)
pbb$cl_afs_return_lag     <- plm::lag(pbb$cl_afs_return,1)
pbb$cl_tradrat_lag        <- plm::lag(pbb$cl_tradrat,1)
      
      
pbb$cl_nco_firstlien_lag    <- plm::lag(pbb$cl_nco_firstlien,1)
pbb$cl_nco_jrlien_lag       <- plm::lag(pbb$cl_nco_jrlien,1)
pbb$cl_nco_heloc_lag        <- plm::lag(pbb$cl_nco_heloc,1)
pbb$cl_nco_const_lag        <- plm::lag(pbb$cl_nco_const,1)
pbb$cl_nco_multi_lag        <- plm::lag(pbb$cl_nco_multi,1)
pbb$cl_nco_nfnr_lag         <- plm::lag(pbb$cl_nco_nfnr,1)
pbb$cl_nco_othre_lag        <- plm::lag(pbb$cl_nco_othre,1)
pbb$cl_nco_ci_lag           <- plm::lag(pbb$cl_nco_ci,1)
pbb$cl_nco_cc_lag           <- plm::lag(pbb$cl_nco_cc,1)
pbb$cl_nco_othcons_lag      <- plm::lag(pbb$cl_nco_othcons,1)
pbb$cl_nco_agr_lag          <- plm::lag(pbb$cl_nco_agr,1)
pbb$cl_nco_lease_lag        <- plm::lag(pbb$cl_nco_lease,1)
pbb$cl_nco_oth_lag          <- plm::lag(pbb$cl_nco_oth,1)
      
      
      
      
##########################################################################
##########################################################################
##########################################################################
      
# Part 5 -- estimate the individual models and assemble the coefficients into a single dataframe
      
      
cc      <- data.frame(pbb)
      
xnm           <- "nim"
      
xstr          <-  paste("cl_", xnm, " ~ cl_", xnm, "_lag + spread + cpr_t3m + rat_rre + rat_cre + rat_ci + rat_cc + rat_trad + rat_inv + shr_assets + time",sep="")
xwt           <- "shr_intearn_assets"
xxr           <- cc[which(cc$ddt >= 1996.25 & cc$date <= date_q0 & cc$entity !=  9999999),]
      
      
xxr_eq    	  <- lm(formula(xstr),xxr,weights=xxr[[xwt]],na.action=na.exclude)
      
xxr_sum       <- summary(xxr_eq)
      
xxr_stats     <- data.frame(mdl=xnm, 
                            nobs=(xxr_sum$df[1]+xxr_sum$df[2]),
                            fstat=xxr_sum$fstatistic,
                            adjr2=xxr_sum$adj.r.squared)
      
      
xxr_coef      <- data.frame(summary(xxr_eq)$coef,stringsAsFactors=F)
      
xxr_mat       <- data.frame(mdl=xnm,
                            vars=attributes(xxr_coef)$row.names,
                            coef=xxr_coef$Estimate,
                            stderr=xxr_coef$Std..Error,
                            tstat=xxr_coef$t.value,
                            pvalue=xxr_coef$Pr...t..,
                            stringsAsFactors=F)
      
      
coef_matrix     <- xxr_mat
      
      
      
##########################################################################
      
xnm           <- "nintrat"              
xstr          <-  paste("cl_", xnm, " ~ cl_", xnm, "_lag + stockgrowth + rat_rre + rat_cre + rat_ci + rat_cc + rat_trad + rat_inv + shr_assets",sep="")
      
xwt           <- "shr_assets"
xxr           <- cc[which(cc$ddt >= 1996.25 & cc$date <= date_q0 & cc$entity !=  9999999),]
      
      
xxr_eq    	  <- lm(formula(xstr),xxr,weights=xxr[[xwt]],na.action=na.exclude)
      
xxr_sum       <- summary(xxr_eq)
      
xxr_stats     <- data.frame(mdl=xnm, 
                            nobs=(xxr_sum$df[1]+xxr_sum$df[2]),
                            fstat=xxr_sum$fstatistic,
                            adjr2=xxr_sum$adj.r.squared)
      
xxr_coef      <- data.frame(summary(xxr_eq)$coef,stringsAsFactors=F)
      
xxr_mat       <- data.frame(mdl=xnm,
                            vars=attributes(xxr_coef)$row.names,
                            coef=xxr_coef$Estimate,
                            stderr=xxr_coef$Std..Error,
                            tstat=xxr_coef$t.value,
                            pvalue=xxr_coef$Pr...t..,
                            stringsAsFactors=F)
      
coef_matrix     <- rbind(coef_matrix,xxr_mat)
      
      
      
##########################################################################
      
xnm           <- "nie_comp"              
xstr          <-  paste("cl_", xnm, " ~ cl_", xnm, "_lag + stockgrowth + rat_rre + rat_cre + rat_ci + rat_cc + rat_trad + rat_inv + shr_assets",sep="")
      
xwt           <- "shr_assets"
xxr           <- cc[which(cc$ddt >= 1996.25 & cc$date <= date_q0 & cc$entity !=  9999999),]
      
      
xxr_eq    	  <- lm(formula(xstr),xxr,weights=xxr[[xwt]],na.action=na.exclude)
      
xxr_sum       <- summary(xxr_eq)
      
xxr_stats     <- data.frame(mdl=xnm, 
                            nobs=(xxr_sum$df[1]+xxr_sum$df[2]),
                            fstat=xxr_sum$fstatistic,
                            adjr2=xxr_sum$adj.r.squared)
      
xxr_coef      <- data.frame(summary(xxr_eq)$coef,stringsAsFactors=F)
      
xxr_mat       <- data.frame(mdl=xnm,
                            vars=attributes(xxr_coef)$row.names,
                            coef=xxr_coef$Estimate,
                            stderr=xxr_coef$Std..Error,
                            tstat=xxr_coef$t.value,
                            pvalue=xxr_coef$Pr...t..,
                            stringsAsFactors=F)
      
coef_matrix     <- rbind(coef_matrix,xxr_mat)
      
      
##########################################################################
      
xnm           <- "nie_fass"              
xstr          <-  paste("cl_", xnm, " ~ cl_", xnm, "_lag + growth + rat_rre + rat_cre + rat_ci + rat_cc + rat_trad + rat_inv + shr_assets + time",sep="")
      
xwt           <- "shr_assets"
xxr           <- cc[which(cc$ddt >= 1996.25 & cc$date <= date_q0 & cc$entity !=  9999999),]
      
      
xxr_eq    	  <- lm(formula(xstr),xxr,weights=xxr[[xwt]],na.action=na.exclude)
      
xxr_sum       <- summary(xxr_eq)
      
xxr_stats     <- data.frame(mdl=xnm, 
                            nobs=(xxr_sum$df[1]+xxr_sum$df[2]),
                            fstat=xxr_sum$fstatistic,
                            adjr2=xxr_sum$adj.r.squared)
      
xxr_coef      <- data.frame(summary(xxr_eq)$coef,stringsAsFactors=F)
      
xxr_mat       <- data.frame(mdl=xnm,
                            vars=attributes(xxr_coef)$row.names,
                            coef=xxr_coef$Estimate,
                            stderr=xxr_coef$Std..Error,
                            tstat=xxr_coef$t.value,
                            pvalue=xxr_coef$Pr...t..,
                            stringsAsFactors=F)
      
coef_matrix     <- rbind(coef_matrix,xxr_mat)
      
      
##########################################################################
      
xnm           <- "nie_allother"              
xstr          <-  paste("cl_", xnm, " ~ cl_", xnm, "_lag + dbspread + rat_rre + rat_cre + rat_ci + rat_cc + rat_trad + rat_inv + shr_assets",sep="")
      
xwt           <- "shr_assets"
xxr           <- cc[which(cc$ddt >= 1996.25 & cc$date <= date_q0 & cc$entity !=  9999999),]
      
      
xxr_eq    	  <- lm(formula(xstr),xxr,weights=xxr[[xwt]],na.action=na.exclude)
      
xxr_sum       <- summary(xxr_eq)
      
xxr_stats     <- data.frame(mdl=xnm, 
                            nobs=(xxr_sum$df[1]+xxr_sum$df[2]),
                            fstat=xxr_sum$fstatistic,
                            adjr2=xxr_sum$adj.r.squared)
      
xxr_coef      <- data.frame(summary(xxr_eq)$coef,stringsAsFactors=F)
      
xxr_mat       <- data.frame(mdl=xnm,
                            vars=attributes(xxr_coef)$row.names,
                            coef=xxr_coef$Estimate,
                            stderr=xxr_coef$Std..Error,
                            tstat=xxr_coef$t.value,
                            pvalue=xxr_coef$Pr...t..,
                            stringsAsFactors=F)
      
coef_matrix     <- rbind(coef_matrix,xxr_mat)
      
      
      
#############################################################################
      
xnm           <- "afs_return"              
xstr          <-  paste("cl_", xnm, " ~ cl_", xnm, "_lag + d10y + unsafe_dbspread_pos",sep="")
      
xwt           <- "shr_afssec"
xxr           <- cc[which(cc$ddt >= 1997.25 & cc$date <= date_q0 & cc$entity !=  9999999),]
      
      
xxr_eq    	  <- lm(formula(xstr),xxr,weights=xxr[[xwt]],na.action=na.exclude)
      
xxr_sum       <- summary(xxr_eq)
      
xxr_stats     <- data.frame(mdl=xnm, 
                            nobs=(xxr_sum$df[1]+xxr_sum$df[2]),
                            fstat=xxr_sum$fstatistic,
                            adjr2=xxr_sum$adj.r.squared)
      
xxr_coef      <- data.frame(summary(xxr_eq)$coef,stringsAsFactors=F)
      
xxr_mat       <- data.frame(mdl=xnm,
                            vars=attributes(xxr_coef)$row.names,
                            coef=xxr_coef$Estimate,
                            stderr=xxr_coef$Std..Error,
                            tstat=xxr_coef$t.value,
                            pvalue=xxr_coef$Pr...t..,
                            stringsAsFactors=F)
      
coef_matrix     <- rbind(coef_matrix,xxr_mat)
      
      
      
#############################################################################
      
xnm           <- "tradrat"              
xstr          <-  paste("cl_", xnm, " ~ cl_", xnm, "_lag + dbspread + dbspread_pos",sep="")
      
xxr           <- cc[which(cc$ddt >= 1997 & cc$date <= date_q0 & cc$entity == 9999999),]
      
      
xxr_eq    	  <- lm(formula(xstr),xxr,na.action=na.exclude)
      
xxr_sum       <- summary(xxr_eq)
      
xxr_stats     <- data.frame(mdl=xnm, 
                            nobs=(xxr_sum$df[1]+xxr_sum$df[2]),
                            fstat=xxr_sum$fstatistic,
                            adjr2=xxr_sum$adj.r.squared)
      
xxr_coef      <- data.frame(summary(xxr_eq)$coef,stringsAsFactors=F)
      
xxr_mat       <- data.frame(mdl=xnm,
                            vars=attributes(xxr_coef)$row.names,
                            coef=xxr_coef$Estimate,
                            stderr=xxr_coef$Std..Error,
                            tstat=xxr_coef$t.value,
                            pvalue=xxr_coef$Pr...t..,
                            stringsAsFactors=F)
      
coef_matrix     <- rbind(coef_matrix,xxr_mat)
      
      
      
#############################################################################
      
      
for (k in c("nco_firstlien","nco_jrlien","nco_heloc")) {
      
  xnm           <- k
  xstr          <-  paste("cl_", xnm, " ~ cl_", xnm, "_lag + yyhpi + yyhpi_zero",sep="")
      
  xxr           <- cc[which(cc$ddt >= 1996.25 & cc$date <= date_q0 & cc$entity == 9999999),]
      
  xxr_eq    	  <- lm(formula(xstr),xxr,na.action=na.exclude)
      
  xxr_sum       <- summary(xxr_eq)
      
  xxr_stats     <- data.frame(mdl=xnm, 
                              nobs=(xxr_sum$df[1]+xxr_sum$df[2]),
                              fstat=xxr_sum$fstatistic,
                              adjr2=xxr_sum$adj.r.squared)
      
  xxr_coef      <- data.frame(summary(xxr_eq)$coef,stringsAsFactors=F)
      
  xxr_mat       <- data.frame(mdl=xnm,
                              vars=attributes(xxr_coef)$row.names,
                              coef=xxr_coef$Estimate,
                              stderr=xxr_coef$Std..Error,
                              tstat=xxr_coef$t.value,
                              pvalue=xxr_coef$Pr...t..,
                              stringsAsFactors=F)
      
  coef_matrix     <- rbind(coef_matrix,xxr_mat)
}
      
      
      
#############################################################################
      
      
for (k in c("nco_const","nco_multi","nco_nfnr")) {
        
  xnm           <- k
  xstr          <-  paste("cl_", xnm, " ~ cl_", xnm, "_lag + yycppi_zero",sep="")
        
  xxr           <- cc[which(cc$ddt >= 1996.26 & cc$date <= date_q0 & cc$entity == 9999999),]
        
  xxr_eq    	  <- lm(formula(xstr),xxr,na.action=na.exclude)
        
  xxr_sum       <- summary(xxr_eq)
        
  xxr_stats     <- data.frame(mdl=xnm, 
                              nobs=(xxr_sum$df[1]+xxr_sum$df[2]),
                              fstat=xxr_sum$fstatistic,
                              adjr2=xxr_sum$adj.r.squared)
        
  xxr_coef      <- data.frame(summary(xxr_eq)$coef,stringsAsFactors=F)
        
  xxr_mat       <- data.frame(mdl=xnm,
                              vars=attributes(xxr_coef)$row.names,
                              coef=xxr_coef$Estimate,
                              stderr=xxr_coef$Std..Error,
                              tstat=xxr_coef$t.value,
                              pvalue=xxr_coef$Pr...t..,
                              stringsAsFactors=F)
        
  coef_matrix     <- rbind(coef_matrix,xxr_mat)
}
      
      
#############################################################################
      
      
for (k in c("nco_ci","nco_cc","nco_agr","nco_lease","nco_oth")) {
        
  xnm           <- k
  xstr          <-  paste("cl_", xnm, " ~ cl_", xnm, "_lag + unchange",sep="")
        
  xxr           <- cc[which(cc$ddt >= 1996.25 & cc$date <= date_q0 & cc$entity == 9999999),]
        
  xxr_eq    	  <- lm(formula(xstr),xxr,na.action=na.exclude)
        
  xxr_sum       <- summary(xxr_eq)
        
  xxr_stats     <- data.frame(mdl=xnm, 
                              nobs=(xxr_sum$df[1]+xxr_sum$df[2]),
                              fstat=xxr_sum$fstatistic,
                              adjr2=xxr_sum$adj.r.squared)
        
  xxr_coef      <- data.frame(summary(xxr_eq)$coef,stringsAsFactors=F)
        
  xxr_mat       <- data.frame(mdl=xnm,
                              vars=attributes(xxr_coef)$row.names,
                              coef=xxr_coef$Estimate,
                              stderr=xxr_coef$Std..Error,
                              tstat=xxr_coef$t.value,
                              pvalue=xxr_coef$Pr...t..,
                              stringsAsFactors=F)
        
  coef_matrix     <- rbind(coef_matrix,xxr_mat)
}
      
      
      
      
#############################################################################
      
      
xnm           <- "nco_othre"
xstr          <-  paste("cl_", xnm, " ~ cl_", xnm, "_lag + yycppi + yycppi_zero",sep="")
        
xxr           <- cc[which(cc$ddt >= 1996.25 & cc$date <= date_q0 & cc$entity == 9999999),]
        
xxr_eq    	  <- lm(formula(xstr),xxr,na.action=na.exclude)
        
xxr_sum       <- summary(xxr_eq)
        
xxr_stats     <- data.frame(mdl=xnm, 
                            nobs=(xxr_sum$df[1]+xxr_sum$df[2]),
                            fstat=xxr_sum$fstatistic,
                            adjr2=xxr_sum$adj.r.squared)
        
xxr_coef      <- data.frame(summary(xxr_eq)$coef,stringsAsFactors=F)
        
xxr_mat       <- data.frame(mdl=xnm,
                            vars=attributes(xxr_coef)$row.names,
                            coef=xxr_coef$Estimate,
                            stderr=xxr_coef$Std..Error,
                            tstat=xxr_coef$t.value,
                            pvalue=xxr_coef$Pr...t..,
                            stringsAsFactors=F)
        
coef_matrix     <- rbind(coef_matrix,xxr_mat)
      
      
      
#############################################################################
      
      
xnm           <- "nco_othcons"
xstr          <-  paste("cl_", xnm, " ~ cl_", xnm, "_lag + unchange + time",sep="")
      
xxr           <- cc[which(cc$ddt >= 1996.25 & cc$date <= date_q0 & cc$entity == 9999999),]
      
xxr_eq    	  <- lm(formula(xstr),xxr,na.action=na.exclude)
      
xxr_sum       <- summary(xxr_eq)
      
xxr_stats     <- data.frame(mdl=xnm, 
                            nobs=(xxr_sum$df[1]+xxr_sum$df[2]),
                            fstat=xxr_sum$fstatistic,
                            adjr2=xxr_sum$adj.r.squared)
      
xxr_coef      <- data.frame(summary(xxr_eq)$coef,stringsAsFactors=F)
      
xxr_mat       <- data.frame(mdl=xnm,
                            vars=attributes(xxr_coef)$row.names,
                            coef=xxr_coef$Estimate,
                            stderr=xxr_coef$Std..Error,
                            tstat=xxr_coef$t.value,
                            pvalue=xxr_coef$Pr...t..,
                            stringsAsFactors=F)
      
coef_matrix     <- rbind(coef_matrix,xxr_mat)
      
coef_matrix$vars[which(coef_matrix$vars=="(Intercept)")] <- "constant"
      
      
      
      
##########################################################################
##########################################################################
##########################################################################
      
# Part 6 -- create additional variables needed for the bank-level projections and attach the
#           ops risk and global market shock losses
      
      
      
dd        <- cc
dd        <- dd[which(dd$entity != 9999999),]
dd        <- dd[which(dd$date == date_q0),]
      
dd        <- merge(banks,dd,by="entity",all.x=T)
      
      
  
if(exclude_div==1) {
  dd$cl_dividend      <- 0
  } else {
    dd$cl_dividend        <- dd$qprefdividend + dd$qcommdividend
}
      
      
dd$qnonint_inc_nt     <- dd$qnonint_inc - dd$qtradrev_inc
dd$pretax             <- dd$cl_ppnr - dd$qllprov + dd$qafssecur_inc + dd$qhtmsecur_inc
      
dd$qnetxoff           <- dd$qnetxoff_total
dd$ln_tot             <- dd$cl_loans
dd$cet1_ratio         <- dd$cet1 / dd$rwa * 100
      
      
if(do_special==1) {
  dd$loss_ops           <- dd$cl_assets / sum(dd$cl_assets) * official_ops / 9
  } else {
  dd$loss_ops           <- 0
}
      
      
      
dd      <- merge(dd,official_gms,by="entity",all.x=T)
      
if (do_special == 1) {
  dd$loss_gms[is.na(dd$loss_gms)] <- 0
  } else {
    dd$loss_gms <- 0
}
      
      
dd        <- merge(dd,gsib,by="entity",all.x=T)
dd$charge[is.na(dd$charge)] <- 0
      
dd$req    <- ifelse(dd$charge == 0, 0.045, 0.045 + dd$charge/100)
      
      
      
##########################################################################
##########################################################################
##########################################################################
      
# Part 7 -- generate all of the bank-level projections using the estimated coefficients
      
      
dd0      <- dd
      
      
      
dd0      <- dd0[,c("entity","ddt","time","rat_cc","rat_ci","rat_cre","rat_inv","rat_rre","rat_trad",
                  "shr_afssec","shr_assets","shr_intearn_assets","unsafe_ratio",
                  "cl_nim", "cl_tradrat", "cl_nintrat", "cl_nie_comp", "cl_nie_fass",
                  "cl_nie_allother",  "cl_afs_return",
                  "cl_nco_heloc", "cl_nco_firstlien", "cl_nco_jrlien", "cl_nco_const", "cl_nco_multi", "cl_nco_nfnr", "cl_nco_othre",
                  "cl_nco_agr", "cl_nco_oth","cl_nco_ci","cl_nco_cc","cl_nco_othcons","cl_nco_lease",
                  "cl_intearn_ass","cl_assets","cl_afssec", "cl_trading_ass", "avg_assets",
                  "ln_firstlien", "ln_jrlien", "ln_heloc", "ln_const", "ln_multi", "ln_nfnr",
                  "ln_agr","ln_cc","ln_ci", "ln_othcons", "ln_oth", "ln_othre", "ln_lease",
                  "qint_inc_net","qnonint_inc","qtradrev_inc","qnonint_exp_comp","qnonint_exp_fass","qnonint_exp_allother",
                  "qafssecur_inc","qhtmsecur_inc",
                  "qnetxoff_firstlien","qnetxoff_jrlien","qnetxoff_heloc","qnetxoff_const","qnetxoff_multi",
                  "qnetxoff_nfnr","qnetxoff_agr","qnetxoff_cc","qnetxoff_ci","qnetxoff_othcons","qnetxoff_oth",
                  "qnetxoff_othre","qnetxoff_lease",
                  "llres", "loss_ops", "loss_gms",
                  "qminor_int", "cl_dividend", "equity", "cet1", "rwa","req"
                  )]
      
                             
ee      <- dd0
      
idt_q0 <- min(ee$ddt)
      
      
      
### generate 13 quarterly projections      
for (tt in 1:13) {
      
  ff              <- ee[which(ee$ddt==idt_q0 + (tt-1)*0.25),]
  ff_names        <- names(ff)
          
  ff$ddt          <- ff$ddt + 0.25
  ff$constant     <- 1
                
          
  keydt     <- min(ff$ddt)
      
  ppt       <- pp[which(pp$ddt==keydt),]
      
  ff        <- merge(ff,ppt,by="ddt",all.x=T)
      
  ff$unsafe_dbspread      <- ff$unsafe_ratio * ff$dbspread
  ff$unsafe_dbspread_pos  <- ifelse(ff$dbspread > 0, ff$unsafe_dbspread, 0)
        
        
  for (m in c("nim","nintrat","tradrat","afs_return", 
              "nie_comp","nie_fass","nie_allother",
              "nco_firstlien","nco_jrlien","nco_heloc","nco_const","nco_multi",
              "nco_nfnr","nco_othre","nco_ci","nco_cc","nco_othcons",
              "nco_agr","nco_lease","nco_oth")) {
                    
    ff[[paste("cl_",m,"_lag",sep="")]]  <- ff[[paste("cl_",m,sep="")]]
    beta        <- coef_matrix[which(coef_matrix==m),c("vars","coef")]
    xx          <- ff[,beta$vars]
    yy          <- as.matrix(xx) %*% as.matrix(beta$coef)
       
      
    ff[[paste("cl_",m,sep="")]]           <- yy
    ff[[paste("cl_",m,"_lag",sep="")]]    <- NULL
      
  }
      
        
      
      
  ff                  <- ff[,ff_names]
        
        
  ff$qint_inc_net           <- (ff$cl_nim / 400) * ff$cl_intearn_ass
  ff$qnonint_inc            <- (ff$cl_nintrat / 400) * ff$cl_assets  
  ff$qtradrev_inc           <- (ff$cl_tradrat / 400) * ff$cl_trading_ass
  ff$qnonint_exp_comp       <- (ff$cl_nie_comp  / 400) * ff$cl_assets
  ff$qnonint_exp_fass       <- (ff$cl_nie_fass / 400) * ff$cl_assets
  ff$qnonint_exp_allother   <- (ff$cl_nie_allother / 400) * ff$cl_assets
        
  ff$qafssecur_inc          <- (ff$cl_afs_return / 400) * ff$cl_afssec
  ff$qhtmsecur_inc          <- 0
        
  ff$qnetxoff_agr           <- (ff$cl_nco_agr / 400) * ff$ln_agr
  ff$qnetxoff_cc            <- (ff$cl_nco_cc / 400) * ff$ln_cc
  ff$qnetxoff_ci            <- (ff$cl_nco_ci / 400) * ff$ln_ci
  ff$qnetxoff_const         <- (ff$cl_nco_const / 400) * ff$ln_const
  ff$qnetxoff_firstlien     <- (ff$cl_nco_firstlien / 400) * ff$ln_firstlien
  ff$qnetxoff_heloc         <- (ff$cl_nco_heloc / 400) * ff$ln_heloc
  ff$qnetxoff_jrlien        <- (ff$cl_nco_jrlien / 400) * ff$ln_jrlien
  ff$qnetxoff_lease         <- (ff$cl_nco_lease / 400) * ff$ln_lease
  ff$qnetxoff_multi         <- (ff$cl_nco_multi / 400) * ff$ln_multi
  ff$qnetxoff_nfnr          <- (ff$cl_nco_nfnr / 400) * ff$ln_nfnr
  ff$qnetxoff_oth           <- (ff$cl_nco_oth / 400) * ff$ln_oth
  ff$qnetxoff_othcons       <- (ff$cl_nco_othcons / 400) * ff$ln_othcons
  ff$qnetxoff_othre         <- (ff$cl_nco_othre / 400) * ff$ln_othre
      
      
  ee          <- rbind(ee,ff) 
          
}
      
      
##########################################################################
##########################################################################
##########################################################################
      
# Part 8 -- generate the provision expense projection
      
      
ee              <- ee[order(ee$entity,ee$ddt),]
      
ee$loss_gms     <- ifelse(ee$ddt==idt_q0 + 0.25, ee$loss_gms, 0)
      
ee$extra_items  <- 0
      
      
      
# create the rolling forward 4-qtr sum of NCOs
ee      <- ee %>% arrange(entity,ddt) %>% group_by(entity) %>% mutate(q4_agr = zoo::rollapply(qnetxoff_agr, list(seq(1,4)), sum, fill=0))
ee      <- ee %>% arrange(entity,ddt) %>% group_by(entity) %>% mutate(q4_cc = zoo::rollapply(qnetxoff_cc, list(seq(1,4)), sum, fill=0))
ee      <- ee %>% arrange(entity,ddt) %>% group_by(entity) %>% mutate(q4_ci= zoo::rollapply(qnetxoff_ci, list(seq(1,4)), sum, fill=0))
ee      <- ee %>% arrange(entity,ddt) %>% group_by(entity) %>% mutate(q4_const = zoo::rollapply(qnetxoff_const, list(seq(1,4)), sum, fill=0))
ee      <- ee %>% arrange(entity,ddt) %>% group_by(entity) %>% mutate(q4_firstlien = zoo::rollapply(qnetxoff_firstlien, list(seq(1,4)), sum, fill=0))
ee      <- ee %>% arrange(entity,ddt) %>% group_by(entity) %>% mutate(q4_heloc = zoo::rollapply(qnetxoff_jrlien, list(seq(1,4)), sum, fill=0))
ee      <- ee %>% arrange(entity,ddt) %>% group_by(entity) %>% mutate(q4_jrlien = zoo::rollapply(qnetxoff_heloc, list(seq(1,4)), sum, fill=0))
ee      <- ee %>% arrange(entity,ddt) %>% group_by(entity) %>% mutate(q4_lease = zoo::rollapply(qnetxoff_lease, list(seq(1,4)), sum, fill=0))
ee      <- ee %>% arrange(entity,ddt) %>% group_by(entity) %>% mutate(q4_multi = zoo::rollapply(qnetxoff_multi, list(seq(1,4)), sum, fill=0))
ee      <- ee %>% arrange(entity,ddt) %>% group_by(entity) %>% mutate(q4_nfnr = zoo::rollapply(qnetxoff_nfnr, list(seq(1,4)), sum, fill=0))
ee      <- ee %>% arrange(entity,ddt) %>% group_by(entity) %>% mutate(q4_oth = zoo::rollapply(qnetxoff_oth, list(seq(1,4)), sum, fill=0))
ee      <- ee %>% arrange(entity,ddt) %>% group_by(entity) %>% mutate(q4_othcons = zoo::rollapply(qnetxoff_othcons, list(seq(1,4)), sum, fill=0))
ee      <- ee %>% arrange(entity,ddt) %>% group_by(entity) %>% mutate(q4_othre = zoo::rollapply(qnetxoff_othre, list(seq(1,4)), sum, fill=0))
      
      
nv            <- names(ee)[which(substring(names(ee),1,3)=="q4_")]
ee$q4_total   <- rowSums(ee[,nv],na.rm=TRUE)
      
      
nv                  <- names(ee)[which(substring(names(ee),1,9)=="qnetxoff_")]
ee$qnetxoff_total   <- rowSums(ee[,nv],na.rm=TRUE)
      
      
gg            <- ee[which(ee$ddt==idt_q0),c("entity","llres","q4_total")]
gg$trueup     <- (gg$q4_total - gg$llres)
      
      
ee            <- merge(ee,gg[,c("entity","trueup")],by="entity",all.x=T)
ee            <- ee[order(ee$entity,ee$ddt),]
      
ee$stepN      <- 1- ((ee$ddt - idt_q0)/0.25)/ true_horiz
ee$stepN      <- ifelse(ee$stepN < 0, 0, ee$stepN)
      
      
ee      <- ee %>% arrange(entity,ddt) %>% group_by(entity) %>% mutate(b_lower = q4_total - trueup * stepN)
ee      <- ee %>% arrange(entity,ddt) %>% group_by(entity) %>% mutate(b_upper = 2.5 * b_lower)
      

###
### function definition
###

calc_res <- function(x_res, x_lower, x_upper, x_adj, x_horiz) {
  for (tt in 2:length(x_res)) {
          
    x_res[tt] <- x_res[tt-1] + ifelse(tt <= x_horiz + 1, x_adj[tt], 0)
    x_res[tt] <- ifelse(x_res[tt] >=  x_upper[tt], x_upper[tt], x_res[tt])
    x_res[tt] <- ifelse(x_res[tt] <=  x_lower[tt], x_lower[tt], x_res[tt])
          
  }
        
  x_res
        
}

###
###
###

ee              <- ee %>% arrange(entity,ddt) %>% group_by(entity) %>% mutate(llresX = calc_res(llres, b_lower, b_upper, trueup / true_horiz, true_horiz))
      
ee              <- ee %>% group_by(entity) %>% mutate(qprov = qnetxoff_total + llresX - dplyr::lag(llresX,n=1,order_by=ddt))
      
      

      
##########################################################################
##########################################################################
##########################################################################
      
# Part 9 -- assemble the necessary pieces of capital and produce the final projections
      
      
      
ee$cl_dividend  <- ifelse(ee$cl_dividend < 0, 0, ee$cl_dividend)
ee$cl_dividend[which(ee$ddt==idt_q0)] <- 0
      
      
      
ee$ppnr         <- ee$qint_inc_net + ee$qnonint_inc + ee$qtradrev_inc - ee$qnonint_exp_comp - ee$qnonint_exp_fass - ee$qnonint_exp_allother
      
ee$pretax       <- ee$ppnr + ee$qafssecur_inc + ee$qhtmsecur_inc -ee$qprov - ee$loss_ops - ee$loss_gms
      
      
ee$taxes        <- ee$pretax * tax_rate
      
ee$netinc       <- ee$pretax - ee$taxes + ee$qminor_int + ee$extra_items
      
ee$netinc[which(is.na(ee$netinc))] <- 0
      
ee              <- ee %>% arrange(entity,ddt) %>% group_by(entity) %>% mutate(netincX = cumsum(netinc))
ee              <- ee %>% arrange(entity,ddt) %>% group_by(entity) %>% mutate(dividendX = cumsum(cl_dividend))
      
      
ee$equityX      <- ee$equity + ee$netincX - ee$dividendX 
      
ee$taxesX       <- ifelse(ee$taxes > 0, 0, ee$taxes) 
      
ee$taxesX[which(ee$ddt==idt_q0)] <- 0
      
ee              <- ee %>% arrange(entity,ddt) %>% group_by(entity) %>% mutate(taxesX2 = cumsum(taxesX))
      
      
ee$cet1X        <- ee$cet1 + ee$netincX - ee$dividendX + ee$taxesX2
      
ee$req.cet1     <- ee$req * ee$rwa
      
      
gg              <- ee[which(ee$ddt <= idt_q0 + 9*0.25),c("ddt","entity","cet1X","rwa","avg_assets","req.cet1")]
      
gg$cet1         <- gg$cet1X
gg$cet1X        <- NULL

gg$breach       <- ifelse(gg$cet1 < gg$req.cet1, 1, 0)
gg$cet1         <- gg$cet1 / 1e6
gg$rwa          <- gg$rwa / 1e6
gg$avg_assets   <- gg$avg_assets / 1e6

gg_q0           <- gg[which(gg$ddt==idt_q0),c("entity","cet1")]
names(gg_q0)    <- c("entity","cet1_q0")
gg              <- merge(gg,gg_q0,by="entity",all.x=T,all.y=T)


gg              <- merge(banks,gg,by="entity",all.x=T,all.y=T)
gg$cap_ratio    <- gg$cet1 / gg$rwa * 100
gg$lev_ratio    <- gg$cet1 / gg$avg_assets * 100
gg$cet1_drop    <- gg$cet1 - gg$cet1_q0
gg              <- gg[,c("entity","tkr","ddt","cet1","rwa","avg_assets","cap_ratio",
                         "lev_ratio","cet1_drop","breach")]
                         

      
kk              <- gg %>% group_by(entity) %>% summarize(breach=sum(breach))
kk$breach       <- ifelse(kk$breach > 0, 1, 0)
kk              <- merge(banks,kk,by="entity")      
      

hh              <- gg %>% group_by(ddt) %>% summarize(
                    cet1=sum(cet1), 
                    rwa=sum(rwa), 
                    avg_assets=sum(avg_assets)
                    )


hh              <- data.frame(hh)
hh$cap_ratio    <- hh$cet1 / hh$rwa * 100
hh$lev_ratio    <- hh$cet1 / hh$avg_assets * 100
hh$cet1_drop    <- hh$cet1 - hh$cet1[which(hh$ddt==idt_q0)]
 

     
##########################################################################
##########################################################################
##########################################################################

# Part 10 -- adjust the date variable and store the final results in the output folder

rdt         <- data.frame(date=as.numeric(levels(pp$date)),ddt=0,stringsAsFactors=F)
rdt         <- rdt[which(rdt$date >= date_q0),]
rdt$x       <- c(1:nrow(rdt))
rdt$ddt     <- (rdt$x-1) * 0.25 + idt_q0



hh          <- merge(rdt[,c("date","ddt")],hh,by="ddt",all.x=F,all.y=F)
hh$ddt      <- NULL

gg          <- merge(rdt[,c("date","ddt")],gg,by="ddt",all.x=F,all.y=F)
gg$ddt      <- NULL
gg          <- gg[order(gg$entity,gg$date),]


write.csv(coef_matrix, paste(wksp,"output\\estimated_model_coefficients.csv",sep=""),row.names=F)
write.csv(hh, paste(wksp,"output\\aggregate_projections.csv",sep=""),row.names=F)
write.csv(gg, paste(wksp,"output\\firm_projections.csv",sep=""),row.names=F)
write.csv(kk, paste(wksp,"output\\firm_breaches.csv",sep=""),row.names=F)




