Dynamic moral hazard problem

This notebook solves for the optimal contract in Make-Whole Clauses as Skin in the Game in the dynamic moral hazard case.

# The usual suspects, whether we need them or not

import numpy as np 
import matplotlib.pyplot as plt 

from scipy.linalg import expm 
from scipy.stats import lognorm
from scipy.optimize import fsolve
from scipy.optimize import minimize
from scipy.integrate import quad
import scipy.optimize as opt
from scipy.interpolate import interp1d

import pandas as pd 
import math

# starts parameters

A=7 # agent productivity
y=3 # output when productive
c=0 # prepayment penalty
R=1/0.9 # interest rate
m=1 # payment
bet=0.9 # discount rate

# we also need a distribution for outside option V
# We will start with a log-normal 
loc=0 # location of the distribution 
# (Be very careful here, this is not the same as the loc parameter in scipy! 
# Scipy is very confusing in this respect, location is best handled via the scale parameter, see below
disp=2 # dispersion of the distribution

# let's plot the histogram for a good visual

fig,ax=plt.subplots(1,1)
s=disp;
scale=np.exp(loc)

trunc=lognorm.cdf(1,s,0, scale);

x = np.linspace(0,5, 1000)
ax.plot(x, lognorm.pdf(x, s,0,scale),
       'r-', lw=2, alpha=0.5, label='lognorm pdf')

# plt.title("Distribution of project quality (p) at birth")
plt.xlabel(r'$V$')
plt.ylabel('density')

plt.show()

# First, for warm-up we solve for total surplus
# Integrand when computing TS

def Vintegr(x):
    return x*lognorm.pdf(x,s,0, scale)   

# first-order condition for effort

TS=0
def condA(e):
    return lognorm.cdf(e*(y+bet*TS), s,0,scale)*(y+bet*TS)-2*A*e

# solve for optimal effort
root=fsolve(condA,1)

it=1
gap=5;
TSnow=0;
itmax=51
Yiter=np.zeros(itmax)
effiter=np.zeros(itmax)
Yiter[0]=0


while it<itmax:
    TS=TSnow
    # print(TSnow)
    effnow=min(fsolve(condA,1),1)
    # print(effnow)
    TSnow=effnow*(y+ bet*TS)*lognorm.cdf(effnow*(y+bet*TS), s,0,scale) + quad(Vintegr,effnow*(y+bet*TS),np.inf)[0]- A * effnow**2
    
    gap=TS-TSnow
    
    Yiter[it]=TSnow
    effiter[it-1]=effnow
    it+=1

plt.plot(Yiter)
plt.title('Total surplus as a function of time horizon')
#plt.plot(effiter)

#This was this easy stuff, now things get complicated

Text(0.5, 1.0, 'Total surplus as a function of time horizon')

# Now we do value function iteration
# thoughout, we use the fact that m=y is optimal (at least weakly) which gets rid of a variable
# we choose a grid for the state space, initialize value functions, and make grids for policy functions

grid_size=100
Bstate=np.linspace(0, 1.2*TSnow, num=grid_size).squeeze()
Optout=np.zeros(grid_size) # will use this one to record whether the optimization was successful


# this will record policy functions, for each possible state

xnow=np.zeros((4,grid_size))
# Throughout we will input the contract arguments in the following order:
# x[0] is e
# x[1] is c
# x[2] is By
# x[3] is B0


# now we reduce the time it takes to compute the conditional value of the exit option
# by computing it for a lot of lower bounds
# and then using interpolation below
# this greatly reduces the number of integration steps needed 
# while scipy is optimizing

# we assume that the bound can be no higher optimal surplus (to be checked for robustness later)
lowerbound=np.linspace(0, 3.2*TSnow, num=2*grid_size).squeeze()

#the integrand
def Vintegr(x):
    return x*lognorm.pdf(x,s,0, scale) 

condvalue=[quad(Vintegr,lb,np.inf)[0] for lb in lowerbound] # gross, conditional values of outside options for all possible lb

# now we interpolate once and for all

cv_int=interp1d(lowerbound,condvalue,fill_value='extrapolate')

tol=0.005 # to help scipy a bit, so it doesn't overoptimize, we add mini-slack to our equality constraints
#tol is the precision we are tolerating in equality constraints, to help things optimize


Viter=1; # counts the value function iterate
Ynow=0*TSnow*np.ones(grid_size) #one initialization among many possible ones
Ynext=np.zeros(grid_size)
YnextNC=np.zeros(grid_size) # this here is to make sure that the monotonicty correction does not cause issues

Ynow_int=interp1d(Bstate,Ynow,fill_value='extrapolate')

while(Viter<30):
    
    Biter=0

    while Biter<grid_size:
    
        Bnow=Bstate[Biter]


        # Now we write a function for the RHS of the Bellman equation, this is the objective to maximize
        # Note we will input the contract arguments in the following order:
        # x[0] is e
        # x[1] is c
        # x[2] is By
        # x[3] is B0


        def bellman(x):
            e=x[0]
            c=x[1]
            By=x[2]
            B0=x[3]
    
            Yp=Ynow_int(By)
    
            # prob of staying with the project
            Vbar=e*bet*By+(1-e)*bet*B0+c*R
            pnow=lognorm.cdf(Vbar, s,0,scale)
    
            #print(pnow)
            #print(Vbar)
    
            return -(pnow*e*(y+bet*Yp) +cv_int(Vbar)  - A*e**2) #note that this is minus the value function since we are minimizing
    
        def constraints(x):
            e=x[0]
            c=x[1]
            By=x[2]
            B0=x[3]
    
            Yp=Ynow_int(By)
    
            # prob of staying with the project # redudant calculation here, but safer to call it again
            Vbar=e*bet*By+(1-e)*bet*B0+c*R
            pnow=lognorm.cdf(Vbar, s,0,scale)
    
            h=np.zeros(5) # this will have our 5 "inequality" constraints
    
            #foc for effort first
            foc=pnow*(bet*(By-B0))-2*A*e
    
            if e==1 and foc>=0:  # handles the boundary
                h[0]=0
                h[1]=0 # so optimize know that we are good in this case
            else:
                h[0]=foc+tol
                h[1]=-foc+tol #foc must be approximately zero if e is interior
        
            h[2]=pnow*bet*(e*By+(1-e)*B0)-(1-pnow)*c*R+max(cv_int(Vbar),0)-Bnow+tol-A*e**2 #due to interpolation out of bounds cv_int can go negative, we cap that at zero
            h[3]=-h[2]+2*tol
    
            h[4]=By-B0 # no tolerance needed here, this is a feature the solution must have, no need to look elsewhere
    
            # print('x')
            # print(x)
            # print('h')
            # print(h)
    
            return h

        con={'type':'ineq','fun':constraints}
    
        bnds=((0,1),(0,25*TSnow),(0,2*TSnow),(0,TSnow))    

        x0=[0.5,TSnow/5,TSnow/5,TSnow/5]

        # results=minimize(bellman,x0,method='SLSQP',bounds=bnds,constraints=con,options={'disp': True})
        results=minimize(bellman,x0,method='SLSQP',bounds=bnds,constraints=con)
        
        xnow[:,Biter]=results.x
        Ynext[Biter]=-bellman(results.x)
        YnextNC[Biter]=Ynext[Biter]
        
    
        # this next condition handles the fact that some B targets may be out of bounds in early iterations while Y is low
        # due to short horizon
        # this can cause spurious non-monotonicities early
        if Biter>0 and Ynext[Biter]<Ynext[Biter-1]:
            Ynext[Biter]=Ynext[Biter-1]
        # below we make sure that this step does not cause a spurious convergence problem    
    
        Optout[Biter]=results.status
    
        Biter+=1
    
    # Now we update the value function
    
    # Ynow[0]=Ynow[1]
    # xnow[:,0]=xnow[:,1] # we no longer do this
    
    error=np.max(YnextNC-Ynow) #note that the convergence error is measured without the monotonicty correction
    Ynow=Ynext
    Ynow_int=interp1d(Bstate,Ynow,fill_value='extrapolate')
    
    # and we send it back up
    Viter+=1

# now we collect and plot the results of our efforts
# first value functions

Lnow=Ynext-Bstate # this the lenders value

l=np.linspace(0,Bstate[grid_size-1],num=10) # things we need below to annotate chart

tsx=[-0.5,TSnow]
tsy=[TSnow,TSnow]
tsx2=[TSnow,TSnow]
tsy2=[0,TSnow]
fig, axs  = plt.subplots(1,2,figsize=(12, 8)) 

axs[0].plot(Bstate,Ynext)
axs[0].set_title("Total surplus")
axs[0].set_xlabel(r'Borrower promise ($B$)')
axs[0].set_ylabel(r'$Y$')
# axs[0].plot(l,l,'r')
axs[0].plot(tsx,tsy,'g:')
axs[0].plot(tsx2,tsy2,'g:')
axs[0].set_xlim([0, Bstate[grid_size-1]])
axs[0].set_ylim([0, Bstate[grid_size-5]])

axs[0].annotate(r'$Y^*$',
            xy=(0, TSnow), xycoords='data',
            xytext=(2, 10.5), textcoords='data',
            size=14, va="center", ha="center",
            arrowprops=dict(arrowstyle='->',
                            connectionstyle="arc3,rad=-0.2"),
            )

axs[1].plot(Bstate,Ynext-Bstate)
axs[1].set_title(r'Lender surplus ($L=Y-B$)')
axs[1].set_xlabel(r'Borrower promise ($B$)')
axs[1].set_ylabel(r'$L$')
axs[1].set_xlim([0, Bstate[grid_size-1]])
axs[1].set_ylim([-2, 3.8])
liney=[-2,Bstate[grid_size-1]]
linex=[Bstate[np.argmax(Ynext-Bstate)],Bstate[np.argmax(Ynext-Bstate)]]
axs[1].axvline(x=Bstate[np.argmax(Ynext-Bstate)],color='r',linestyle=':')
axs[1].text(2.5, 2, 'Renegotation proof', fontsize=8)
axs[1].annotate('', xy=(5,1.8), xytext=(3,1.8), arrowprops=dict(arrowstyle='->'))
axs[1].axvspan(0, Bstate[np.argmax(Ynext-Bstate)], alpha=0.05,color='r')
axs[1].axvspan(Bstate[np.argmax(Ynext-Bstate)],Bstate[grid_size-1],alpha=0.05,color='g')

<matplotlib.patches.Polygon at 0x1dedae258b0>

# now we plot policy functions

fig, axsn  = plt.subplots(2,2,figsize=(12, 8)) 

axsn[0][0].plot(Bstate,xnow[0])
axsn[0][0].set_title("Effort")
# axsn[0][0].set_xlabel(r'Borrower promise ($B$)')
# axs[0].set_ylabel(r'$Y$')
# axs[0].plot(l,l,'r')

axsn[0][1].plot(Bstate,xnow[2],label=r'$B_y$')
axsn[0][1].set_title(r'Continuation equity ($B_y$)')
# axsn[0][1].set_xlabel(r'Borrower promise ($B$)')
axsn[0][1].plot(Bstate,Bstate,'r:',label=r'$45^{\circ}$ line')
# axs[0].set_ylabel(r'$Y$')
# axs[0].plot(l,l,'r')
axsn[0][1].legend()

axsn[1][0].plot(Bstate[1:],xnow[1,1:],label=r'$c$')
axsn[1][0].set_title(r'Prepayment penalty')
axsn[1][0].set_xlabel(r'Borrower promise ($B$)')


# now we need to calculate the MV value of loans
# for that we first need to interpolate Lnow

Lnow_int=interp1d(Bstate,Lnow,fill_value='extrapolate')


LMV= 1/R*xnow[0]*[y+bet*Lnow_int(i) for i in xnow[2]] #size of obligation if continuation
#axs[1].plot(Bstate,Bstate,'r:')
# axs[0].set_ylabel(r'$Y$')
# axs[0].plot(l,l,'r')

axsn[1][0].plot(Bstate[1:],LMV[1:],label=r'$\frac{e(m+\beta L_y) + (1-e) \beta L_0}{R}$')
axsn[1][0].legend()
axsn[1][0].fill_between(Bstate[1:],LMV[1:],xnow[1,1:],alpha=0.05,color='red')

taunowr=(xnow[1,1:]*R-R*LMV[1:])/(R*LMV[1:])


axsn[1][1].set_title(r'Relative moneyness $\left( \frac{\tau}{e(m+\beta L_y) + (1-e) \beta L_0}\right)$')
axsn[1][1].set_xlabel(r'Borrower promise ($B$)')
axsn[1][1].plot(Bstate[1:],-taunowr)


#now we try to draw the path of By starting for lender optimal B0

#for that we're going to need, once again, some interpolation


By_int=interp1d(Bstate,xnow[2],fill_value='extrapolate')

it=5
pathB=np.zeros(it)
# pathB[0]=Bstate[np.argmax(Ynext-Bstate)]
pathB[0]=Bstate[2]
i=1
while i<it:
    pathB[i]=By_int(pathB[i-1])
    i+=1
    
# axsn[0][1].arrow(pathB[0],pathB[0],0,pathB[1]-pathB[0])
axsn[0][1].annotate("", xytext=(pathB[0],pathB[0]), xy=(pathB[0], pathB[1]),arrowprops=dict(arrowstyle="->",color='green',linestyle=':'))  
axsn[0][1].annotate("", xytext=(pathB[0],pathB[1]), xy=(pathB[1], pathB[1]),arrowprops=dict(arrowstyle="->",color='green',linestyle=':')) 
axsn[0][1].annotate("", xytext=(pathB[1],pathB[1]), xy=(pathB[1], pathB[2]),arrowprops=dict(arrowstyle="->",color='green',linestyle=':')) 
axsn[0][1].annotate("", xytext=(pathB[1],pathB[2]), xy=(pathB[2], pathB[2]),arrowprops=dict(arrowstyle="->",color='green',linestyle=':')) 
axsn[0][1].annotate("", xytext=(pathB[2],pathB[2]), xy=(pathB[2], pathB[3]),arrowprops=dict(arrowstyle="->",color='green',linestyle=':')) 
axsn[0][1].annotate("", xytext=(pathB[2],pathB[3]), xy=(pathB[3], pathB[3]),arrowprops=dict(arrowstyle="->",color='green',linestyle=':')) 
axsn[0][1].annotate("", xytext=(pathB[3],pathB[3]), xy=(pathB[3], pathB[4]),arrowprops=dict(arrowstyle="->",color='green',linestyle=':')) 
axsn[0][1].annotate("", xytext=(pathB[3],pathB[4]), xy=(12,pathB[4]),arrowprops=dict(arrowstyle="->",color='green',linestyle=':')) 

Text(12.110342885323842, 15.998771006361995, '')