FLM2_Bertrand_step2_InferenceStep.R

##########################################################################
##	Project: FLM2 -- Farrell, Liang, Misra: arXiv:2010.14694 
##	Purpose: Analyze loan data from Bertrand et al, QJE 2010 125 (1):263–306
##            -- Step 2: Semiparametric inference
##	Author: Sanjog Misra & Max H. Farrell
##	Date: 2025-04-21
##
## Second stage
##   - Inference on mu = E[H()] (see step 0 file for H functions)
##   - main examples are marginal effects and optimal profits (using fixed point)
##   - completes 50-fold cross fitting (i) loads DNN results from step 1 (ii) estimates Lambda(x)
##
##	Data:
##		- adcontentworth_qjecsv.csv from Bertrand et al, QJE 2010 125 (1):263–306
##		- DNN results for each fold, saved as out.dnn_X.Rd, X=1,...,50 (see step 1)
##
##	Output:
##		- results and figures for the paper
##
##########################################################################

rm(list=ls(all.names = TRUE))

# Load libraries
library(torch)
library(pbmcapply)
library(numDeriv)


##########################################################################
## Load results of the first stage deep net

# each element of out.dnn contains torch model
out.dnn=list()
Nsplits=50
for(j in 1:Nsplits) {
  fname = paste0("output/out.dnn_",j,".Rd")
  out.dnn[[j]] = list()
  out.dnn[[j]]$mod = torch_load(path=fname)
  cat(j,"\r")
}

##########################################################################
## Data preparation
# same as step 1

# Read CSV File
adc = read.csv("adcontentworth_qjecsv.csv")
# Only using Wave>1 and High Risk types
adc = adc[adc$wave>1 & adc$risk=="HIGH",]

# List of variables to be used
vlist = c("offer4","speak_trt","stripany",
          "dphoto_none", "dphoto_black", "dphoto_female",
          "prize", "oneln_trt",
          "use_any", "intshown", "comploss_n", "comp_n")

# Demographics/pre-treatment Covariates
dems = c("waved3","dormancy","trcount","female", "race", "nspeakeligible")

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

# Construct Data
# applied (we assume all applications are approved)
y1 = as.matrix(as.integer(adc$applied),ncol=1) 
table(y1)
# Default
y2 = as.matrix(as.integer((adc$badacct_last)),ncol=1)
table(y2)
# Covariates
x=  model.matrix(~.,data=adc[,dems])
# Remove Intercept
x= x[,-1]
# Treatments (Interest rate offers) in pct
z = model.matrix(~offer4-1,data=adc[,vlist])
# Recode NA as 0 (for loss below)
y2[is.na(y2)]=0

# The Parametric Loss function
parm.loss = function(theta){
  u1 = theta[1]+z*theta[2]
  u2 = theta[3]+z*theta[4]
  
  pr1 = 1/(1+exp(-u1))
  pr2 = 1/(1+exp(-u2))
  
  -sum(y1*log(pr1)+(1-y1)*log(1-pr1)+y1*y2*log(pr2)+y1*(1-y2)*log(1-pr2))
  
}

set.seed(12345)
res = optim(runif(4),parm.loss) #starting values
res = optim(res$par,parm.loss,method="BFGS",hessian=TRUE)
cfs = res$par
se = sqrt(diag(solve(res$hessian)))
tstat = cfs/se
cbind(cfs,se,tstat)

# Equivalently (not used later)
# cfs1 = coef(glm(y1~z,family=binomial))
# cfs2 = coef(glm(y2[y1==1]~z[y1==1],family=binomial))
# cfs=c(cfs1,cfs2)

# Demonstration
profit.par = function(ir){
  term = 4 # Loan duration in Data
  rec = 0 # Recovery percentage 0 means defult is total loss
  u1 = cfs[1]+ir*cfs[2]
  u2 = cfs[3]+ir*cfs[4]
  
  pr1 = 1/(1+exp(-u1))
  pr2 = 1/(1+exp(-u2))
  sum(1000*(pr1*(term*(1-pr2)*ir-pr2)+(1-pr1)*term*.01))
}

profit.par  = Vectorize(profit.par)
curve(profit.par,0,.25)
iopt = optimize(profit.par,maximum = TRUE,c(0,.25))
abline(v=iopt$maximum,lty=2,col='red')





##########################################################################
## Inference steps
#(i) recreate the data and cross fitting flag
#(ii) obtain Lambda(x)^{-1}
#(iii) inference on mu = E[H()]

# Recreate the data from the first stage
# match cross fitting folds by setting seed
y = as.matrix(as.integer(adc$applied),ncol=1)
set.seed(12345)
Nsplits = 50
idx = sample(1:Nsplits,size=length(y),replace = TRUE,prob=rep(1/Nsplits,Nsplits))
table(idx)
# The extra part that a gets (ad attributes)
xa=  model.matrix(~.,data=adc[,c(vlist)])
xa= xa[,-c(1,2)]

# The common parts (demographics)
xb = model.matrix(~.,data=adc[,c(dems)])
xb = xb[,-1]
# z = model.matrix(as.formula(paste0("~",paste0(vlist,sep="",collapse="+"),"-1")),data=adc)                 
z = model.matrix(~offer4-1,data=adc[,vlist])

## looping step
regmods = list()
#apply first step results from DNN, then get Lambda(x)
for(iter in 1:Nsplits) {
  cat("Split ",iter," processing...\n")
  tsmp = (idx!=iter)
  mdat = list(Xa=xa[tsmp,],Xb=xb[tsmp,],Y=as.matrix(y[tsmp],ncol=1),Z=as.matrix(z[tsmp,],ncol=ncol(z)))
  Y <- as.matrix(mdat$Y)
  Z <- as.matrix(mdat$Z)
  Xa <- as.matrix(mdat$Xa)
  Xb <- as.matrix(mdat$Xb)
  n = length(Y)
  
  # Get Model (for a,b )
  model = out.dnn[[iter]]$mod
  #alpha(xa,xb)
  a_out <- model$a_layers(torch_cat(list(Xa, Xb), dim = 2))
  #beta(xb)
  b_out <- model$b_layers(Xb)
  
  ahat=as.matrix(a_out)
  bhat=as.matrix(b_out)
  ab = cbind(ahat,bhat)
  vz = ahat+rowSums(bhat*Z)
  G = 1/(1+exp(-vz))
  Gdot = G*(1-G)
  
  # Construct prediction targets for Lambda(x), elements of Gdot* Trt Trt'  
  l11 = Gdot # Gdot*1
  l12 = l21 = Gdot*Z
  l22 = Gdot*Z^2
  
  # Project on all "Xs"
  xall = cbind(Xa,Xb)
  xall_df <- as.data.frame(xall)
  colnames(xall_df) <- paste0("x", 1:ncol(xall_df))
  
  # Linear Model for each entry of Lambda(x), fit on all other folds
    frml = as.formula(paste0("y~(",paste0("x",1:19,collapse="+"),")*(",paste0("x",1:19,collapse="+"),")"))
    
    newd = data.frame(y=l11,xall_df)
    mmx = model.matrix(frml,newd)[,-1]
    r11 = lm(frml,data = newd)
  
    newd = data.frame(y=l12,xall_df)
    mmx = model.matrix(frml,newd)[,-1]
    r12 = lm(frml,data = newd)
  
    newd = data.frame(y=l22,xall_df)
    mmx = model.matrix(frml,newd)[,-1]
    r22 = lm(frml,data = newd)
    
  regmods[[iter]]=list(r11=r11,r12=r12,r22=r22)
}



# Now on to inference

# if needed
# Rule Function to track is defined globally
 rule.fn <- function(i, Xb_local) {
   sum(Xb_local[i, 5:7] == 1)
 }

# Define function for get IF etc for a given H
getInfObj = function(H){
  for(iter in 1:Nsplits) {
    cat("Split ",iter," processing...\n")
    # Compute influence function adjustment term (apart from H_theta) 
    # Selected the held-out data
    tsmp = (idx==iter)
    mdat = list(Xa=xa[tsmp,],Xb=xb[tsmp,],Y=as.matrix(y[tsmp],ncol=1),Z=as.matrix(z[tsmp,],ncol=ncol(z)))
    Xa = mdat$Xa
    Xb = mdat$Xb
    Z = mdat$Z
    Y = mdat$Y
    # Get Model (for a,b ) for new data
    # Fixing Xa = 0
    model = out.dnn[[iter]]$mod
    a_out_split <- model$a_layers(torch_cat(list(Xa*0, Xb), dim = 2))
    # Compute b(xb)
    b_out_split <- model$b_layers(Xb)
    ahat=as.matrix(a_out_split)
    bhat=as.matrix(b_out_split)
    ab = cbind(ahat,bhat)
    
    # Retain only obs with negative betas
    cnd = ab[,2]< -.1
    Xa = mdat$Xa[cnd,]
    Xb = mdat$Xb[cnd,]
    Z = mdat$Z[cnd,]
    Y = mdat$Y[cnd]
    ab = ab[cnd,]
    ahat=ab[,1]
    bhat=ab[,2]
    xall=cbind(Xa*0,Xb)
    n = length(Y)
    
    # Load Lambda models
    r11 <- regmods[[iter]]$r11
    r12 <- regmods[[iter]]$r12
    r22 <- regmods[[iter]]$r22
    
    # Compute Lam components
    newxall <- as.data.frame(xall)
    colnames(newxall) <- paste0("x", 1:ncol(newxall))
    
    l11_hat = predict(r11,newdata=newxall)
    l12_hat = predict(r12,newdata=newxall)
    l22_hat = predict(r22,newdata=newxall)
    
    # Combine and invert each E(Lam). Use regularized inverse to avoid small denominator problem
    lmat = cbind(l11_hat,l12_hat,l12_hat,l22_hat)
    L = lapply(1:nrow(lmat),FUN=function(i) matrix(lmat[i,],2,2))  
    linv = function(L,lambda=1E-08){
      solve(t(L)%*%L + lambda*diag(nrow(L)))%*%t(L)
    }
    Lami = lapply(L,linv)
    # Get Gradient
    vz = ahat+bhat*Z
    G = 1/(1+exp(-vz))
    
    IFadj = lapply(1:length(Z),function(i) Lami[[i]]%*%((Y[i]-G[i])*c(1,Z[i])))
    
    # Bind ahat and bhat (and any other relevant data if needed)  
    abx = ab
    
    # Apply function to ahat and bhat estimates
    tst = pbmclapply(data.frame(t(abx)),H)
    Hx = matrix(unlist(tst),nrow=length(tst),byrow=TRUE)
    
    # Index update
    its = rep(iter,length(tst))
    
    # Gradient
    hab = function(ab) jacobian(H,ab)
    tst = pbmclapply(data.frame(t(abx)),hab,mc.cores = detectCores()-2)
    habx = tst
    Hd = habx
    
    IFi = function(i){
      c(Hd[[i]]%*%IFadj[[i]])
    }
    
    IF = (sapply(1:length(Y),FUN=IFi))
    if(!is.null(dim(IF))) IF=t(IF)
    
    if(!is.null(rule.fn))  rule = sapply(1:length(Y),FUN=rule.fn,Xb_local=Xb)
    
    if(iter==1) {
      N.stk =Xa.mn =Xb.mn=Z.mn=Y.mn=ab.stack = it.stack= NULL
      rule.stack  =Hx.stack = IF.stack = list()
    } 
    
    Hx.stack[[iter]] = Hx
    IF.stack[[iter]] = IF
    it.stack = c(it.stack,its)
    ab.stack = rbind(ab.stack,ab)
    Xa.mn = rbind(Xa.mn,colMeans(Xa))
    Xb.mn = rbind(Xb.mn,colMeans(Xb))
    Z.mn = c(Z.mn,mean(Z))
    Y.mn = c(Y.mn,mean(Y))
    N.stk = c(N.stk,n)
    rule.stack[[iter]]=rule
  }
  res = list(Hx.stack=Hx.stack,IF.stack=IF.stack,
             it.stack=it.stack,ab.stack=ab.stack,
             Xa.mn=Xa.mn,Xb.mn=Xb.mn,Z.mn = Z.mn,Y.mn=Y.mn,N.stk=N.stk,rule.stack=rule.stack)
}





##########################################################################
## Final inference computations
## results and plots for the paper

# code below is copy and pasted for different H functions. For new results
# define a new H function right here or get from FLM2_Bertrand_step0_H_functions.R




## 1 -- Marginal effect
H = function(abi){
  istar = 0.08392
  pr1 = 1/(1+exp(-abi[1]-abi[2]*istar))
  abi[2]*(1-pr1)*pr1
}   

# Do computation for this H function
res = getInfObj(H)
ab.stack = res$ab.stack
Hx.stack = res$Hx.stack
IF.stack = res$IF.stack
it.stack = res$it.stack
rule.stack = res$rule.stack

# Construct estimate and confidence interval
Hest.plug  = sapply(1:Nsplits,function(i) mean(Hx.stack[[i]]))
Hest  = sapply(1:Nsplits,function(i) mean(Hx.stack[[i]]+IF.stack[[i]]) )
Hvar  = sapply(1:Nsplits,function(i) var(Hx.stack[[i]]+IF.stack[[i]])/length(Hx.stack[[i]]) )
Hse = sqrt(mean(Hvar))
CI = c(mean(Hest)-1.96*Hse,mean(Hest)+1.96*Hse) 
cat(c(Hest=mean(Hest),CI=CI))

# Plotting
# include plugin density, mean, and IF-Adjusted Mean and CI
  talpha=0.01
  Hk = do.call(rbind,Hx.stack)
  Ik = do.call(c,IF.stack)
  qk = quantile(Hk,c(talpha/2,1-talpha/2))
  densplot=function(w,clr='white',bclr='black',add=FALSE,do.log=FALSE,...){
    dw = density(w)
    if(do.log){
      y <- log(w)
      g <- density(y)
      xgrid <- exp(g$x)
      g$y <- c(0, g$y/xgrid)
      g$x <- c(0, xgrid)
      dw=g
    }
    if(!add) plot(dw,col=bclr,...)
    polygon(dw,col = clr,border=bclr)
  }
  
pdf("output/Figure_Bertrand_MarginalEffect.pdf", width = 9, height = 6)
  xlab = "Marginal Effect of Interest Rate"
  densplot(Hk[Hk<qk[2] & Hk>qk[1]],clr="#00000030",bclr = 1,
           main="",xlab=xlab)
  abline(v=mean(Hest.plug),col='black',lty=1,lwd=1.5)
  segments(CI[1],0,CI[2],0,col='black',lwd=5)  
  polygon(x=c(CI[1],CI[2],CI[2],CI[1]),y=c(0,0,1000,1000),col="#00000050",border = F)  
  cat(c(Hest=mean(Hest),CI=CI))
  abline(v=mean(Hest),lty=2,lwd=1.5,col='black')
  
  Hse = sqrt(sapply(Hx.stack,function(z) var(z)/length(z)))
  CI=c(mean(Hk)-1.96*mean(Hse),mean(Hk)+1.96*mean(Hse))
  segments(CI[1],0,CI[2],0,col='black',lwd=5)  
  polygon(x=c(CI[1],CI[2],CI[2],CI[1]),y=c(0,0,1000,1000),col="#00000050",border = F)  
  legend('topleft',lty=c(1,2),legend=c("Plug-in Mean","Adjusted Mean"))
dev.off()






## 2 -- Optimal Price
H = function(abi){ 
  L=1000
  term = 4
  rec=0.0
  pi.i = function(ir){
    pr1 = 1/(1+exp(-abi[1]-abi[2]*ir))
    u2 = cfs[3]+ir*cfs[4] # This is from wayup in the begining
    pr2 = 1/(1+exp(-u2))
    sum((pr1*(term*(1-pr2)*ir-(1-rec)*pr2)+(1-pr1)*term*.01))
  }
  set.seed(12345)
  istar = (optimize(pi.i,c(0,.5),maximum = TRUE))$max
  istar
}

# Do computation for this H function
res = getInfObj(H)
ab.stack = res$ab.stack
Hx.stack = res$Hx.stack
IF.stack = res$IF.stack
it.stack = res$it.stack
rule.stack = res$rule.stack

# Construct estimate and confidence interval
Hest.plug  = sapply(1:Nsplits,function(i) mean(Hx.stack[[i]]))
Hest  = sapply(1:Nsplits,function(i) mean(Hx.stack[[i]]+IF.stack[[i]]) )
Hvar  = sapply(1:Nsplits,function(i) var(Hx.stack[[i]]+IF.stack[[i]])/length(Hx.stack[[i]]) )
Hse = sqrt(mean(Hvar))
CI = c(mean(Hest)-1.96*Hse,mean(Hest)+1.96*Hse) 
cat(c(Hest=mean(Hest),CI=CI))

# Plotting
# Pretty Plot for Fixed points
pdf("output/Figure_Bertrand_OptimalPrice.pdf", width = 9, height = 6)
L=1000
term = 4
rec=0.0
irmax = 0.25 # The max IR
irmin = 0.01 # min
abi = ab[1,] # Placeholder
prr = function(ir){
  pr1 = 1/(1+exp(-abi[1]-abi[2]*ir))
  u2 = cfs[3]+ir*cfs[4] # This is from wayup in the begining
  pr2 = 1/(1+exp(-u2))
  L*((pr1*(term*(1-pr2)*ir-(1-rec)*pr2)+(1-pr1)*term*.01))
}

pp = optimize(prr,c(irmin,irmax),maximum = TRUE)

f2 = function(ir) {
  .001*grad(prr,ir)+ir
}

par(mar=c(5, 4, 4, 4) + 0.1)
abi=ab[5,]
curve(f2,from = 0,to=irmax+.01,ylim=c(-.5,.5),type="n",xlim=c(0.02,.25),
      xlab="Interest Rate (LHS)",ylab="Interest Rate (RHS)")

pp = optimize(prr,c(0,irmax),maximum = TRUE)
pstk = c(pp$maximum)

# Fixed Point Plots
set.seed(123)
splt = 50
NN =  100 #max = length(Hx.stack[[splt]])
idx5 = idx[idx==splt]
ab = ab.stack[it.stack==splt,]
for(i in 1:NN) {
  abi=ab[i,]
  prr = function(ir){
    pr1 = 1/(1+exp(-abi[1]-abi[2]*ir))
    u2 = cfs[3]+ir*cfs[4] # This is from way up in the begining
    pr2 = 1/(1+exp(-u2))
    L*((pr1*(term*(1-pr2)*ir-(1-rec)*pr2)+(1-pr1)*term*.01))
  }
  curve(f2,0*irmin+.01,irmax+.01,add=TRUE,col='#0000003A')
  set.seed(12345)
  pp = optimize(prr,c(0,.5),maximum = TRUE)
  points(pp$maximum,pp$maximum,pch=19,col='#000000A0')
  pstk = c(pstk,pp$maximum)
  segments(x0 =pp$maximum,x1 =  pp$maximum,y0=-.4,y1=pp$maximum,col="#0000000A")
}

abline(0,1,col='black')

pd = density(unlist(pstk))
pdy.og = pd$y
pd$y = .2*pd$y/(max(pd$y))-.4
lines(pd$x,pd$y,type='l',col='grey')
polygon(pd,col = "grey",border="black")
rug(pstk,col='grey')

# Add right-side Y-axis
axis(4, at = c(-.4,-35,-.3,-.3,-25,-.2), labels = c(0,5,10,15,20,25))  # Right side Y-axis
mtext(expression("Density (plugin)"), side = 4, line = 2)  # Label
abline(h=-.4)

Hk=do.call(c,Hx.stack)
#pd = density(unlist(Hk))
#pd$y = .2*pd$y/(max(pd$y))-.4
# lines(pd$x,pd$y,type='l',col='black',lty=2)
rug(Hk,col='#00000003')
# Construct Estimator
Hest.plug  = sapply(1:Nsplits,function(i) mean(Hx.stack[[i]]))
Hest  = sapply(1:Nsplits,function(i) mean(Hx.stack[[i]]+IF.stack[[i]]) )
Hvar  = sapply(1:Nsplits,function(i) var(Hx.stack[[i]]+IF.stack[[i]])/length(Hx.stack[[i]]) )
Hse = sqrt(mean(Hvar))
CI = c(mean(Hest)-1.96*Hse,mean(Hest)+1.96*Hse) 
segments(CI[1],-.4,CI[2],-.4,col='black',lwd=5,lty=1)  
polygon(x=c(CI[1],CI[2],CI[2],CI[1]),y=c(-.4,-.4,1000,1000),col="#00000020",border = F)  
segments(mean(Hest.plug),-.4,mean(Hest.plug),4,col='black',lty=1)  
segments(mean(Hest),-.4,mean(Hest),4,col='black',lty=2)  
legend('topright',lty=c(1,2),legend=c("Plug-in Mean","Adjusted Mean"))

Hse = sqrt(sapply(Hx.stack,function(z) var(z)/length(z)))
CI=c(mean(Hk)-1.96*mean(Hse),mean(Hk)+1.96*mean(Hse))
segments(CI[1],-.4,CI[2],-.4,col='black',lwd=5)  
polygon(x=c(CI[1],CI[2],CI[2],CI[1]),y=c(-.4,-.4,1000,1000),col="#00000050",border = F)  

dev.off()






## 3 -- Expected Profits from Personalization
H = function(abi){ 
  L=1000
  term = 4
  rec=0.0
  pi.i = function(ir){
    pr1 = 1/(1+exp(-abi[1]-abi[2]*ir))
    u2 = cfs[3]+ir*cfs[4] # This is from wayup in the begining
    pr2 = 1/(1+exp(-u2))
    sum((pr1*(term*(1-pr2)*ir-(1-rec)*pr2)+(1-pr1)*term*.01))
  }
  set.seed(12345)
  istar = (optimize(pi.i,c(0,.5),maximum = TRUE))$max
  pi.i(istar)
}

# Do computation for this H function
res = getInfObj(H)
ab.stack = res$ab.stack
Hx.stack = res$Hx.stack
IF.stack = res$IF.stack
it.stack = res$it.stack
rule.stack = res$rule.stack

# Construct estimate and confidence interval
Hest.plug  = sapply(1:Nsplits,function(i) mean(Hx.stack[[i]]))
Hest  = sapply(1:Nsplits,function(i) mean(Hx.stack[[i]]+IF.stack[[i]]) )
Hvar  = sapply(1:Nsplits,function(i) var(Hx.stack[[i]]+IF.stack[[i]])/length(Hx.stack[[i]]) )
Hse = sqrt(mean(Hvar))
CI = c(mean(Hest)-1.96*Hse,mean(Hest)+1.96*Hse) 
cat(c(Hest=mean(Hest),CI=CI))

# Plotting
# include plugin density, mean, and IF-Adjusted Mean and CI
talpha=0.01
Hk = do.call(rbind,Hx.stack)
Ik = do.call(c,IF.stack)
qk = quantile(Hk,c(talpha/2,1-talpha/2))
densplot=function(w,clr='white',bclr='black',add=FALSE,do.log=FALSE,...){
  dw = density(w)
  if(do.log){
    y <- log(w)
    g <- density(y)
    xgrid <- exp(g$x)
    g$y <- c(0, g$y/xgrid)
    g$x <- c(0, xgrid)
    dw=g
  }
  if(!add) plot(dw,col=bclr,...)
  polygon(dw,col = clr,border=bclr)
}

pdf("output/Figure_Bertrand_ExpectedProfits.pdf", width = 9, height = 6)
  xlab = "Expected Profits"
  densplot(Hk[Hk<qk[2] & Hk>qk[1]],clr="#00000030",bclr = 1,
           main="",xlab=xlab)
  abline(v=mean(Hest.plug),col='black',lty=1,lwd=1.5)
  segments(CI[1],0,CI[2],0,col='black',lwd=5)  
  polygon(x=c(CI[1],CI[2],CI[2],CI[1]),y=c(0,0,1000,1000),col="#00000050",border = F)  
  cat(c(Hest=mean(Hest),CI=CI))
  abline(v=mean(Hest),lty=2,lwd=1.5,col='black')
  
  Hse = sqrt(sapply(Hx.stack,function(z) var(z)/length(z)))
  CI=c(mean(Hk)-1.96*mean(Hse),mean(Hk)+1.96*mean(Hse))
  segments(CI[1],0,CI[2],0,col='black',lwd=5)  
  polygon(x=c(CI[1],CI[2],CI[2],CI[1]),y=c(0,0,1000,1000),col="#00000050",border = F)  
  legend('topright',lty=c(1,2),legend=c("Plug-in Mean","Adjusted Mean"))
dev.off()