Assignment1AB.Rmd

---
title : "Assignment1A"
author: "Anisha Vijayan (UIN: 662618335)"
        "Preethi Srinivasan (UIN: 663981973)"
date  : "2/17/2022"
output: html_document
---

```{r setup}
library(tidyverse)
library(lubridate)
library(dplyr)
library(magrittr)
library(pROC)
library(rpart)
library(caret)
library(ROCR)
library(C50)
library(ranger)
library(broom)
library(rpart.plot)
lcdf <- read_csv('lcDataSample.csv')
options(digits = 4)  # to set the total digits to display as 4  (decimals will be displayed upto 2 decimal places)

```



```{r Question 3a }
str(lcdf)

#for extracting numeric columns
data_num <- select_if(lcdf, is.numeric)
num_col <- colnames(data_num)


#for converting char cols to date format
lcdf$issue_d # already in date format


lcdf$last_credit_pull_d 
lcdf$last_credit_pull_d  <- paste(lcdf$last_credit_pull_d , "-01", sep = "")
lcdf$last_credit_pull_d  <- parse_date_time(lcdf$last_credit_pull_d , "myd")

lcdf$earliest_cr_line    
lcdf$earliest_cr_line <- paste(lcdf$earliest_cr_line, "-01", sep = "")
lcdf$earliest_cr_line <- parse_date_time(lcdf$earliest_cr_line, "myd")

lcdf$next_pymnt_d
lcdf$next_pymnt_d <- paste(lcdf$next_pymnt_d, "-01", sep = "")
lcdf$next_pymnt_d <- parse_date_time(lcdf$next_pymnt_d, "myd")



#converting character columns to Category
lcdf <- lcdf %>% mutate_if(is.character, as.factor)

data_category <- select_if(lcdf, is.factor)
fac_col <- colnames(data_category)



#for logical columns
data_log <- select_if(lcdf, is.logical)
log_col <- colnames(data_log)


```



```{r Question 3c}

#To find the count of values in the loan_status column
lcdf %>% group_by(loan_status) %>% summarise(count=n())
#From the results, majority of the customers come under Fully Paid category and 15,377 
#come under Charged Off category

#Graph showing count of people with different loans
ggplot(lcdf, aes(loan_status))+geom_bar(stat="count",color="blue",size=0.25)


ggplot(lcdf, aes(loan_amnt)) + geom_histogram(aes(fill=purpose))+geom_bar(stat="count")

ggplot(lcdf, aes(x=purpose,fill=grade))+geom_bar(stat="count",color="blue",size=0.25)+ theme(axis.text.x=element_text(size=15, angle=90,hjust=0.95,vjust=0.2)) #The graph shows number of people taking loans for different purpose in different grades

# Annual income vs Loan amount
ggplot(lcdf, aes(x=loan_status,y=annual_inc))+geom_point(size=0.25)+ theme(axis.text.x=element_text(size=15, angle=90,hjust=0.95,vjust=0.2)) 

#Total funded amount till date
amount_date <-lcdf%>%group_by(last_pymnt_d) %>% summarise(total=sum(funded_amnt)) 
# Highest funded was in the year 2012 October

#Grade wise Interest Rates
ggplot(lcdf, aes(x=int_rate,y=grade,fill=grade))+geom_point(color="blue",size=0.25)
#  Grade F G has highest Interest Rate

boxplot(loan_amnt~annual_inc, data=lcdf)

```


```{r Question 3d i}
#Total number of fully-paid and charged-off loans
lcdf %>% group_by(loan_status) %>% tally() 
#Fully Paid - 94567 , Charged-off - 15377	

#Removing loan_status other than "Fully Paid' and "Charged Off"
lcdf <- lcdf %>% filter(loan_status == "Fully Paid" | loan_status == "Charged Off")

#summary of data as per grade and subgrade
lcdf %>% group_by(grade,sub_grade) %>% summarise(nLoans=n(), defaults=sum(loan_status=="Charged Off"), defaultRate=defaults/nLoans, avgInterest= mean(int_rate), stdInterest=sd(int_rate), avgLoanAMt=mean(loan_amnt), avgPmnt=mean(total_pymnt))
#Shows Average Interest rate, Standard Interest Rate , Average Loan amount, Average Payment grouped by Grade and SubGrade

```


```{r Question 3d ii}

#How does number of loans, loan amount, interest rate vary by grade
lcdf %>% group_by(grade) %>% tally() #Grades A,B,C,D have more number of loans

#Sum of loan amount for different Grades
lcdf %>% group_by(grade) %>% summarise(sum(loan_amnt))  
# Loan amount for Grade B is the highest


#Average interest rate for Every Grade and Subgrade
lcdf %>% group_by(grade,sub_grade) %>% summarise(mean(int_rate)) 
#From the results, Grade B and Sub Grade B5  has the highest Average interest rate

#summary of data according to grade and subgrade
lcdf %>% group_by(grade,sub_grade) %>% summarise(nLoans=n(), avgInterest= mean(int_rate), stdInterest=sd(int_rate), minInterest = min(int_rate), maxInterest = max(int_rate))
#From the results, Average interest ranges from 5.6% to 12.36 %, Min interest range from 5.3% to 6 % and Max interest Lies between 6.03 % to 14.09%


```


```{r Question 3d iii}

#Term of the loan is the duration between the last-payment-date and the loan issue-date
#Checking the format of these two columns with date values
head(lcdf[, c("last_pymnt_d", "issue_d")])

#issue_d is a date variable (of type date), while last_pymnt_d is of type character (like "Dec-2018", having month-year but no date). 
#We need to change the character type to date:
lcdf$last_pymnt_d <- paste(lcdf$last_pymnt_d, "-01", sep = "")
lcdf$last_pymnt_d <- parse_date_time(lcdf$last_pymnt_d,  "myd")

#to reconfirm the format
head(lcdf[, c("last_pymnt_d", "issue_d")])

#To get the duration in years
x<- as.duration(lcdf$issue_d  %--% lcdf$last_pymnt_d)/dyears(1)
head(x)

#to calculate actual term
lcdf$actualTerm <- ifelse(lcdf$loan_status == "Fully Paid", as.duration(lcdf$issue_d %--% lcdf$last_pymnt_d)/dyears(1), 3)

# plotting grade vs actual term
ggplot(lcdf, aes(x = grade)) + geom_boxplot(aes(fill = actualTerm))

```


```{r Question 3d iv}

#Average recoveries for each loan status
lcdf %>% group_by(loan_status) %>%summarise(avgRec=mean(recoveries))
   #shows that recoveries are there only for the Charged-Off loans


#There are different variables for recoveries, to find the total amount of recoveries 
lcdf %>% group_by(loan_status) %>%summarise(avgRec=mean(recoveries), avgPmnt=mean(total_pymnt), mean(total_rec_prncp), mean(total_rec_int), mean(total_rec_late_fee))
   #so we find that 'recoveries' has the total of recoveries on principal, on interest, and late-fees

```


```{r Question 3d v }

#calculate performance (returns) from a loan
lcdf %>% select(loan_status, int_rate, funded_amnt, total_pymnt) %>% head()

#calculate the annualized percentage return
lcdf$annRet <- ((lcdf$total_pymnt - lcdf$funded_amnt)/lcdf$funded_amnt)*(12/36)*100

lcdf %>% select(loan_status, term,int_rate, funded_amnt, total_pymnt, annRet) %>% head()

#summarize by grade the loan details
lcdf %>% group_by(grade) %>% summarise(nLoans = n(), defaults = sum(loan_status == "Charged Off"), avgInterest = mean(int_rate), stdInterest = sd(int_rate), avgLoanAMt = mean(loan_amnt), avgPmnt = mean(total_pymnt), avgRet = mean(annRet), stdRet = sd(annRet), minRet = min(annRet), maxRet = max(annRet))

# filter out loans where annual return is less than 0
lcdf %>% select(loan_status, int_rate, funded_amnt, total_pymnt, annRet) %>% filter(annRet < 0) %>% head()


# to check if all the above returned loans are from Charged off status
lcdf %>% select(loan_status, int_rate, funded_amnt, total_pymnt, annRet) %>% filter(annRet < 0) %>% count(loan_status)
lcdf %>% select(loan_status, int_rate, funded_amnt, total_pymnt, annRet) %>% filter(annRet < 0 & loan_status != "Charged Off") %>% count()  # returns 0 as there is no lines other than charged off

#summarize by grade for Fully Paid back loans
lcdf %>% filter( loan_status == "Fully Paid") %>% group_by(grade) %>% summarise(nLoans=n(), avgInterest= mean(int_rate), 
avgLoanAmt=mean(loan_amnt), avgPmnt=mean(total_pymnt), avgRet=mean(annRet), minRet=min(annRet), maxRet=max(annRet))


#summarize by grade for Charged odd loans
lcdf %>% filter( loan_status == "Charged Off") %>% group_by(grade) %>% summarise(nLoans=n(), avgInterest= mean(int_rate), 
avgLoanAmt=mean(loan_amnt), avgPmnt=mean(total_pymnt), avgRet = mean(annRet), minRet=min(annRet), maxRet=max(annRet))

#In case the loans were paid back early then actual loan term = issue date - last payment date
# we have already computed this as actualTerm in Q3 iii before


# For charged off loans, they are not paid back by 3 yr term so we can set the actual term as 3 for those
#lcdf$actualTerm <- ifelse(lcdf$loan_status == "Fully Paid", as.duration(lcdf$issue_d %--% lcdf$last_pymnt_d)/dyears(1), 3)

#Then, considering this actual term, the actual annual return is
lcdf$actualReturn <- ifelse( lcdf$actualTerm > 0,
((lcdf$total_pymnt - lcdf$funded_amnt)/lcdf$funded_amnt)*(1/lcdf$actualTerm)*100, 0)


#take a look these variables for the first few rows of data 
lcdf %>% select(loan_status, int_rate, funded_amnt, total_pymnt, annRet, actualTerm, actualReturn) %>% head()


# Loan details for Charged off loans
lcdf %>% select(loan_status, int_rate, funded_amnt, total_pymnt, annRet, actualTerm, actualReturn) %>% filter(loan_status =="Charged Off") %>% head()


# summary of interst rate, actual return (total returns, actual annualized return), actual term
lcdf %>% filter(loan_status %in% c("Charged Off","Fully Paid")) %>% group_by(loan_status) %>% summarise( intRate = mean(int_rate), totRet = mean((total_pymnt - funded_amnt)/funded_amnt), avgActRet = mean(actualReturn),avgActTerm = mean(actualTerm) )

#drill down to see the performance by loan grade and status
lcdf %>% filter(loan_status %in% c("Charged Off","Fully Paid")) %>% group_by(loan_status, grade) %>% summarise( intRate = mean(int_rate), totRet = mean((total_pymnt - funded_amnt)/funded_amnt), avgActRet = mean(actualReturn),avgActTerm = mean(actualTerm) )

#drill down to see the performance by loan grade, subgrade and status for Fully Paid
lcdf %>% filter(loan_status %in% c("Fully Paid")) %>% group_by(loan_status, grade, sub_grade) %>% summarise( intRate = mean(int_rate), totRet = mean((total_pymnt - funded_amnt)/funded_amnt), avgActRet = mean(actualReturn),avgActTerm = mean(actualTerm) )


```




```{r Question 3d vi }

#purpose for the money borrowed
lcdf %>% group_by(purpose) %>% tally()

#Summary of data based on the purpose of the loan arranged in desc order of number of loans
lcdf %>% group_by(purpose) %>% summarise(nLoans = n(), defaults = sum(loan_status == "Charged Off"), defaultRate = defaults/nLoans, avgIntRate = mean(int_rate),  avgLoanAmt = mean(loan_amnt),  avgActRet = mean(actualReturn), avgActTerm = mean(actualTerm)) %>% arrange(desc(nLoans))

#Summary of data based on the purpose of the loan arranged in desc order of avg loan amount
lcdf %>% group_by(purpose) %>% summarise(nLoans = n(), defaults = sum(loan_status == "Charged Off"), defaultRate = defaults/nLoans, avgIntRate = mean(int_rate),  avgLoanAmt = mean(loan_amnt),  avgActRet = mean(actualReturn), avgActTerm = mean(actualTerm)) %>% arrange(desc(avgLoanAmt))


#some of category levels have very few examples
# so we have recoded such categories with very few cases to "other"
lcdf$purpose <- fct_recode(lcdf$purpose, other = "wedding", other = "renewable_energy")


#Plot of loan amount by purpose
boxplot(lcdf$loan_amnt ~ lcdf$purpose, main = "Loan Amount vs Purpose",xlab = "Purpose", ylab = "Loan Amount")

# Defaults vary by purpose
lcdf %>% group_by(purpose) %>% summarise(nLoans = n(),defaults = sum(loan_status == "Charged Off"), defaultRate = defaults/nLoans) %>% arrange(desc(defaultRate))


# loan-grade vary by purpose
table(lcdf$purpose, lcdf$grade) # mostly in B and C grades for each purpose
lcdf %>% select(grade, purpose) %>% table() %>% barplot(legend = rownames(.) , beside = T)


```


```{r Question 3d vii }

#convert emp_length to factor -- with factor levels ordered in a meaningful way
lcdf$emp_length <- factor(lcdf$emp_length, levels = c("n/a", "< 1 year","1 year","2 years", "3 years" ,  "4 years",   "5 years",   "6 years",   "7 years" ,  "8 years", "9 years", "10+ years" ))


#summary by emp_length
lcdf %>% group_by(emp_length) %>% summarise(nLoans=n(), avgAnnlInc = mean(annual_inc), defaults=sum(loan_status=="Charged Off"), defaultRate=defaults/nLoans, avgIntRate=mean(int_rate),  avgLoanAmt=mean(loan_amnt),  avgActRet = mean(actualReturn), avgActTerm=mean(actualTerm)) 


#loan-grade assigned by LC vary vs emp_length
table(lcdf$grade, lcdf$emp_length)
lcdf %>% select(grade, emp_length) %>% table() %>% barplot(main = "Emp Length vs Grade", xlab = "Emp Length" ,beside = T, xpd = F)

#Loan Status vs the emp length
df <- lcdf %>% filter(loan_status %in% c("Fully Paid", "Charged Off"))
df %>% group_by(loan_status, emp_length) %>% tally()

#loan purpose vs emp_length
table(lcdf$purpose, lcdf$emp_length)

```



```{r Question 3d viii}

#Derived attributes
#Proportion of satisfactory bankcard accounts 
lcdf$propSatisBCAccts <- ifelse(lcdf$num_bc_tl > 0, lcdf$num_bc_sats/lcdf$num_bc_tl, 0)

#length of borrower's history with LC i.e time between earliest_cr_line and issue_d
lcdf$borrHistory <- as.duration(lcdf$earliest_cr_line %--% lcdf$issue_d  ) / dyears(1)


#Ratio of openAccounts to totalAccounts
lcdf$ratioOpenAccount <- ifelse(lcdf$open_acc > 0, lcdf$open_acc/lcdf$total_acc, 0)

```

 
```{r Question 3d2}
# percentage of loans on loan statuses
lcdf %>% group_by(loan_status) %>% tally()


#Charged Off 13.98%
#Fully Paid 86.01%

```



```{r Question 3e Missing Values}

# to handle missing values
lcdf <- lcdf %>% select_if(function(x){ ! all(is.na(x)) } )
#  all(is.na(x)) will evaluate to True if all the values in x are missing.
#  So we keep those variables x which do NOT have all values missing
 

#Of the columns remaining, names of columns with missing values
names(lcdf)[colSums(is.na(lcdf)) > 0]
# colSums ( is.na( lcdf ) ) returns the total number of True (i.e. NA) values in each column of lcdf
# We then get the names of these columns


#missing value proportions in each column
colMeans(is.na(lcdf))
# or, get only those columns where there are missing values
colMeans(is.na(lcdf))[colMeans(is.na(lcdf))>0]


#Consider open_acc_6m, which has 97% missing
summary(as.factor(lcdf$open_acc_6m))    # shows the counts by different values of the variable
table(lcdf$open_acc_6m)  #gives the same output  -- but it does not show the NAs

# We can replace missing values in a variable with
#      replace_na( variable, "value for missing")     
table( replace_na( lcdf$open_acc_6m, "missing") )   # shows the 'missing' values
table( lcdf$loan_status, replace_na( lcdf$open_acc_6m, "missing") ) 
# shows counts by loan_status at different values of the variable

#to get a bar-plot of these
cc<-table( lcdf$loan_status, replace_na( lcdf$open_acc_6m, "missing") )
barplot(cc, col=c("darkblue","red"),legend = rownames(cc))  # here, one bar dominates others
# For a better display, we can get proportion of ChargedOff as cc[1,]/(cc[2,]+cc[1,]).  Then to plot this..
barplot(cc[1,]/(cc[2,]+cc[1,]), legend = rownames(cc), ylab = "prop ChargedOff", main="Prop ChargedOff by open_acc_6m")




#Consider the "mths_since_" variables -- what do they represent (see data dictionary)
# to check if these are  missing values here due to zeros; or due to no known values in the period considered (then the actual value would be larger than the max value)

#  Variable mths_since_last_record has more than 80% values missing
cc<-table( lcdf$loan_status, replace_na( lcdf$mths_since_last_record, "missing") )
cc[1,]/(cc[2,]+cc[1,])


#For mths_since_last_delinq, which has around 50% values missing 
cc<-table( lcdf$loan_status, replace_na( lcdf$mths_since_last_delinq, "missing") )
cc[1,]/(cc[2,]+cc[1,])
   

#For mths_since_recent_inq, which has around 10% values missing
cc<-table( lcdf$loan_status, replace_na( lcdf$mths_since_recent_inq, "missing") )
cc[1,]/(cc[2,]+cc[1,])

#Suppose you decide to remove variables which have more than 60% missing values
nm<-names(lcdf)[colMeans(is.na(lcdf))>0.6]
lcdf <- lcdf %>% select(-all_of(nm))

#Impute missing values for remaining variables which have missing values
# - first get the columns with missing values
colMeans(is.na(lcdf))[colMeans(is.na(lcdf))>0]

#summary of data in these columns
nm<- names(lcdf)[colSums(is.na(lcdf))>0]
summary(lcdf[, nm])


#For bc_open_to_buy, suppose we want to replace the missing values by the median
#  -- we will try this out and put results in a temporary dataset lcx, with the attributes that have missing values
lcx <-lcdf[, c(nm)]
lcx <- lcx %>% replace_na(list(bc_open_to_buy=median(lcx$bc_open_to_buy, na.rm=TRUE)))


#Similarly for the other variables
#After trying this out on the temporary dataframe lcx, if we are sure this is what we want, we can now  replace the missing values on the lcdf dataset

lcdf<- lcdf %>% replace_na(list(mths_since_last_delinq=-500, bc_open_to_buy=median(lcdf$bc_open_to_buy, na.rm=TRUE), mo_sin_old_il_acct=1000, mths_since_recent_bc=1000, mths_since_recent_inq=50, num_tl_120dpd_2m = median(lcdf$num_tl_120dpd_2m, na.rm=TRUE),percent_bc_gt_75 = median(lcdf$percent_bc_gt_75, na.rm=TRUE), bc_util=median(lcdf$bc_util, na.rm=TRUE) ))

#Lets check if this has addressed all missing values
colMeans(is.na(lcdf))[colMeans(is.na(lcdf))>0]

#Variables with missing values
nm<-names(lcdf)[colMeans(is.na(lcdf))>0]
glimpse(lcdf %>% select(nm))
   #we notice that these are all numeric variables  -- replace by the median values

#To replace the few missing values in a column by the column median values
# -- Try this before making changes in lcdf 
lcx <- lcdf  #copy to lcx
lcx <- lcx %>% mutate_if(is.numeric,  ~ifelse(is.na(.x), median(.x, na.rm = TRUE), .x))
     # if any column has  missing values, replace with median value in that column

#If this works, do the same in lcdf
lcdf<- lcdf %>% mutate_if(is.numeric,  ~ifelse(is.na(.x), median(.x, na.rm = TRUE), .x))


dim(lcdf)  #how many variables left 


```


```{r Question 4 Leakage}

#Dropping some variable/columns which are not useful or which we will not use in developing predictive models
#Also dropping variables those which will cause 'leakage' 

lcdf <- lcdf %>% mutate_if(is.character, as.factor)


#Identify the variables to remove
varsToRemove = c('funded_amnt_inv', 'term', 'emp_title', 'pymnt_plan', 'earliest_cr_line', 'title', 'zip_code', 'addr_state', 'out_prncp', 'out_prncp_inv', 'total_pymnt_inv', 'total_rec_prncp', 'total_rec_int', 'total_rec_late_fee', 'recoveries', 'collection_recovery_fee', 'last_credit_pull_d', 'policy_code', 'disbursement_method', 'debt_settlement_flag',  'application_type')
 

#Drop them from the lcdf data-frame
lcdf <- lcdf %>% select(-all_of(varsToRemove))  


#Drop all the variables with names starting with "hardship" 
lcdf <- lcdf %>% select(-starts_with("hardship"))

#similarly, all variable starting with "settlement"
lcdf <- lcdf %>% select(-starts_with("settlement"))

#some additional variables to drop which we will not use in following analyses
varsToRemove2 <- c("last_pymnt_d", "last_pymnt_amnt", "issue_d")
lcdf <- lcdf %>% select(-all_of(varsToRemove2))

```


```{r Question 5 Univariate Anlayses}

# Ordering some of the factor variables before plotting ROC
lcdf$emp_length <- ordered(lcdf$emp_length, levels = c("n/a", "< 1 year", "1 year", "2 years", "3 years", "4 years", "5 years", "6 years", "7 years", "8 years", "9 years", "10+ years"))
lcdf$grade <- as.ordered(lcdf$grade)
lcdf$sub_grade <- as.ordered(lcdf$sub_grade)
lcdf$loan_amnt <- ordered(lcdf$loan_amnt)
lcdf$loan_status <- factor(lcdf$loan_status, levels = c("Fully Paid","Charged Off"))
lcdf$loan_amnt <- as.numeric(lcdf$loan_amnt)


# returns the auc value for loan_amt as the single predictor for loan_status
roc_obj <- roc(lcdf$loan_status,lcdf$loan_amnt )
auc(roc_obj)  #Area under the curve: 0.515

# For a factor variable, we can consider the factor levels as numbers:
roc_obj1 <- roc(lcdf$loan_status,as.numeric(lcdf$emp_length))
auc(roc_obj1)  #Area under the curve: 0.529


#For both numeric and factor variables:
aucAll <- sapply(lcdf %>% mutate_if(is.factor, as.numeric) %>% select_if(is.numeric), auc, response = lcdf$loan_status) 


#To determine which variables have auc > 0.5
aucAll[aucAll > 0.5]


tidy(aucAll[aucAll > 0.5]) %>% view()

#transfer the data to a csv file
#write.csv(tidy(aucAll[aucAll > 0.5]), file = "text.csv", row.names = FALSE)

#chi square test for purpose
chisq.test(lcdf$loan_status,lcdf$purpose)

```


```{r Question 6 DT}

#Splitting Training Data set and Testing Data Set
lcdf$loan_status <- factor(lcdf$loan_status, levels=c("Fully Paid", "Charged Off"))
class(lcdf$loan_status)

# For training set 50%
TRNPROP = 0.5  #proportion of examples in the training sample
nr <- nrow(lcdf)
trnIndex <- sample(1:nr, size = round(TRNPROP * nr), replace=FALSE)

lcdfTrn <- lcdf[trnIndex, ]
lcdfTst <- lcdf[-trnIndex, ]


#DT models using rpart

#Take a look at teh data
glimpse(lcdf)


#variables like actualTerm, actualReturn, 'annRet', 'total_pymnt'must be removed to avoid leakage
#These will be useful in performance assessment, but should not be used in building the model.
varsOmit <- c('actualTerm', 'actualReturn', 'annRet', 'total_pymnt', 'emp_length', 'delinq_2yrs',
'inq_last_6mths',
'mths_since_last_delinq',
'pub_rec',
'revol_bal',
'revol_util',
'total_acc',
'acc_now_delinq',
'tot_coll_amt',
'tot_cur_bal',
'total_rev_hi_lim',
'acc_open_past_24mths',
'avg_cur_bal',
'bc_open_to_buy',
'bc_util',
'chargeoff_within_12_mths',
'delinq_amnt',
'mo_sin_old_il_acct',
'mo_sin_old_rev_tl_op',
'mo_sin_rcnt_rev_tl_op',
'mo_sin_rcnt_tl',
'mort_acc',
'mths_since_recent_bc',
'mths_since_recent_inq',
'num_tl_30dpd',
'pct_tl_nvr_dlq',
'percent_bc_gt_75',
'pub_rec_bankruptcies',
'num_bc_sats',
'tax_liens')



#Check of the target, loan_status, is a factor variable -- if not, convert to  a factor variable
lcdfTrn$loan_status <- factor(lcdfTrn$loan_status, levels=c("Fully Paid", "Charged Off"))
lcdfTst$loan_status <- factor(lcdfTst$loan_status, levels=c("Fully Paid", "Charged Off"))


dim(lcdf)
dim(lcdfTrn)
dim(lcdfTst)



#MODEL 1

#Build a DT with training data, removing all the variables that are not required with min split = 30
lcDT1 <- rpart(loan_status ~., data=lcdfTrn %>% select(-all_of(varsOmit)), method="class", parms = list(split = "information"), control = rpart.control(cp=0.0001, minsplit = 30))



printcp(lcDT1)

#variable importance
lcDT1$variable.importance


#Pruned tree with cp = 0.00015
lcDT1p <- prune.rpart(lcDT1, cp = 0.00015) 
printcp(lcDT1p)
lcDT1p$variable.importance

#Evaluate performance

# Performance of unpruned decision tree on Training data
predTrn = predict(lcDT1,lcdfTrn, type='class')
table(pred = predTrn, true = lcdfTrn$loan_status)
mean(predTrn == lcdfTrn$loan_status)

# Performance of unpruned decision tree on Test data
predTst = predict(lcDT1,lcdfTst, type='class')
table(pred = predTst, true=lcdfTst$loan_status)
mean(predTst == lcdfTst$loan_status)



# Performance of pruned decision tree on Training data
predTrnP = predict(lcDT1p,lcdfTrn, type='class')
table(pred = predTrnP, true = lcdfTrn$loan_status)
mean(predTrnP == lcdfTrn$loan_status)

# Performance of pruned decision tree on Test data
predTstP = predict(lcDT1p,lcdfTst, type='class')
table(pred = predTstP, true=lcdfTst$loan_status)
mean(predTstP == lcdfTst$loan_status)


#confusionMatrix 

confusionMatrix(predTrn, lcdfTrn$loan_status)
confusionMatrix(predTrnP, lcdfTrn$loan_status)
confusionMatrix(predTst, lcdfTst$loan_status)
confusionMatrix(predTstP, lcdfTst$loan_status)


#ROC plot for MODEL 1

score=predict(lcDT1,lcdfTst, type="prob")[,"Charged Off"]
pred=prediction(score, lcdfTst$loan_status, label.ordering = c("Fully Paid", "Charged Off"))
    #label.ordering here specifies the 'negative', 'positive' class labels   

#ROC curve
aucPerf <-performance(pred, "tpr", "fpr")
plot(aucPerf) 
abline(a=0, b= 1)

#AUC value
aucPerf=performance(pred, "auc")
aucPerf@y.values


#Lift curve
liftPerf <-performance(pred, "lift", "rpp")
plot(liftPerf)




#MODEL 2
#For Training Set = 50% using minsplit of 70
#Build a DT with training data, removing all the variables that are not required with min split = 70
lcDT2 <- rpart(loan_status ~., data=lcdfTrn %>% select(-all_of(varsOmit)), method="class", parms = list(split = "information"), control = rpart.control(cp=0.0001, minsplit = 70))

printcp(lcDT2)

#variable importance
lcDT2$variable.importance


#Pruned tree with cp = 0.00015
lcDT2p <- prune.rpart(lcDT2, cp = 0.00015) 
printcp(lcDT2p)


#Evaluate performance

# Performance of unpruned decision tree on Training data
pred2Trn = predict(lcDT2,lcdfTrn, type='class')
table(pred = pred2Trn, true = lcdfTrn$loan_status)
mean(pred2Trn == lcdfTrn$loan_status)

# Performance of unpruned decision tree on Test data
pred2Tst = predict(lcDT2,lcdfTst, type='class')
table(pred = pred2Tst, true=lcdfTst$loan_status)
mean(pred2Tst == lcdfTst$loan_status)



# Performance of pruned decision tree on Training data
pred2TrnP = predict(lcDT2p,lcdfTrn, type='class')
table(pred = pred2TrnP, true = lcdfTrn$loan_status)
mean(pred2TrnP == lcdfTrn$loan_status)

# Performance of pruned decision tree on Test data
pred2TstP = predict(lcDT2p,lcdfTst, type='class')
table(pred = pred2TstP, true=lcdfTst$loan_status)
mean(pred2TstP == lcdfTst$loan_status)


#confusionMatrix 

confusionMatrix(pred2Trn, lcdfTrn$loan_status)
confusionMatrix(pred2TrnP, lcdfTrn$loan_status)
confusionMatrix(pred2Tst, lcdfTst$loan_status)
confusionMatrix(pred2TstP, lcdfTst$loan_status)

#ROC plot for MODEL 2

score=predict(lcDT2,lcdfTst, type="prob")[,"Charged Off"]
pred=prediction(score, lcdfTst$loan_status, label.ordering = c("Fully Paid", "Charged Off"))
    #label.ordering here specifies the 'negative', 'positive' class labels   

#ROC curve
aucPerf <-performance(pred, "tpr", "fpr")
plot(aucPerf) 
abline(a=0, b= 1)

#AUC value
aucPerf=performance(pred, "auc")
aucPerf@y.values


#Lift curve
liftPerf <-performance(pred, "lift", "rpp")
plot(liftPerf)




#Building a decision tree using C.50

#C50 : Decision Tree
#build a tree model
c5_DT1 <- c5_DT1 <- C5.0(loan_status ~., data=lcdfTrn %>%  select(-all_of(varsOmit)),  control=C5.0Control(minCases=30,earlyStopping =FALSE))

#model details
summary(c5_DT1)



#You may find that the tree has only one root node due to the class imbalance in the data
lcdfTrn %>% group_by(loan_status) %>% tally()
   #show about 6 times more 'Fully Paid' than 'Charged Off' loans


#To consider a more balanced data for building the tree, C5.0 has a 'weights' parameter - this can specify a vector of weights for each example
#Suppose we want to weight the 'Charged Off' examples as 6, and 'Fully Paid' examples as 1
caseWeights <- ifelse(lcdfTrn$loan_status=="Charged Off", 6, 1)


#Then use these caseWeights in the C5.0 function
c5_DT1 <- C5.0(loan_status ~., data=lcdfTrn %>%  select(-all_of(varsOmit)), weights = caseWeights, control=C5.0Control(minCases=30,earlyStopping =FALSE))

summary(c5_DT1)

predc5Trn <- predict(c5_DT1, lcdfTrn, type='prob')
head(predc5Trn)
   #this show two columns,  with scores ('prob') for each class label

CTHRESH=0.5
table(pred = predc5Trn[,'Fully Paid' ] > CTHRESH, true=lcdfTrn$loan_status)

predc5Tst <- predict(c5_DT1, lcdfTst, type='prob')
table(pred = predc5Tst[,'Fully Paid' ] > CTHRESH, true=lcdfTst$loan_status)


#Rules
c5_rules1 <- C5.0(loan_status ~., data=lcdfTrn %>%  select(-all_of(varsOmit)), weights = caseWeights, rules=TRUE, control=C5.0Control(minCases=30))

summary(c5_rules1)


predTrnc5 <- predict(c5_DT1, lcdfTrn, type='class')
confusionMatrix(predTrnc5, lcdfTrn$loan_status)

predTstc5 <- predict(c5_DT1, lcdfTst, type='class')
confusionMatrix(predTstc5, lcdfTst$loan_status)



#ROC plot for C5.0

score=predict(c5_DT1,lcdfTst, type="prob")[,"Charged Off"]
pred=prediction(score, lcdfTst$loan_status, label.ordering = c("Fully Paid", "Charged Off"))
    #label.ordering here specifies the 'negative', 'positive' class labels   

#ROC curve
aucPerf <-performance(pred, "tpr", "fpr")
plot(aucPerf) 
abline(a=0, b= 1)

#AUC value
aucPerf=performance(pred, "auc")
aucPerf@y.values


#Lift curve
liftPerf <-performance(pred, "lift", "rpp")
plot(liftPerf)


#Model 1 lcDT1p is better in comparison to others

```


```{r Functions for Evaluating Performance}

#ROC curve and AUC value
fnROCPerformance <- function(scores, dat) {  #Note the label-ordering - so, scores should be prob of 'Fully Paid'
    pred=prediction(scores, dat$loan_status, label.ordering = c("Charged Off", "Fully Paid" ))

  #ROC curve
  aucPerf <-performance(pred, "tpr", "fpr")
  plot(aucPerf)
  abline(a=0, b= 1)

  #AUC value
  aucPerf=performance(pred, "auc")
  sprintf("AUC: %f", aucPerf@y.values)
        
}


#decile lift performance
#   the 'score' parameter should give 'prob' of loan_status == 'Charged Off'
fnDecileLiftsPerformance_defaults  <- function( scores, dat) {  #score is for loan_status=='Charged Off'
  totDefRate= sum(dat$loan_status=="Charged Off")/nrow(dat)
  decPerf <- data.frame(scores)
  decPerf <- cbind(decPerf, status=dat$loan_status, grade=dat$grade)
  decPerf <- decPerf %>% mutate(decile = ntile(-scores, 10))
  decPerf<-  decPerf  %>% group_by(decile) %>% summarise ( 
    count=n(), numDefaults=sum(status=="Charged Off"), defaultRate=numDefaults/count,
    totA=sum(grade=="A"),totB=sum(grade=="B" ), totC=sum(grade=="C"), totD=sum(grade=="D"),
    totE=sum(grade=="E"),totF=sum(grade=="F") )
  decPerf$cumDefaults=cumsum(decPerf$numDefaults)                      
  decPerf$cumDefaultRate=decPerf$cumDefaults/cumsum(decPerf$count)                      
  decPerf$cumDefaultLift<- decPerf$cumDefaultRate/(sum(decPerf$numDefaults)/sum(decPerf$count))
  
  print(decPerf)
}


#Returns performance by deciles
fnDecileReturnsPerformance <- function( scores, dat) {
  decRetPerf <- data.frame(scores)
  decRetPerf <- cbind(decRetPerf, status=dat$loan_status, grade=dat$grade, actRet=dat$actualReturn, actTerm = dat$actualTerm)
  decRetPerf <- decRetPerf %>% mutate(decile = ntile(-scores, 10))
  decRetPerf %>% group_by(decile) %>% summarise (
    count=n(), numDefaults=sum(status=="Charged Off"), avgActRet=mean(actRet), minRet=min(actRet), maxRet=max(actRet),
    avgTer=mean(actTerm), totA=sum(grade=="A"), totB=sum(grade=="B" ), totC=sum(grade=="C"), totD=sum(grade=="D"),
    totE=sum(grade=="E"), totF=sum(grade=="F") )
}


  
```


```{r Question 7 }

#Model1
rfModel1 <- ranger(loan_status ~., data=lcdfTrn %>%  select(-all_of(varsOmit)), num.trees = 200, importance='permutation', probability = TRUE)

#variable importance
vimp_rfGp<-importance(rfModel1)
vimp_rfGp %>% view()

#Get the predictions -- look into the returned object
scoreTrn <- predict(rfModel1,lcdfTrn)
head(scoreTrn$predictions)

#classification performance , at specific threshold 
table(pred = scoreTrn$predictions[, "Fully Paid"] > 0.7, actual=lcdfTrn$loan_status)

scoreTst <- predict(rfModel1,lcdfTst)
table(pred = scoreTst$predictions[, "Fully Paid"] > 0.7, actual=lcdfTst$loan_status)

#ROC curve, AUC on Training data
pred=prediction(scoreTrn$predictions[, "Fully Paid"], lcdfTrn$loan_status, label.ordering = c("Charged Off","Fully Paid" ))  #ROC curve
aucPerf <-performance(pred, "tpr", "fpr")
plot(aucPerf , col="red")
abline(a=0, b= 1)


#ROC curve, AUC on Test data
pred=prediction(scoreTst$predictions[, "Fully Paid"], lcdfTst$loan_status, label.ordering = c("Charged Off","Fully Paid" ))  #ROC curve
aucPerf <-performance(pred, "tpr", "fpr")
plot(aucPerf , col="red")
abline(a=0, b= 1)

#AUC value
aucPerf=performance(pred, "auc")
sprintf("AUC: %f", aucPerf@y.values)


#Or call the performance function defined above
fnROCPerformance(predict(rfModel1,lcdfTst)$predictions[,"Fully Paid"], dat=lcdfTst)

#for decile defaults-lift performance
fnDecileLiftsPerformance_defaults( predict(rfModel1,lcdfTrn)$predictions[,"Charged Off"], lcdfTrn  )

#for decile returns performance
fnDecileReturnsPerformance( predict(rfModel1,lcdfTrn)$predictions[,"Fully Paid"], lcdfTrn  )
     



# Model 2
#Different parameters for random forest - for example, if the default model is seen to overfit
rfModel2 <- ranger(loan_status ~., data=lcdfTrn %>%  select(-all_of(varsOmit)),
                   num.trees =500, probability = TRUE, min.node.size = 50, importance='permutation', max.depth = 15)


#variable importance
vimp_rf2<-importance(rfModel2)
vimp_rf2 %>% view()

#Get the predictions -- look into the returned object
scoreTrn2 <- predict(rfModel2,lcdfTrn)
head(scoreTrn2$predictions)

#classification performance , at specific threshold 
#on training data
table(pred = scoreTrn2$predictions[, "Fully Paid"] > 0.7, actual=lcdfTrn$loan_status)

#on test data
scoreTst2 <- predict(rfModel2,lcdfTst)
table(pred = scoreTst2$predictions[, "Fully Paid"] > 0.7, actual=lcdfTst$loan_status)

#ROC curve, AUC on Training data for Model 2 - looks good
pred=prediction(scoreTrn2$predictions[, "Fully Paid"], lcdfTrn$loan_status, label.ordering = c("Charged Off","Fully Paid" ))  #ROC curve
aucPerf <-performance(pred, "tpr", "fpr")
plot(aucPerf , col="green")
abline(a=0, b= 1)

#ROC curve, AUC on Test data for Model 2 - looks good
pred=prediction(scoreTst2$predictions[, "Fully Paid"], lcdfTst$loan_status, label.ordering = c("Charged Off","Fully Paid" ))  #ROC curve
aucPerf <-performance(pred, "tpr", "fpr")
plot(aucPerf , col="green")
abline(a=0, b= 1)


#AUC value
aucPerf=performance(pred, "auc")
sprintf("AUC: %f", aucPerf@y.values)


#####

# ROC curve for model1 and model 2 on test data
pred_rf200=predict(rfModel1,lcdfTst)$predictions
perfROC_rf200=performance(prediction(predict(rfModel1,lcdfTst)$predictions[,2], lcdfTst$loan_status,label.ordering=c("Fully Paid", "Charged Off")), "tpr", "fpr")

pred_rf500=predict(rfModel2,lcdfTst)$predictions
perfROC_rf500=performance(prediction(predict(rfModel2,lcdfTst)$predictions[,2], lcdfTst$loan_status,label.ordering=c("Fully Paid", "Charged Off")), "tpr", "fpr")

plot(perfROC_rf200,  col="red")
plot(perfROC_rf500, add=TRUE, col="green")
abline(a=0, b= 1)
 # model 2(green) is slightly higher than model1(red)


#lift curve for the random forest model 1 and 2 on test data

perfLift_rfTst200=performance(prediction(predict(rfModel1,lcdfTst)$predictions[,2], lcdfTst$loan_status,label.ordering=c("Fully Paid", "Charged Off")), "lift", "rpp")
perfLift_rfTst500=performance(prediction(predict(rfModel2,lcdfTst)$predictions[,2], lcdfTst$loan_status,label.ordering=c("Fully Paid", "Charged Off")), "lift", "rpp")

#Comparing Lift Curves of two  models

plot(perfLift_rfTst200, col="red")
plot(perfLift_rfTst500, add=TRUE, col="green")
 # for model 2(green) the lift curve looks better than for model1(red)


#Confusion Matrix for both models


pred200 = ifelse(pred_rf200[, "Fully Paid"] >= 0.5, "Fully Paid", "Charged Off")
pred200 <- factor(pred200, levels=c("Fully Paid", "Charged Off"))
ct200 = table(pred = pred200, true = lcdfTst$loan_status)

pred500 = ifelse(pred_rf500[, "Fully Paid"] >= 0.5, "Fully Paid", "Charged Off")
pred500 <- factor(pred500, levels=c("Fully Paid", "Charged Off"))
ct500 = table(pred = pred500, true = lcdfTst$loan_status)

ct200 # confusion matrix for model1
ct500 # confusion matrix for model2

#confusion matrix for both tables along with other measures
confusionMatrix(pred200,lcdfTst$loan_status, positive="Charged Off")
confusionMatrix(pred500,lcdfTst$loan_status, positive="Charged Off")

```


```{r Question 8 }

# to check the avg interest rate for fully paid loans 
lcdf %>% filter(loan_status %in% c("Fully Paid")) %>% summarise( intRate = mean(int_rate) )
  # 11.8

# therefore, we can take the profit value on FP loans = (100+3*11.8) = 135.4 so we take 35 as profit value

# to check the avg interest rate for Charged off loans 
lcdf %>% filter(loan_status %in% c("Charged Off")) %>% summarise( intRate = mean(int_rate) )
  # 13.86

# therefore, we can take the profit value on CO loans = (100+3*13.86) = 141.58 so we take 42 as loss value


pred_DT1p = predict(lcDT1p,lcdfTst)
pred_DT1p_1 = ifelse(pred_DT1p[, "Fully Paid"] >= 0.5, "Fully Paid", "Charged Off")

pred_DT1p_1 <- factor(pred_DT1p_1, levels=c("Fully Paid", "Charged Off"))
CMDT <- table(pred=pred_DT1p_1, true=lcdfTst$loan_status)


#Profit/Loss table
PROFITVAL <- 35
LOSSVAL <- -42

PLtable <- matrix(c(PROFITVAL,6,LOSSVAL ,6),ncol=2,byrow=FALSE)
colnames(PLtable) <- c("Act Full Paid","Act Charged Off")
rownames(PLtable) <- c("Pred Fully Paid","Pred Charged Off")
PLtable <- as.table(PLtable)

#Calculate Profit
#Decision tree
CMDT
sum(CMDT*PLtable)

#Random Forest
ct500
sum(ct500*PLtable)


```


```{r Question 8 a}

#decision tree and random forest comparison on Test data
#for decision tree 1
perfROC_dt1Tst=performance(prediction(predict(lcDT1,lcdfTst, type="prob")[,2], lcdfTst$loan_status), "tpr", "fpr")
#for decision tree 1 pruned
perfROC_dt2Tst=performance(prediction(predict(lcDT1p,lcdfTst, type="prob")[,2], lcdfTst$loan_status), "tpr", "fpr")
#for random forest model2 with 500 trees
perfROC_rf500=performance(prediction(predict(rfModel2,lcdfTst)$predictions[,2], lcdfTst$loan_status,label.ordering=c("Fully Paid", "Charged Off")), "tpr", "fpr")


plot(perfROC_dt1Tst, col='red')
plot(perfROC_dt2Tst, col='blue', add=TRUE)
plot(perfROC_rf500, col='green', add=TRUE)
legend('bottomright', c('DecisionTree-30', 'Pruned DecisionTree', 'RandomForest'), lty=1, col=c('red', 'blue', 'green'))  
 # Random forest looks better compared to other DT models

```


```{r Question 8 b}

#Performance with profit.loss
#Incorporating profits & costs
PROFITVAL <- 35
LOSSVAL <- -42

#performance for Decision Tree model
scoreDT <- predict(lcDT1p,lcdfTst, type="prob") [,"Fully Paid"]  

#Note- we want to identify those loans with high prob for being FullyPaid for DT
prPerf <- data.frame(scoreDT)
prPerf <- cbind(prPerf, status=lcdfTst$loan_status)
prPerf <- prPerf[order(-scoreDT) ,]  #sort in desc order of  prob(fully_paid)
prPerf$profit <- ifelse(prPerf$status == 'Fully Paid', PROFITVAL, LOSSVAL)
prPerf$cumProfit <- cumsum(prPerf$profit)

#what is the max profit for DT model
max( prPerf$cumProfit ) 

prPerf$cumProfit[which.max(prPerf$cumProfit)]

#score for that max profit 
prPerf$scoreDT[which.max(prPerf$cumProfit)]



#Performance for Random Forest
scoreRF <- predict(rfModel2,lcdfTst, type="response")$predictions [,"Fully Paid"]

prPerfRF <- data.frame(scoreRF)
prPerfRF <- cbind(prPerfRF, status=lcdfTst$loan_status)
prPerfRF <- prPerfRF[order(-scoreRF) ,]  #sort in desc order of  prob(fully_paid)
prPerfRF$profit <- ifelse(prPerfRF$status == 'Fully Paid', PROFITVAL, LOSSVAL)
prPerfRF$cumProfit <- cumsum(prPerfRF$profit)

#Max cum profit 
prPerfRF$cumProfit[which.max(prPerfRF$cumProfit)]

#score for that max profit 
prPerfRF$scoreRF[which.max(prPerfRF$cumProfit)]


# Decile Lift performance by default 
# Decision Tree
fnDecileLiftsPerformance_defaults( predict(lcDT1p, lcdfTrn)[,"Fully Paid"], lcdfTrn  ) #on training data
fnDecileLiftsPerformance_defaults( predict(lcDT1p, lcdfTst)[,"Fully Paid"], lcdfTst  ) #on test data

# Random Forest
fnDecileLiftsPerformance_defaults( predict(rfModel2, lcdfTrn, type="response")$predictions[,"Fully Paid"], lcdfTrn  ) #on training data
fnDecileLiftsPerformance_defaults( predict(rfModel2, lcdfTst, type="response")$predictions[,"Fully Paid"], lcdfTst  ) #on test data


# Returns performance by decile
# Decision Tree
fnDecileReturnsPerformance( predict(lcDT1p, lcdfTrn)[,"Fully Paid"], lcdfTrn  ) #on training data
fnDecileReturnsPerformance( predict(lcDT1p, lcdfTst)[,"Fully Paid"], lcdfTst  ) #on test data

# Random Forest
fnDecileReturnsPerformance( predict(rfModel2, lcdfTrn, type="response")$predictions[,"Fully Paid"], lcdfTrn  ) #on training data
fnDecileReturnsPerformance( predict(rfModel2, lcdfTst, type="response")$predictions[,"Fully Paid"], lcdfTst  ) #on test data



#to compare against the default approach of investing in CD with 2% int (i.e. $6 profit out of $100 in 3 years)
prPerfRF$cdRet <- 35
prPerfRF$cumCDRet <- cumsum(prPerfRF$cdRet)
plot(prPerfRF$cumProfit,col='green',lwd = 0.5)
lines(prPerf$cumProfit,col='red',lwd = 0.5)
lines(prPerfRF$cumCDRet, col='yellow',lwd = 0.5)
legend('bottomright', c('cum profit RF','cum profit DT', 'profit'), lty=1, col=c('green','red', 'yellow'))



```