FLM2_Bertrand_valueOfStructure.R

##########################################################################
##	Project: FLM2 -- Farrell, Liang, Misra: arXiv:2010.14694 
##	Purpose: Demonstrate value of structural models using Bertrand data
##	Author: Sanjog Misra & Max H. Farrell
##	Date: 2025-04-21
##
## Simple example showing the value of structure
##   - estimate demand functions using forests, deep nets, & splines
##   - plot results
##   - show that optimization fails
##   - use a simple structural logit and show that it works
##
##	Data:
##		- adcontentworth_qjecsv.csv from Bertrand et al, QJE 2010 125 (1):263–306
##
##########################################################################

rm(list=ls(all.names = TRUE))
library(ranger)
library(binsreg)
library(ggplot2)
library(torch)




##########################################################################
## DNN function for least squares  regression: y = f(x) + e
linDNN <- function(x=x,y=y,arch=c(20,20),NEPOCH=1000,.seed=1234){
  # input dimension (number of input features)
  d_in <- ncol(x)
  # Number of outcomes
  d_out <- ncol(y) 
  # number of observations in full data set
  n <- nrow(x)
  
  # Set a seed (mostly for replication)
  set.seed(.seed)
  torch_manual_seed(.seed)
  
  # Model constructor
  model <- nn_sequential(nn_linear(d_in, arch[1]),nn_relu())
  j=1 # in case single layer
  # Loop through architecture 
  if(length(arch)>1){
    for(j in 2:length(arch)){
      model$add_module(name = paste("layer_",j-1),module=nn_linear(arch[j-1],arch[j]))  
      model$add_module(name = paste0("relu_",j-1),nn_relu())  
    }
  }
  
  # (Parameter) Output layer
  model$add_module("xb",nn_linear(arch[j],d_out))  
  
  # Learning framework
  # for ADAM, need to choose a reasonable learning rate
  learning_rate <- 0.01
  optimizer <- optim_adam(model$parameters, lr = learning_rate)
  
  # Initialize
  for(i in seq_along(model$parameters))
  {
    nn_init_normal_(model$parameters[[i]],0,.1)
  }
  
  # Training Loop
  intv = NEPOCH/100
  cat("\n Begining training...\n")
  pb <- txtProgressBar(min=1,max=100,style=3)
  pct=0
  loss.stor = NULL
  st = proc.time()
  
  # Training loop
  for (t in 1:NEPOCH) {
    ### -------- Forward pass -------- 
    ### Causal Model 
    y_pred <- model(x)
    loss <- torch_mean((y_pred-y)^2)
    # nnf_mse_loss(y_pred, y, reduction = "mean")
    
    # -------- Backpropagation -------- 
    # Need to zero out the gradients before the backward pass
    optimizer$zero_grad()
    
    # gradients are computed on the loss tensor
    loss$backward()
    
    # -------- Update weights -------- 
    # use the optimizer to update model parameters
    optimizer$step()
    
    # progress
    if(t%%intv==0) {pct=pct+1; setTxtProgressBar(pb, pct)}
    
    loss.stor = c(loss.stor,as.numeric(loss)) 
  }
  
  list(model=model,yhat=y_pred,loss.stor=loss.stor)
  
}





##########################################################################
## Data preparation

# Read CSV File
adc <- read.csv("adcontentworth_qjecsv.csv")
# Only using Wave>1
adc <- adc[adc$wave>1 & adc$risk=="HIGH",]
# Loans applied
table(adc$applied)

# Change offer to be in (0,1)
adc$offer4 = adc$offer4/100
term = 4 # loan term
rec = 0 # recovery % on default 

# collapse data for plotting, average of application 0/1 for each value of interest rate
df2 = aggregate(adc$applied~adc$offer4,FUN=mean)


png('output/Fig0_SimpleEG_dataOnly.png')
  plot(x=df2$`adc$offer4`,y=df2$`adc$applied`, pch=19, xlim=c(0.025,.125), xlab="Interest Rate", ylab="Application Probability")
  legend("topright", legend = c("Average Outcome"), lty=c(NA), lwd=3, col=c("black"), pch=c(19))
dev.off()





##########################################################################
## Demand Estimation

## Random forests
# default (appears undersmoothed)
r1_u = ranger(applied~offer4,data=adc)
p1_u = predict(r1_u,data=adc)
# add more smoothing
r1 = ranger(applied~offer4,data=adc,num.trees = 1000,max.depth = 1)
p1 = predict(r1,data=adc)

png('output/Fig1a_SimpleEG_forests.png')
  plot(x=df2$`adc$offer4`,y=df2$`adc$applied`, pch=19, xlim=c(0.025,.125), xlab="Interest Rate", ylab="Application Probability")
  lines(p1_u$predictions[order(adc$offer4)]~adc$offer4[order(adc$offer4)],col='darkorange',lwd=3, lty="dotted")
  lines(p1$predictions[order(adc$offer4)]~adc$offer4[order(adc$offer4)],col='blue',lwd=3,lty="dashed")
  legend("topright", legend = c("Average Outcome", "Random Forest - Default", "Random Forest - Smoothed"), lty=c(NA,"dotted","dashed"), lwd=3, col=c("black","darkorange", "blue"), pch=c(19,NA,NA))
dev.off()


## Deep Nets
d1 = linDNN(y=as.matrix(adc$applied),x = as.matrix(adc$offer4),arch=c(10,10),NEPOCH = 5000)
dp1 = as.numeric(d1$model(as.matrix(adc$offer4)))
d2 = linDNN(y=as.matrix(adc$applied),x = as.matrix(adc$offer4),arch=c(20, 50, 20),NEPOCH = 5000)
dp2 = as.numeric(d2$model(as.matrix(adc$offer4)))
png('output/Fig1b_SimpleEG_DNN.png')
  plot(x=df2$`adc$offer4`,y=df2$`adc$applied`, pch=19, xlim=c(0.025,.125), xlab="Interest Rate", ylab="Application Probability")
  lines(dp1[order(adc$offer4)]~adc$offer4[order(adc$offer4)],col='blue',lwd=3, lty="dashed")
  lines(dp2[order(adc$offer4)]~adc$offer4[order(adc$offer4)],col='darkorange',lwd=3, lty="dotted")
  legend("topright", legend = c("Average Outcome", "Neural Network (20,50,20)", "Neural Network (10,10)"), lty=c(NA,"dotted","dashed"), lwd=3, col=c("black","darkorange", "blue"), pch=c(19,NA,NA))
dev.off()





# Logit
plot(x=df2$`adc$offer4`,y=df2$`adc$applied`,pch=19,xlim=c(0.025,.125),xlab="Interest Rate",ylab="Application Probability")
out = glm(applied~offer4,data=adc,family=binomial)
ahat = predict(out,type='response')
lines(ahat[order(adc$offer4)]~adc$offer4[order(adc$offer4)],col='red',lwd=2,type='l')
png('output/Fig1c_SimpleEG_logit.png')
  plot(x=df2$`adc$offer4`,y=df2$`adc$applied`, pch=19, xlim=c(0.025,.125), xlab="Interest Rate", ylab="Application Probability")
  lines(ahat[order(adc$offer4)]~adc$offer4[order(adc$offer4)],col='blue',lwd=2,lty=2)
  legend("topright", legend = c("Average Outcome", "Structural Logit"), lty=c(NA,2), lwd=3, col=c("black", "blue"), pch=c(19,NA))
dev.off()









##########################################################################
## Demand Estimation Use models above to extrapolate to find expected DEMAND by interest rate


## First extrapolate just a little bit, to 20% interest
  
  #New interest rate offers at which to obtain predictions
  limit <- .2  # x100 percent, e.g., 2=200% interest
  dd = adc[sample(1:nrow(adc),size=200),] #only need a few rows for plotting purposes
  dd$offer4 = .1 + runif(nrow(dd))*1.9
  dd$offer4 = runif(nrow(dd))*limit
  
  
  ## Random forests
  rf_extrap = predict(r1,data=dd)
  ## Deep Nets
  dnn_extrap = as.numeric(d1$model(as.matrix(dd$offer4)))
  # Logitic regression
  logit_extrap = predict(out,type='response',newdata = dd)
  
  #extrapolated demand functions
  png('output/Fig1d_SimpleEG_extrapolated_demand_20pct.png')
    plot((logit_extrap)[order(dd$offer4)]~dd$offer4[order(dd$offer4)], type='l', ylab="Application Probability", xlab="Interest Rate", lwd=3, col="black", lty=1)
    lines((rf_extrap$predictions)[order(dd$offer4)]~dd$offer4[order(dd$offer4)], lwd=3, col="darkorange", lty=2)
    lines((dnn_extrap)[order(dd$offer4)]~dd$offer4[order(dd$offer4)], lwd=3, col="darkgreen", lty=3)
    legend("topright", legend = c("Random Forest", "Neural Network", "Structural Logit"), lty=c(2,3,1), lwd=3, col=c("darkorange", "darkgreen","black"))
  dev.off()

  
  
  ## Second extrapolate a lot, to 200%
  
  #New interest rate offers at which to obtain predictions
  limit <- 2  # x100 percent, e.g., 2=200% interest
  dd = adc[sample(1:nrow(adc),size=200),] #only need a few rows for plotting purposes
  dd$offer4 = .1 + runif(nrow(dd))*1.9
  dd$offer4 = runif(nrow(dd))*limit
  
  ## Random forests
  rf_extrap = predict(r1,data=dd)
  ## Deep Nets
  dnn_extrap = as.numeric(d1$model(as.matrix(dd$offer4)))
  # Logitic regression
  logit_extrap = predict(out,type='response',newdata = dd)
  
  #extrapolated demand functions
  png('output/Fig1e_SimpleEG_extrapolated_demand_200pct.png')
    plot((logit_extrap)[order(dd$offer4)]~dd$offer4[order(dd$offer4)], type='l', ylab="Application Probability", xlab="Interest Rate", lwd=3, col="black", lty=1)
    lines((rf_extrap$predictions)[order(dd$offer4)]~dd$offer4[order(dd$offer4)], lwd=3, col="darkorange", lty=2)
    lines((dnn_extrap)[order(dd$offer4)]~dd$offer4[order(dd$offer4)], lwd=3, col="darkgreen", lty=3)
    legend("topright", legend = c("Random Forest", "Neural Network", "Structural Logit"), lty=c(2,3,1), lwd=3, col=c("darkorange", "darkgreen","black"))
  dev.off()
  
  
  
  
###################
## Use models above to extrapolate to find expected REVENUE by interest rate

#rev = demand * (interest rate) * 4 months

png('output/Fig1f_SimpleEG_extrapolated_revenue.png')
  plot((rf_extrap$predictions*4*dd$offer4)[order(dd$offer4)]~dd$offer4[order(dd$offer4)], type='l', lwd=3, col="darkorange", lty=2, ylab="Expected Revenue (per 1$)", xlab="Interest Rate")
  lines((dnn_extrap*4*dd$offer4)[order(dd$offer4)]~dd$offer4[order(dd$offer4)], lwd=3, col="darkgreen", lty=3)
  lines((logit_extrap*4*dd$offer4)[order(dd$offer4)]~dd$offer4[order(dd$offer4)], lwd=3, col="black", lty=1)
  abline(h=0, lty="dotted", col="lightgrey")
  legend("topleft", legend = c("Random Forest", "Neural Network", "Structural Logit"), lty=c(2,3,1), lwd=3, col=c("darkorange", "darkgreen","black"))
dev.off()



##### Confidence band
# above problems are not related to statistical uncertainty

#test monotonicity - not rejected of course
summary(  binstest(y=applied, x=offer4, data=adc, testshapel=0, deriv=1, randcut=1, nsims=5000, simsgrid = 100) )
summary(  binstest(y=applied, x=offer4, data=adc, estmethod="glm", family="binomial", testshapel=0, deriv=1, randcut=1, nsims=5000, simsgrid = 100) )


# Confidence band with the three earlier estimates
plotting.data <- data.frame(x=adc$offer4, rf=p1$predictions, nn=dp1, logit=ahat)
#using least squares (linear probability model)
spline.fit <- binsreg(y=applied, x=offer4, data=adc, cb=TRUE, randcut=1, nsims=5000, simsgrid = 100, dots = c(2,2), dotsgrid = 0, dotsgridmean = 0, plotxrange = range(plotting.data$x))
#logistic regression
spline.fit <- binsglm(y=applied, x=offer4, data=adc, family=binomial(), cb=TRUE, randcut=1, nsims=5000, simsgrid = 100, dots = c(2,2), dotsgrid = 0, dotsgridmean = 0, plotxrange = range(plotting.data$x))

fig <- spline.fit$bins_plot + 
  labs(x='Interest Rate', y='Application Probability') + 
  theme(axis.text=element_text(size=16), 
        axis.title=element_text(size=16), 
        plot.background = element_blank(), 
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        legend.position=c(0.85, 0.9),
        legend.background = element_blank(),
        legend.box.background = element_rect(colour = "black")) + 
  geom_line(data=plotting.data, aes(x=x, y=rf, color="Random Forest", linetype="Random Forest"), size=1.5) +  
  geom_line(data=plotting.data, aes(x=x, y=nn, color="Neural Net", linetype="Neural Net"), size=1.5) +  
  geom_line(data=plotting.data, aes(x=x, y=logit, color="Structural Logit", linetype="Structural Logit"), size=1.5) +  
  scale_color_manual(name='Model',
                     breaks=c('Random Forest', 'Neural Net', 'Structural Logit'),
                     values=c('Random Forest'='darkorange', 
                              'Neural Net'='darkgreen', 
                              'Structural Logit'='black')) +
  scale_linetype_manual(name='Model',
                        breaks=c('Random Forest', 'Neural Net', 'Structural Logit'),
                        values=c('Random Forest'='dashed', 
                                 'Neural Net'='dotted', 
                                 'Structural Logit'='solid')) +
  guides(color = guide_legend("Model"), linetype = guide_legend("Model"))  # Combine legends



png('output/Fig0_SimpleEG_confidence_band.png')
  plot(fig)
dev.off()