Project v5.Rmd

---
title: "Final Project"
author: "Warren Speth, Jeff McKee, Verdon Tsang, Jacob Bennett"
output:
  pdf_document: default
  html_document: default
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

```{r message=FALSE}
#load required packages
library(caret)
library(data.table)
library(kableExtra)
library(glmnet)
library(rpart)
library(rpart.plot)
library(randomForest)
library(xgboost)
library(gbm)
library(class)

log_prob <- function(pred){
  1/(1+exp(-pred))
}
```

## Code Outline

This code progresses through these major sections:
1. Load and prepare data for analysis (including train/test split)
2. Generate models (using train data)
3. Use models to calculate the number of calls that the Bank should make to customers to maximize profits (on the test data)
4. Assess how profitable each model's strategy would have actually been (on the test data)
5. Compare models
6. k Nearest Neighbors

**Important Note**
In several of the following sections, we have commented out the code that we used to generate a model. This is for time-saving purposes when running this file: you can instead directly load the models that were generated by the code. 


# Section 1: Load and prepare data for analysis (including train/test split)
## Load Data

```{r}
#load data
df <- read.csv("bank_additional_full.csv")

#fix dashes in var names
levels(df$job)[2] <- "bluecollar"
levels(df$job)[7] <- "selfemployed"

#Create outcome variable and model matrix categorical variables 
df <- as.data.table(model.matrix(~.-1,df))
df$Class <- ifelse(df$ytrue==1,TRUE,FALSE)

```

## Data cleaning, create test and train sets

```{r}
#create ID variable
df[,ID := .I]

#split data into test and train
set.seed(99)
training_size <- floor(0.80 * nrow(df))
train_ind <- sample(seq_len(nrow(df)), size = training_size)

df_train <- df[train_ind, ]

df_test <- df[-train_ind, ]

#for purposes of our analysis we removed outliers with more than 5 campaign treatments
#set final outcome to false and campaign to 5 if more than 5 campaign treatments 
df_test[df_test$campaign>5,'Class']<-FALSE
df_test$campaign[df_test$campaign>5]<-5
#turn campaign into dummy variables
df_test[,campaign2 := ifelse(campaign == 2, 1,0)]
df_test[,campaign3 := ifelse(campaign == 3, 1,0)]
df_test[,campaign4 := ifelse(campaign == 4, 1,0)]
df_test[,campaign5 := ifelse(campaign == 5, 1,0)]

#expand train data and set outcome to false if not the final campaign in the data
full_df <- df_train[rep(1:.N,campaign)][,campaign_no:=1:.N,by=ID]
full_df[campaign != campaign_no,'Class'] <- FALSE
#turn campaign into dummy variables
full_df[,campaign2 := ifelse(campaign_no == 2, 1,0)]
full_df[,campaign3 := ifelse(campaign_no == 3, 1,0)]
full_df[,campaign4 := ifelse(campaign_no == 4, 1,0)]
full_df[,campaign5 := ifelse(campaign_no == 5, 1,0)]
#create model matrix with only variables used for model building
model_df<-full_df[campaign_no<=5,c(56,1:30,47:49, 59:62)]
x_matrix <- model.matrix(Class~.*(campaign2+campaign3+campaign4+campaign5)-campaign2*(campaign3+campaign4+campaign5)-campaign3*(campaign4+campaign5)-campaign4*campaign5+campaign2+campaign3+campaign4+campaign5-1,data=model_df)

```

# Section 2: Generate models (using train data)
# Models have been saved in same root folder as code for your convenience

## Logit

```{r}
# test_reg   <-glm(Class~(age+jobadmin.+jobbluecollar+jobentrepreneur+jobhousemaid+jobmanagement+jobretired+jobselfemployed+jobservices+jobstudent+jobtechnician+jobunemployed+jobunknown) *( campaign2+campaign3)+campaign4+campaign5+maritalsingle+maritalunknown+educationbasic.6y+educationbasic.9y+educationhigh.school+educationilliterate+educationprofessional.course+educationuniversity.degree+educationunknown+defaultunknown+defaultyes+housingunknown+housingyes+loanunknown+loanyes,family='binomial',data=model_df)
# save(test_reg, file="test_reg.rda")

load(file="test_reg.rda")
```

## Regularized Regression

```{r}
# # 10 fold cross validation to determine lambda. Test alphas of 1,0.5, and 0 (i.e. LASSO, Enet, and Ridge)
# foldid=sample(1:10,size=length(model_df$Class),replace=TRUE)
# cv_ridge <- cv.glmnet(x=x_matrix, y=model_df$Class, foldid=foldid, lambda = exp(seq(from = -10, to = 0, by = .1)), family="binomial",alpha=0)
# cv_enet <- cv.glmnet(x=x_matrix, y=model_df$Class, foldid=foldid, lambda = exp(seq(from = -10, to = 0, by = .1)), family="binomial",alpha=.5)
# cv_lasso <- cv.glmnet(x=x_matrix, y=model_df$Class, foldid=foldid, lambda = exp(seq(from = -10, to = 0, by = .1)), family="binomial",alpha=1)
# 
# save(cv_ridge, file="cv_ridge.rda")
# save(cv_lasso, file="cv_lasso.rda")
# save(cv_enet, file="cv_enet.rda")

load(file="cv_ridge.rda")
load(file="cv_lasso.rda")
load(file="cv_enet.rda")

#Plot binomial deviance of three models to select best performing
plot(log(cv_lasso$lambda),cv_lasso$cvm,pch=19,col="red",xlab="log(Lambda)",ylab=cv_lasso$name)
points(log(cv_enet$lambda),cv_enet$cvm,pch=19,col="grey")
points(log(cv_ridge$lambda),cv_ridge$cvm,pch=19,col="blue")
legend("topleft",legend=c("alpha= 1","alpha= .5","alpha 0"),pch=19,col=c("red","grey","blue"))

#Output non-zero coefficients to understand what drives LASSO performance
predict(cv_lasso,type="nonzero")
nonzero <- c(1,8,11,25,31,32,34)
lasso_coefficients <- predict(cv_lasso,type='coefficients')
nonzero_coeffs <- row.names(lasso_coefficients)[nonzero]
kable(lasso_coefficients[nonzero_coeffs,])
```

## Regression Tree (with pruning)

```{r}
set.seed(99)

#fit model
#alt: data=model_df
#alt: data=df_train
fit1 = rpart(Class~., data=model_df, control=rpart.control(minsplit=5, cp=0.0001, xval=10), method="class")

nbig = length(unique(fit1$where))
cat('Size of big tree: ',nbig,'\n')

cat('RMSE=',sqrt(mean((df_test$Class-predict(fit1,df_test,type="prob")[,2])^2)),'\n')

#Prune model
cptable = printcp(fit1)
bestcp = cptable[ which.min(cptable[,"xerror"]), "CP" ]
fit1B = prune(fit1,cp=bestcp)
nbigB = length(unique(fit1B$where))
cat('Size of pruned tree: ',nbigB,'\n')

# cat('RMSE=',sqrt(mean((df_test$Class-predict(fit1B,df_test,type="prob")[,2])^2)),'\n')
```

## Random Forest

```{r tree all vars, echo=FALSE}
set.seed(99)
model_df2 <- model_df
model_df2$Class <- as.factor(model_df2$Class)

# fit2 <- randomForest(Class~., data=model_df2, type="classification")
# save(fit2, file="random_forest.rda")
load(file="random_forest.rda")

# cat('RMSE is ',sqrt(mean((df_test$Class-predict(fit2, newdata= df_test, type="prob")[,2])^2)),'\n')
```

## Boosting

```{r boosting}
set.seed(99)
# Inputs, Regression Setting
fboost <- gbm(Class~.,           # regression model
           data=model_df,              # data set
           distribution="bernoulli", # logistic regression for 0-1 outcomes
           n.trees=500,             # Total number of trees/iterations
           interaction.depth = 1,   #1 means additive, 2 means 2-way interaction, etc
           shrinkage=0.1           # Shrinkage parameter, weak predictor
           )

#Output variable importance
varImp(fboost, numTrees = 500)

# Returned Values, Returned predictor
# Making prediction
# cat('RMSE is',sqrt(mean((df_test$Class-log_prob(predict(fboost,newdata=df_test,n.trees=500)))^2)),'\n')
```

# Section 3: Use models to calculate the number of calls that the Bank should make to customers to maximize profits (on the test data)
## Assign Profitiability Parameters 

```{r}
#Relevant Outputs: fail_cost (cost of call resulting in failure), Success_profit (net profit of call resulting in success)

#ASSUMPTIONS:

          wage <-12 #Estimate of hourly wage 
          fail_time <-mean(df$duration[df$Class=="FALSE"]) / 3600
          success_time <-mean(df$duration[df$Class=="TRUE"]) / 3600
          Inter_call_time <-1/60
          fail_cost<-wage*(fail_time+Inter_call_time)
          success_cost <-wage*(success_time+Inter_call_time)
          
#Profit from successful customer: According to a paper, CD rates in the Netherlands average ~60% of Euribor rates (https://onlinelibrary.wiley.com/doi/full/10.1111/irfi.12143)
          Market_rate <-.025 
          CD_rate <-Market_rate*.6
          Balance<-10000 #(simply assumed)
          Term <- .25 #(.25 years, assumed based on rate)
          Revenue <- Balance*Term*(Market_rate-CD_rate)  #(roughly a function of not having to borrow this money at market rates, flat curve assumed)
          Success_profit <-Revenue - success_cost
```

## Function for calculating expected profit per call - LOGIT

```{r}
#Generate estimates of cost per call, benefit (revenue/profit) of time deposit product for bank


Compute_Sample_Profitability<-function(obs){
          sample_person <-df_test[obs]
          
#from regression coefficients + person's info: Simple prediction of success likelihood on the nth call "prob(n)" (conditional on failures for all previous calls)

#Loops through prediction turning on campaigns to calculate probabilities of success on calls 1:5
prob <-rep(NA,5)
    sample_person$campaign2<-0
    sample_person$campaign3<-0
    sample_person$campaign4<-0
    sample_person$campaign5<-0
prob[1] <-log_prob(predict(test_reg,sample_person))
    sample_person$campaign2<-1
prob[2] <-log_prob(predict(test_reg,sample_person))
    sample_person$campaign2<-0
    sample_person$campaign3<-1
prob[3] <-log_prob(predict(test_reg,sample_person))
    sample_person$campaign3<-0
    sample_person$campaign4<-1
prob[4] <-log_prob(predict(test_reg,sample_person))
    sample_person$campaign4<-0
    sample_person$campaign5<-1
prob[5] <-log_prob(predict(test_reg,sample_person))
    sample_person$campaign5<-0


#For Calls 1 through N, compute likelihood of success by the nth call: p_sub(n)
        
p_sub<-rep(NA,5)
  p_sub[1]<-prob[1]
  for (i in 2:5) {
    p_sub[i] <-p_sub[i-1]+prob[i]*(1-p_sub[i-1])
  }

  
#For calls 1 through N, compute lift of call n (this represents unconditional probability of success on call n)
lift<-rep(NA,5)
  lift[1]<-p_sub[1]
  for (i in 2:5) {
    lift[i] <-p_sub[i]-p_sub[i-1]
  }
  
#Compute number of calls to target for a given customer (N.B.: As long as probability of success is declining monotonically, as soon as the nth call is not profitable, the (n+1)th call will also not be profitable.) 
  
  #Test: Identify whether the nth call (if made) would be profitable: Profit(n|All calls prior="no")=[Profit of Success]*p_sub(n) -[Cost of Failure]*(1-p_sub(n)) 
  Profit <- rep(NA,5)
  for (i in 1:5){
    Profit[i] <- Success_profit*prob[i]-fail_cost*(1-prob[i])
  }
  #Decision: for the first unprofitable call "A," choose  A-1 as the target number of calls.
  Calls_to_make<-0
  Calls_to_make <-NROW(Profit[Profit>0])
  #based out this decision (target customer with A-1 calls), also compute: 
    #(1) the expected number of calls with that rule(this is not simply A-1 b/c customer may sign up earlier), and
Expected_calls <- 0
if(Calls_to_make<=1) {
  Expected_calls<-Calls_to_make
}  else {
    for (i in 1:(Calls_to_make-1)){
      Expected_calls<-Expected_calls+i*lift[i]
    }
    Expected_calls<-Expected_calls+Calls_to_make*(1-p_sub[Calls_to_make-1])
  }

    #(2) the expected profit from this strategy (weighted average of profitabilities forscenarios given this targeting plan, weighting each scenario by lift(n) and then adding together)
  
  
  Expected_profit <- 0
  if (Calls_to_make>0){
    for (i in 1:Calls_to_make){
      Expected_profit <-Expected_profit+lift[i]*(Success_profit-fail_cost*(i-1))
    }
    Expected_profit <-Expected_profit-fail_cost*Calls_to_make*(1-p_sub[Calls_to_make])
  } 

  Metrics <-c("Euribor","Max Calls","Expected Calls", "Profit")
  values<-c(Market_rate,Calls_to_make,Expected_calls,Expected_profit)
  results <-data.frame(Metrics,values)
#  print(results)
  
  #Return maximum number of calls to make and success probabilities of each call for decision algorithm 
  return_frame <- data.frame(max_calls = Calls_to_make, p1 = prob[1], p2 = prob[2], p3 = prob[3], p4 = prob[4], p5 = prob[5])
  
  return(return_frame)

}

```

## Function for calculating expected profit per call - Regularized Regression
# Identical to function above but adapted for glm.predict and to take in different models/hyperparameters

```{r}
#Generate estimates of cost per call, benefit (revenue/profit) of time deposit product for bank


Compute_Sample_Profitability_LASSO<-function(obs,model){
          sample_person <-df_test[obs]

#Loops through prediction turning on campaigns to calculate probabilities of success on calls 1:5
prob <-rep(NA,5)
    sample_person$campaign2<-0
    sample_person$campaign3<-0
    sample_person$campaign4<-0
    sample_person$campaign5<-0
    lasso_test<-sample_person[,c(56,1:30,47:49, 58:61)]
    test_matrix <- model.matrix(Class~.*(campaign2+campaign3+campaign4+campaign5)-campaign2*(campaign3+campaign4+campaign5)-campaign3*(campaign4+campaign5)-campaign4*campaign5+campaign2+campaign3+campaign4+campaign5-1,data=lasso_test)
prob[1] <-log_prob(predict(model,newx=test_matrix,s=model$lambda.min))
    sample_person$campaign2<-1
    lasso_test<-sample_person[,c(56,1:30,47:49, 58:61)]
    test_matrix <- model.matrix(Class~.*(campaign2+campaign3+campaign4+campaign5)-campaign2*(campaign3+campaign4+campaign5)-campaign3*(campaign4+campaign5)-campaign4*campaign5+campaign2+campaign3+campaign4+campaign5-1,data=lasso_test)
prob[2] <-log_prob(predict(model,newx=test_matrix,s=model$lambda.min))
    sample_person$campaign2<-0
    sample_person$campaign3<-1
    lasso_test<-sample_person[,c(56,1:30,47:49, 58:61)]
    test_matrix <- model.matrix(Class~.*(campaign2+campaign3+campaign4+campaign5)-campaign2*(campaign3+campaign4+campaign5)-campaign3*(campaign4+campaign5)-campaign4*campaign5+campaign2+campaign3+campaign4+campaign5-1,data=lasso_test)
prob[3] <-log_prob(predict(model,newx=test_matrix,s=model$lambda.min))
    sample_person$campaign3<-0
    sample_person$campaign4<-1
    lasso_test<-sample_person[,c(56,1:30,47:49, 58:61)]
    test_matrix <- model.matrix(Class~.*(campaign2+campaign3+campaign4+campaign5)-campaign2*(campaign3+campaign4+campaign5)-campaign3*(campaign4+campaign5)-campaign4*campaign5+campaign2+campaign3+campaign4+campaign5-1,data=lasso_test)
prob[4] <-log_prob(predict(model,newx=test_matrix,s=model$lambda.min))
    sample_person$campaign4<-0
    sample_person$campaign5<-1
    lasso_test<-sample_person[,c(56,1:30,47:49, 58:61)]
    test_matrix <- model.matrix(Class~.*(campaign2+campaign3+campaign4+campaign5)-campaign2*(campaign3+campaign4+campaign5)-campaign3*(campaign4+campaign5)-campaign4*campaign5+campaign2+campaign3+campaign4+campaign5-1,data=lasso_test)
prob[5] <-log_prob(predict(model,newx=test_matrix,s=model$lambda.min))


    #For Calls 1 through N, compute likelihood of success by the nth call: p_sub(n)
        
p_sub<-rep(NA,5)
  p_sub[1]<-prob[1]
  for (i in 2:5) {
    p_sub[i] <-p_sub[i-1]+prob[i]*(1-p_sub[i-1])
  }

  
  #For calls 1 through N, compute lift of call n (this represents unconditional probability of success on call n)
lift<-rep(NA,5)
  lift[1]<-p_sub[1]
  for (i in 2:5) {
    lift[i] <-p_sub[i]-p_sub[i-1]
  }
  
#Compute number of calls to target for a given customer (N.B.: As long as probability of success is declining monotonically, as soon as the nth call is not profitable, the (n+1)th call will also not be profitable.) 
  #Test: Identify whether the nth call (if made) would be profitable: Profit(n|All calls prior="no")=[Profit of Success]*p_sub(n) -[Cost of Failure]*(1-p_sub(n)) 
  Profit <- rep(NA,5)
  for (i in 1:5){
    Profit[i] <- Success_profit*prob[i]-fail_cost*(1-prob[i])
  }
  #Decision: for the first unprofitable call "A," choose  A-1 as the target number of calls.
  Calls_to_make<-0
  Calls_to_make <-NROW(Profit[Profit>0])
  #based out this decision (target customer with A-1 calls), also compute 
    #(1) the expected number of calls with that rule(this is not simply A-1 b/c customer may sign up earlier), and
Expected_calls <- 0
if(Calls_to_make<=1) {
  Expected_calls<-Calls_to_make
}  else {
    for (i in 1:(Calls_to_make-1)){
      Expected_calls<-Expected_calls+i*lift[i]
    }
    Expected_calls<-Expected_calls+Calls_to_make*(1-p_sub[Calls_to_make-1])
  }

    #(2) the expected profit from this strategy (weighted average of profitabilities forscenarios given this targeting plan, weighting each scenario by lift(n) and then adding together)
  
  
  Expected_profit <- 0
  if (Calls_to_make>0){
    for (i in 1:Calls_to_make){
      Expected_profit <-Expected_profit+lift[i]*(Success_profit-fail_cost*(i-1))
    }
    Expected_profit <-Expected_profit-fail_cost*Calls_to_make*(1-p_sub[Calls_to_make])
  } 

  Metrics <-c("Euribor","Max Calls","Expected Calls", "Profit")
  values<-c(Market_rate,Calls_to_make,Expected_calls,Expected_profit)
  results <-data.frame(Metrics,values)
#  print(results)
  
  return_frame <- data.frame(max_calls = Calls_to_make, p1 = prob[1], p2 = prob[2], p3 = prob[3], p4 = prob[4], p5 = prob[5])
  
  return(return_frame)
}

```

## Function for calculating expected profit per call - Regression Tree

```{r}
sample_tree <-df_test
          
# Loops through prediction turning on campaigns to calculate probabilities of success on calls 1:5

sample_tree$campaign2<-0
sample_tree$campaign3<-0
sample_tree$campaign4<-0
sample_tree$campaign5<-0
tree_test<-sample_tree[,c(56,1:30,47:49, 58:61)]
prob1_tree <- predict(fit1,newdata=tree_test,type="prob")[,2]

tree_test$campaign2<-1
prob2_tree <- predict(fit1,newdata=tree_test,type="prob")[,2]

tree_test$campaign2<-0
tree_test$campaign3<-1
prob3_tree <- predict(fit1,newdata=tree_test,type="prob")[,2]

tree_test$campaign3<-0
tree_test$campaign4<-1
prob4_tree <- predict(fit1,newdata=tree_test,type="prob")[,2]

tree_test$campaign4<-0
tree_test$campaign5<-1
prob5_tree <- predict(fit1,newdata=tree_test,type="prob")[,2]

prob_tree<-cbind(prob1_tree,prob2_tree,prob3_tree,prob4_tree,prob5_tree)

# For Calls 1 through N, compute likelihood of success by the nth call: p_sub(n)
        
p_sub1_tree <- prob_tree[,1]
p_sub2_tree <- p_sub1_tree+prob_tree[,2]*(1-p_sub1_tree)
p_sub3_tree <- p_sub2_tree+prob_tree[,3]*(1-p_sub2_tree)
p_sub4_tree <- p_sub3_tree+prob_tree[,4]*(1-p_sub3_tree)
p_sub5_tree <- p_sub4_tree+prob_tree[,5]*(1-p_sub4_tree)

  
# For calls 1 through N, compute lift of call n: lift(n)= p_sub(n) - p_sub(n-1) (this represents unconditional probability of success on call n)

lift1_tree <- p_sub1_tree
lift2_tree <- p_sub2_tree-p_sub1_tree
lift3_tree <- p_sub3_tree-p_sub2_tree
lift4_tree <- p_sub4_tree-p_sub3_tree
lift5_tree <- p_sub5_tree-p_sub4_tree

#Compute number of calls to target for a given customer (N.B.: As long as probability of success is declining monotonically, as soon as the nth call is not profitable, the (n+1)th call will also not be profitable.) 
  #Test: Identify whether the nth call (if made) would be profitable. THis is simple:
   #Profit(n|All calls prior="no")=[Profit of Success]*p_sub(n) -[Cost of Failure]*(1-p_sub(n)) 

profit1_tree <- Success_profit*prob1_tree - fail_cost*(1-prob1_tree)
profit2_tree <- Success_profit*prob2_tree - fail_cost*(1-prob2_tree)
profit3_tree <- Success_profit*prob3_tree - fail_cost*(1-prob3_tree)
profit4_tree <- Success_profit*prob4_tree - fail_cost*(1-prob4_tree)
profit5_tree <- Success_profit*prob5_tree - fail_cost*(1-prob5_tree)
profit_tree <- cbind(profit1_tree,profit2_tree,profit3_tree,profit4_tree,profit5_tree)

# calls_to_make is the number of calls you plan to make with this strategy
calls_to_make_tree <-rowSums(profit_tree>0)

# success_calls is 0 if you don't call them OR call and don't make a sale; if you experience success on that call, it represents the call number on which you were successful
success_calls_tree <- ifelse(calls_to_make_tree>=df_test$campaign & df_test$Class,df_test$campaign,0)

# Fail_calls is the number of unsuccessful calls you would made during the campaign, if you called each customer based on calls_to_make
fail_calls_tree <- ifelse(calls_to_make_tree>=df_test$campaign,ifelse(!df_test$Class,df_test$campaign,df_test$campaign-1),calls_to_make_tree)
```

## Function for calculating expected profit per call - Random Forest

```{r}
sample_tree <-df_test
          
#Loops through prediction turning on campaigns to calculate probabilities of success on calls 1:5

sample_tree$campaign2<-0
sample_tree$campaign3<-0
sample_tree$campaign4<-0
sample_tree$campaign5<-0
rf_test<-sample_tree[,c(56,1:30,47:49, 58:61)]
prob1_rf <- predict(fit2,newdata=rf_test,type="prob")[,2]

rf_test$campaign2<-1
prob2_rf <- predict(fit2,newdata=rf_test,type="prob")[,2]

rf_test$campaign2<-0
rf_test$campaign3<-1
prob3_rf <- predict(fit2,newdata=rf_test,type="prob")[,2]

rf_test$campaign3<-0
rf_test$campaign4<-1
prob4_rf <- predict(fit2,newdata=rf_test,type="prob")[,2]

rf_test$campaign4<-0
rf_test$campaign5<-1
prob5_rf <- predict(fit2,newdata=rf_test,type="prob")[,2]

prob_rf<-cbind(prob1_rf,prob2_rf,prob3_rf,prob4_rf,prob5_rf)

#For Calls 1 through N, compute likelihood of success by the nth call: p_sub(n)
        
p_sub1_rf <- prob_rf[,1]
p_sub2_rf <- p_sub1_rf+prob_rf[,2]*(1-p_sub1_rf)
p_sub3_rf <- p_sub2_rf+prob_rf[,3]*(1-p_sub2_rf)
p_sub4_rf <- p_sub3_rf+prob_rf[,4]*(1-p_sub3_rf)
p_sub5_rf <- p_sub4_rf+prob_rf[,5]*(1-p_sub4_rf)

  
#For calls 1 through N, compute lift of call n: lift(n)= p_sub(n) - p_sub(n-1) (this represents unconditional probability of success on call n)

lift1_rf <- p_sub1_rf
lift2_rf <- p_sub2_rf-p_sub1_rf
lift3_rf <- p_sub3_rf-p_sub2_rf
lift4_rf <- p_sub4_rf-p_sub3_rf
lift5_rf <- p_sub5_rf-p_sub4_rf

#Compute number of calls to target for a given customer (N.B.: As long as probability of success is declining monotonically, as soon as the nth call is not profitable, the (n+1)th call will also not be profitable.) 
  #Test: Identify whether the nth call (if made) would be profitable: Profit(n|All calls prior="no")=[Profit of Success]*p_sub(n) -[Cost of Failure]*(1-p_sub(n)) 

profit1_rf <- Success_profit*prob1_rf - fail_cost*(1-prob1_rf)
profit2_rf <- Success_profit*prob2_rf - fail_cost*(1-prob2_rf)
profit3_rf <- Success_profit*prob3_rf - fail_cost*(1-prob3_rf)
profit4_rf <- Success_profit*prob4_rf - fail_cost*(1-prob4_rf)
profit5_rf <- Success_profit*prob5_rf - fail_cost*(1-prob5_rf)
profit_rf <- cbind(profit1_rf,profit2_rf,profit3_rf,profit4_rf,profit5_rf)

# calls_to_make is the number of calls you plan to make with this strategy
calls_to_make_rf <-rowSums(profit_rf>0)

# success_calls is 0 if you don't call them OR call and don't make a sale; if you experience success on that call, it represents the call number on which you were successful
success_calls_rf <- ifelse(calls_to_make_rf>=df_test$campaign & df_test$Class,df_test$campaign,0)

# Fail_calls is the number of unsuccessful calls you would made during the campaign, if you called each customer based on calls_to_make
fail_calls_rf <- ifelse(calls_to_make_rf>=df_test$campaign,ifelse(!df_test$Class,df_test$campaign,df_test$campaign-1),calls_to_make_rf)
```

## Function for calculating expected profit per call - Boosted Trees

```{r}
sample_tree <- df_test
          
#Loops through prediction turning on campaigns to calculate probabilities of success on calls 1:5
sample_tree$campaign2<-0
sample_tree$campaign3<-0
sample_tree$campaign4<-0
sample_tree$campaign5<-0
gbm_test<-sample_tree[,c(56,1:30,47:49, 58:61)]
prob1_gbm <- log_prob(predict(fboost,newdata=gbm_test,n.trees=500))

gbm_test$campaign2<-1
prob2_gbm <- log_prob(predict(fboost,newdata=gbm_test,n.trees=500))

gbm_test$campaign2<-0
gbm_test$campaign3<-1
prob3_gbm <- log_prob(predict(fboost,newdata=gbm_test,n.trees=500))

gbm_test$campaign3<-0
gbm_test$campaign4<-1
prob4_gbm <- log_prob(predict(fboost,newdata=gbm_test,n.trees=500))

gbm_test$campaign4<-0
gbm_test$campaign5<-1
prob5_gbm <- log_prob(predict(fboost,newdata=gbm_test,n.trees=500))

prob_gbm <- cbind(prob1_gbm,prob2_gbm,prob3_gbm,prob4_gbm,prob5_gbm)

#For Calls 1 through N, compute likelihood of success by the nth call: p_sub(n)
p_sub1_gbm <- prob_gbm[,1]
p_sub2_gbm <- p_sub1_gbm+prob_gbm[,2]*(1-p_sub1_gbm)
p_sub3_gbm <- p_sub2_gbm+prob_gbm[,3]*(1-p_sub2_gbm)
p_sub4_gbm <- p_sub3_gbm+prob_gbm[,4]*(1-p_sub3_gbm)
p_sub5_gbm <- p_sub4_gbm+prob_gbm[,5]*(1-p_sub4_gbm)

#For calls 1 through N, compute lift of call n: lift(n)= p_sub(n) - p_sub(n-1) (this represents unconditional probability of success on call n)

lift1_gbm <- p_sub1_gbm
lift2_gbm <- p_sub2_gbm-p_sub1_gbm
lift3_gbm <- p_sub3_gbm-p_sub2_gbm
lift4_gbm <- p_sub4_gbm-p_sub3_gbm
lift5_gbm <- p_sub5_gbm-p_sub4_gbm

#Compute number of calls to target for a given customer (N.B.: As long as probability of success is declining monotonically, as soon as the nth call is not profitable, the (n+1)th call will also not be profitable.) 
  #Test: Identify whether the nth call (if made) would be profitable. THis is simple:
   #Profit(n|All calls prior="no")=[Profit of Success]*p_sub(n) -[Cost of Failure]*(1-p_sub(n)) 

profit1_gbm <- Success_profit*prob1_gbm - fail_cost*(1-prob1_gbm)
profit2_gbm <- Success_profit*prob2_gbm - fail_cost*(1-prob2_gbm)
profit3_gbm <- Success_profit*prob3_gbm - fail_cost*(1-prob3_gbm)
profit4_gbm <- Success_profit*prob4_gbm - fail_cost*(1-prob4_gbm)
profit5_gbm <- Success_profit*prob5_gbm - fail_cost*(1-prob5_gbm)
profit_gbm <- cbind(profit1_gbm,profit2_gbm,profit3_gbm,profit4_gbm,profit5_gbm)

# calls_to_make is the number of calls you plan to make with this strategy
calls_to_make_gbm <-rowSums(profit_gbm>0)

# success_calls is 0 if you don't call them OR call and don't make a sale; if you experience success on that call, it represents the call number on which you were successful
success_calls_gbm <- ifelse(calls_to_make_gbm>=df_test$campaign & df_test$Class,df_test$campaign,0)

# Fail_calls is the number of unsuccessful calls you would made during the campaign, if you called each customer based on calls_to_make
fail_calls_gbm <- ifelse(calls_to_make_gbm>=df_test$campaign,ifelse(!df_test$Class,df_test$campaign,df_test$campaign-1),calls_to_make_gbm)
```

# Section 4: Assess how profitable each model's strategy would have actually been (on the test data)

# These blocks of code simply call the functions above for various models
# Some portions of the code are slow to run so outputs have been saved in the folder and simply loaded in

```{r}
# #Predict decisions (max calls to make) for test dataset using lasso
# 
# lasso_results <- data.frame(Max_calls_lasso = rep(0,nrow(df_test)), p1_lasso = rep(0,nrow(df_test)), p2_lasso = rep(0,nrow(df_test)), p3_lasso = rep(0,nrow(df_test)), p4_lasso = rep(0,nrow(df_test)), p5_lasso = rep(0,nrow(df_test)))
# 
# for (x in 1:nrow(df_test)) {
#   lasso_results[x,] <-Compute_Sample_Profitability_LASSO(x, cv_lasso)
#   }
# 
# df_test <- cbind(df_test,lasso_results)
# 
# hist(df_test$Max_calls_lasso)
# nrow(df_test[df_test$Max_calls_lasso<df_test$campaign])/nrow(df_test)

```

```{r}
# #Predict decisions (max calls to make) for test dataset using elastic net
# 
# enet_results <- data.frame(Max_calls_enet = rep(0,nrow(df_test)), p1_enet = rep(0,nrow(df_test)), p2_enet = rep(0,nrow(df_test)), p3_enet = rep(0,nrow(df_test)), p4_enet = rep(0,nrow(df_test)), p5_enet = rep(0,nrow(df_test)))
# 
# for (x in 1:nrow(df_test)) {
#   enet_results[x,] <-Compute_Sample_Profitability_LASSO(x, cv_enet)
#   }
# 
# df_test <- cbind(df_test,enet_results)
# 
# hist(df_test$Max_calls_enet)
# nrow(df_test[df_test$Max_calls_enet<df_test$campaign])/nrow(df_test)

```

```{r}
# #Predict decisions (max calls to make) for test dataset using ridge
# 
# ridge_results <- data.frame(Max_calls_ridge = rep(0,nrow(df_test)), p1_ridge = rep(0,nrow(df_test)), p2_ridge = rep(0,nrow(df_test)), p3_ridge = rep(0,nrow(df_test)), p4_ridge = rep(0,nrow(df_test)), p5_ridge = rep(0,nrow(df_test)))
# 
# for (x in 1:nrow(df_test)) {
#   ridge_results[x,] <-Compute_Sample_Profitability_LASSO(x, cv_ridge)
#   }
# 
# df_test <- cbind(df_test,ridge_results)
# 
# hist(df_test$Max_calls_ridge)
# nrow(df_test[df_test$Max_calls_ridge<df_test$campaign])/nrow(df_test)

```

```{r}
# #make decision for logit (employment) model
# logit_results <- data.frame(Max_calls_logit = rep(0,nrow(df_test)), p1_logit = rep(0,nrow(df_test)), p2_logit = rep(0,nrow(df_test)), p3_logit = rep(0,nrow(df_test)), p4_logit = rep(0,nrow(df_test)), p5_logit = rep(0,nrow(df_test)))
# 
# suppressWarnings( for (x in 1:nrow(df_test)) {
#   logit_results[x,] <-Compute_Sample_Profitability(x)
#   }
# )
# 
# df_test <- cbind(df_test,logit_results)
# 
# hist(df_test$Max_calls_logit)
# nrow(df_test[df_test$Max_calls_logit<df_test$campaign])/nrow(df_test)

```

```{r}
# # Randomly selects Max_calls in range 0:5

# set.seed(99)
# df_test$Max_calls_random <- sample(0:5,nrow(df_test),replace=TRUE)
```

```{r}
# save(df_test,file='df_test.rda')

load(file='df_test.rda')
```

# Compute profitability conditional on max number of calls expected to be profitable prediction given by model/decision algorithm

```{r}
#compute profit of calling customers randomly
random_profit <-0
for (x in 1:nrow(df_test)) {
  if  (df_test$campaign[x]>df_test$Max_calls_random[x]){
       random_profit<-random_profit-df_test$Max_calls_random[x]*fail_cost 
      } else {
                if(df_test$Class[x]=="FALSE"){ 
                  random_profit<-random_profit-df_test$campaign[x]*fail_cost
                  } else { random_profit<-random_profit+Success_profit-(df_test$campaign[x]-1)*fail_cost}
              }
}
random_profit

#compute profit of calling all (as called)
baseline_profit <-0
for (x in 1:nrow(df_test)) {
  if(df_test$Class[x]=="FALSE") { baseline_profit<-baseline_profit-df_test$campaign[x]*fail_cost
                                } else { baseline_profit<-baseline_profit+Success_profit-(df_test$campaign[x]-1)*fail_cost}
  
}
baseline_profit

two_call_profit <-0
for (x in 1:nrow(df_test)) {
  if  (df_test$campaign[x]>2){
      two_call_profit<-two_call_profit-2*fail_cost 
      } else {
                if(df_test$Class[x]=="FALSE"){ 
                  two_call_profit<-two_call_profit-df_test$campaign[x]*fail_cost
                  } else { two_call_profit<-two_call_profit+Success_profit-(df_test$campaign[x]-1)*fail_cost}
              }

}

three_call_profit <-0
for (x in 1:nrow(df_test)) {
  if  (df_test$campaign[x]>3){
      three_call_profit<-three_call_profit-3*fail_cost 
      } else {
                if(df_test$Class[x]=="FALSE"){ 
                  three_call_profit<-three_call_profit-df_test$campaign[x]*fail_cost
                  } else { three_call_profit<-three_call_profit+Success_profit-(df_test$campaign[x]-1)*fail_cost}
              }

}

#Compute profit according to logit model
logit_profit <-0
for (x in 1:nrow(df_test)) {
  if  (df_test$campaign[x]>df_test$Max_calls_logit[x]){
      logit_profit<-logit_profit-df_test$Max_calls_logit[x]*fail_cost 
      } else {
                if(df_test$Class[x]=="FALSE"){ 
                  logit_profit<-logit_profit-df_test$campaign[x]*fail_cost
                  } else { logit_profit<-logit_profit+Success_profit-(df_test$campaign[x]-1)*fail_cost}
              }

}
logit_profit

#Compute profit of calling according to LASSO model
lasso_profit <-0
for (x in 1:nrow(df_test)) {
  if  (df_test$campaign[x]>df_test$Max_calls_lasso[x]){
      lasso_profit<-lasso_profit-df_test$Max_calls_lasso[x]*fail_cost 
      } else {
                if(df_test$Class[x]=="FALSE"){ 
                  lasso_profit<-lasso_profit-df_test$campaign[x]*fail_cost
                  } else { lasso_profit<-lasso_profit+Success_profit-(df_test$campaign[x]-1)*fail_cost}
              }

}
lasso_profit

#Compute profit according to GBM model
gbm_profit <- -fail_cost*sum(fail_calls_gbm) + Success_profit*sum(success_calls_gbm>0)

#compute "theoretical best profit" ->targeting only people who say yes
best_profit <-0
for (x in 1:nrow(df_test)) {
  if(df_test$Class[x]=="TRUE"){  
      best_profit<-best_profit+Success_profit-(df_test$campaign[x]-1)*fail_cost
              }

}
NROW(df_test[df_test$Class=="TRUE"])
best_profit


#Compute profit - ridge regression
ridge_profit <-0
for (x in 1:nrow(df_test)) {
  if  (df_test$campaign[x]>df_test$Max_calls_ridge[x]){
       ridge_profit<-ridge_profit-df_test$Max_calls_ridge[x]*fail_cost 
      } else {
                if(df_test$Class[x]=="FALSE"){ 
                  ridge_profit<-ridge_profit-df_test$campaign[x]*fail_cost
                  } else { ridge_profit<-ridge_profit+Success_profit-(df_test$campaign[x]-1)*fail_cost}
              }

}
ridge_profit

#Compute profit - elastic net
enet_profit <-0
for (x in 1:nrow(df_test)) {
  if  (df_test$campaign[x]>df_test$Max_calls_enet[x]){
      enet_profit<-enet_profit-df_test$Max_calls_enet[x]*fail_cost 
      } else {
                if(df_test$Class[x]=="FALSE"){ 
                  enet_profit<-enet_profit-df_test$campaign[x]*fail_cost
                  } else { enet_profit<-enet_profit+Success_profit-(df_test$campaign[x]-1)*fail_cost}
              }

}
enet_profit

#Compute profit of regression tree model
regression_tree_profit <- -fail_cost*sum(fail_calls_tree) + Success_profit*sum(success_calls_tree>0)

#Compute profit of random forest model
rf_profit <- -fail_cost*sum(fail_calls_rf) + Success_profit*sum(success_calls_rf>0)


```

# Section 5. Compare models

```{r}
#Plot projected profits on test set of applying each of the strategies/models
profit_plot <- barplot(c(random_profit,baseline_profit,two_call_profit,three_call_profit,logit_profit,lasso_profit, gbm_profit), names = c("random","baseline","two-call","three-call","logit","lasso","gbm"),ylim=c(0,10000), main="Profits by Strategy", xlab="Strategy", ylab="Profits", col = c('gray','gray','gray','gray','green','blue','red'))
text(profit_plot,c(random_profit,baseline_profit,two_call_profit,three_call_profit,logit_profit,lasso_profit, gbm_profit)+300, labels=format(c(random_profit,baseline_profit,two_call_profit,three_call_profit,logit_profit,lasso_profit, gbm_profit),digits=0))
```

# Generate profit curve of prediction models for evaluation of targeting
```{r}
#Generate profit curve from 1 to 100% targeted using logit model
model_curve_logit <- data.frame(n = seq(from = .01, to= 1, by = .01), logit_profit = 0)

#Generate probability of conversion for the call prior to the last one observed in the dataset
for (n in 1:nrow(df_test)){
  df_test_profit <- df_test[n, current_prob := get(paste0("p",df_test$campaign[n],"_logit"))]
}
#Rank calls according to their probability of conversion
df_test_profit[,rank := frank(current_prob)]
df_test_profit <- df_test_profit[order(-rank)]

#Calculate profits at each point in the curve using probabilities generated above
for (p in seq(from = .01, to= 1, by = .01)){
  num<-round(nrow(df_test_profit)*p)
  model_curve_logit[p*100,2] <- df_test_profit[1:num, sum(Class)]*Success_profit - (num-df_test_profit[1:num, sum(Class)])*fail_cost
}

#For some reason the loop breaks while running so I had to do this part twice.
for (p in c(.06,.07)){
  num<-round(nrow(df_test_profit)*p)
  model_curve_logit[p*100,2] <- df_test_profit[1:num, sum(Class)]*Success_profit - (num-df_test_profit[1:num, sum(Class)])*fail_cost
}

```

```{r}
#Same as above but for LASSO
model_curve_lasso <- data.frame(n = seq(from = .01, to= 1, by = .01), lasso_profit = 0)

for (n in 1:nrow(df_test)){
  df_test_profit <- df_test[n, current_prob := get(paste0("p",df_test$campaign[n],"_lasso"))]
}
df_test_profit[,rank := frank(current_prob)]
df_test_profit <- df_test_profit[order(-rank)]

for (p in seq(from = .01, to= 1, by = .01)){
  num<-round(nrow(df_test_profit)*p)
  model_curve_lasso[p*100,2] <- df_test_profit[1:num, sum(Class)]*Success_profit - (num-df_test_profit[1:num, sum(Class)])*fail_cost
}

for (p in c(.06,.07)){
  num<-round(nrow(df_test_profit)*p)
  model_curve_lasso[p*100,2] <- df_test_profit[1:num, sum(Class)]*Success_profit - (num-df_test_profit[1:num, sum(Class)])*fail_cost
}
```

```{r}
#Same as above but for GBM
model_curve_gbm <- data.frame(n = seq(from = .01, to= 1, by = .01), gbm_profit = 0)

df_test_profit <- df_test

for (n in 1:nrow(df_test)){
  df_test_profit[n,'current_prob'] <- as.numeric(prob_gbm[n,df_test$campaign[n]])
}
df_test_profit <- df_test_profit[,rank := frank(df_test_profit$current_prob)]
df_test_profit <- df_test_profit[order(-rank)]

for (p in seq(from = .01, to= 1, by = .01)){
  num<-round(nrow(df_test_profit)*p)
  model_curve_gbm[p*100,2] <- df_test_profit[1:num, sum(Class)]*Success_profit - (num-df_test_profit[1:num, sum(Class)])*fail_cost
}

for (p in c(.06,.07)){
  num<-round(nrow(df_test_profit)*p)
  model_curve_gbm[p*100,2] <- df_test_profit[1:num, sum(Class)]*Success_profit - (num-df_test_profit[1:num, sum(Class)])*fail_cost
}
```

```{r}
#This time the dataset is sorted randomly instead of based on a prediction model. Otherwise the code does the same thing. 
model_curve_random <- data.frame(n = seq(from = .01, to= 1, by = .01), random_profit = 0)

random_order <- sample(1:nrow(df_test),size=nrow(df_test),replace=FALSE)
df_test_profit <- df_test[order(random_order)]

for (p in seq(from = .01, to= 1, by = .01)){
  num<-round(nrow(df_test_profit)*p)
  model_curve_random[p*100,2] <- df_test_profit[1:num, sum(Class)]*Success_profit - (num-df_test_profit[1:num, sum(Class)])*fail_cost
}

for (p in c(.06,.07)){
  num<-round(nrow(df_test_profit)*p)
  model_curve_random[p*100,2] <- df_test_profit[1:num, sum(Class)]*Success_profit - (num-df_test_profit[1:num, sum(Class)])*fail_cost
}

```

```{r}
#Combine all the curves above into one dataset
model_curve <- cbind(model_curve_logit,lasso_profit = model_curve_lasso$lasso_profit, gbm_profit = model_curve_gbm$gbm_profit, random_profit = model_curve_random$random_profit)

#Plot the curves on one graph
ggplot(data=model_curve) + geom_line(aes(x=n,y=lasso_profit,col='LASSO')) + geom_line(aes(x=n,y=logit_profit, col = 'Logit')) + geom_line(aes(x=n,y=random_profit, col = 'Random')) + geom_line(aes(x=n,y=gbm_profit, col = 'GBM')) + labs(title="Total Profits by Percent Targeted", x = "Percent Targeted", y = "Profits", color = "Model") + scale_color_manual(values=c("LASSO"="blue", "GBM"="red", "Logit"="green", "Random"="black")) + theme(plot.title = element_text(hjust = 0.5))
```

# Output interface of prioritized list of leads for telemarketers

```{r}
#Load original data with categorical variables instead of model matrix
raw_data <- read.csv("bank_additional_full.csv")
raw_data <- as.data.table(raw_data)

#Add customer IDs so I can merge probabilities later
interface <- raw_data[,cust_id := .I]

#Keep only those customers who have not converted and output their age, job, marital, education, and campaign statuses. These are variables that are likely important for telemarketers
interface <- interface[y=="false",c('cust_id','age','job','marital','education','campaign')]

#Merge propensities from selected LASSO model back onto dataset and sort the data by propensity 
interface <- merge(x=interface, y=df_test[,c('current_prob','Max_calls_lasso','ID')], by.x='cust_id', by.y='ID', all.x = FALSE, all.y = FALSE)
interface[,priority := frank(-current_prob,ties.method="first")]

#Clean up variables names and order
names(interface)[c(6,8)] <- c('campaign_calls_received','maximum_times_to_call')
interface <- interface[order(priority),c(9,1:6)]

#Output final interface for top 10 priority leads
head(interface,10)
```

# Section 6. k Nearest Neighbors
```{r}
### Lookalike modeling - kNN. This section uses the same training and testing samples as other modeling approaches. Neighbor votes from 1 to 20 neighbors, as defined by Euclidean distance, are explored.  We also explored using Mahalanobis (i.e., variance-weighted) distances, but the computing power required to do so made implementation impractical. As shown in our slides and executive summary, kNN is not an effective approach at this scale. However, the approach yieled high specificity--but low target rates--at high neighbor counts, suggesting kNN may be appropriate with larger prospect lists, such as those enabled by social media marketing channels.

TrnX <- df_train[,1:30]
OrigTrnG <- df_train$Class

TstX <- df_test[,1:30]

# nnresult <- data.frame("Neighbors" = NULL, "Total_Calls" = NULL, "Good_Calls" = NULL, "Total_Profit" = NULL)
# for (i in 1:20){
#   nn <- knn(TrnX, TstX, OrigTrnG, k=i, l = 0, prob=FALSE, use.all = FALSE)
#   nn <- cbind(nn,df_test[,56])
# 
#   nn$call <- 0
#   nn$call[nn$nn==TRUE] <- 1
# 
#   nn$match <- nn$nn == nn$Class & nn$Class== TRUE
#   nn$profit <- 0
#   nn$profit[nn$match==TRUE] <- Success_profit
# 
#   nn$failcall <- nn$nn == TRUE & nn$Class == FALSE
#   nn$failcost <-0
#   nn$failcost[nn$failcall==TRUE] <- fail_cost
# 
#   thisresult <- data.frame("Neighbors" = i ,
#                            "Total_Calls" = sum(nn$call),
#                            "Good_Calls" = length(nn$match[nn$match==TRUE]),
#                            "Total_Profit" = sum(nn$profit) - sum(nn$failcost)
#                            )
# 
# 
#   nnresult <- rbind(nnresult, thisresult)
# }
# save(nnresult, file="nnresult.rda")
load(file="nnresult.rda")
write.csv(nnresult, "knn_result.csv")
nnresult

```