Replication: Dynamic Contracting in Continuous Time

Author
Affiliation

Marcelo Ortiz M.

Universitat Pompeu Fabra

Published

February 7, 2025

State Variable

The cumulative output of the project evolves according to: dX_t = A_t dt + \sigma dB_t

The continuation value W_t is the total utility that the principal expects the agent to derive from the future after a given moment of time t:

W_t = \mathbb{E}^A \left[ \int_t^\infty e^{-r(s-t)} \left( u(C_s) - h(A_s) \right) ds \mid \mathcal{F}_t \right],

where: - r is the discount rate, - u(C_s) is the utility from consumption C_s, - h(A_s) is the disutility from effort A_s, - \mathbb{E}^A denotes the expectation under the probability measure induced by the agent’s effort strategy A, - \mathcal{F}_t represents the information available up to time t.

Define the agent’s total expected utility up to time t as:

V_t = r \int_0^t e^{-rs} \left( u(C_s) - h(A_s) \right) ds + e^{-rt} W_t.

This formulation accounts for both the accumulated utility up to time t and the discounted continuation value W_t.

Under the agent’s optimal effort strategy, V_t should be a martingale with respect to the filtration \mathcal{F}_t generated by the observable output process X_t. This implies:

\begin{aligned} dV_t &= r e^{-rt} \left( u(C_t) - h(A_t) \right) dt - r e^{-rt} W_t dt + e^{-rt} dW_t \\ &= 0. \end{aligned}

The Martingale Representation Theorem states that any martingale can be expressed as an integral with respect to a Brownian motion. Thus, there exists a progressively measurable process Y_t such that:

dV_t = r e^{-rt} Y_t \sigma dZ_t,

where: - \sigma is the volatility of the output process, - Z_t is a standard Brownian motion under the probability measure corresponding to the agent’s effort.

Derive the Dynamics of W_t:

Equating the two expressions for dV_t and solving for dW_t, we obtain:

\begin{align} dW_t = r \left( W_t - u(C_t(W_t)) + h(A_t(W_t)) \right) dt + r Y_t(W_t) (dX_t-A(W_t)dt). \end{align}

Here, Y_t represents the sensitivity of the agent’s continuation value to changes in the output process X_t. Notice that dX_t-A(W_t)dt=\sigma dZ_t.

If the agent takes the contracted level of effort, then \mathbb{E}[dX_t-A(W_t)]=0, so in that case, the drift of the continuation value is r \left( W_t - u(C_t(W_t)) + h(A_t(W_t)) \right).

Interpretation:

  • Drift: The continuation value grows at the interest rate r, minus the agent’s utility from consumption and the disutility from effort at time t.
  • Incentives: If the agent deviate from the contracted effort level, the actual effort impacts the drift of its continuation value through the second term in dW_t.
    • The agent’s optimal effort is chosen to maximize the expected change in the continuation value \mathbb{E}[dW_t] minus the cost of effort h(A_t): rY(W_t)A_t-r*h(A_t). As effort is compensated, the principal will set the sensitivity Y(W_t) that the minimum level that induces the agent to choose the contracted effort level.
    • when the cost of effort is differentiable, the optimal effort sensitivity in the contract is \gamma(a)=h'(a).

Principal’s Hamilton-Jacobi-Bellman Equation

The principal’s optimal control problem is: F(W_0)=\max_{a_t>0, c_t} \mathbb{E} \left[ \int_0^\infty e^{-r(u-t)} ( r(a_u - c_u) du \right], such that: \begin{align} dW_t &= r(W_t - u(c_t) + h(a_t)) dt + r \gamma(a_t)( dX_t-a_t dt) \\ C_t &\geq 0 \\ a_t &\in [0,\hat{A}] \\ \gamma(a_t) &\text{ is incentive compatible} \end{align}

We can convert this optimal control problem into a dynamic problem by defining the principal’s profit function F(W) as: F(W_t) = \max_{a > 0, c} \mathbb{E}_t \left[ \int_t^{t+s} e^{-r(u-t)} (a_u - c_u)du+ e^{-rs}F(W_{t+s}) \right] By applying the Ito’s Lemma to e^{-rs}F(W_{t+s}), dividing by s and letting s goes to zero, we obtain the Hamilton-Jacobi-Bellman (HJB) equation for the principal’s profit function F(W):

\begin{align} r F(W) = \max_{a > 0, c} \left\{ r(a - c) + F'(W) r(W - u(c) + h(a)) + \frac{1}{2} F''(W) r^2 \gamma(a)^2 \sigma^2 \right\} \end{align}

  • r(a - c) is the instantaneous profit flow of the principal,
  • F'(W) r(W - u(c) + h(a)) accounts for how the agent’s continuation value changes due to effort and consumption,
  • \frac{1}{2} F''(W) r^2 \gamma(a)^2 \sigma^2 captures the impact of volatility and the sensitivity of effort in the contract.

Parameter values and Functional forms

  • r = 0.1 is the discount rate,
  • \sigma = 1.0 is the volatility of output process,
  • u(C) = \sqrt{C} is the agent’s utility function,
  • h(A) = 0.5 A^2 + 0.4 A is the effort cost function.

Profits from retiring the agent

To retire the agent with value u(c), the principal can offers him constant c and allows him to choose zero effort. The profits from this for retiring are: F_0(u(c))=-c As there is no effort, the continuation value is the same as the agent’s utility function W_t=U_t, meaning that: F_0(W)=-W^2

Boundary Conditions & Smooth Pasting

The function F(W) must satisfy the following conditions:

  1. Lower boundary condition (zero profit at zero continuation value): F(0) = 0

  2. Retirement boundary condition (profits from retiring the agent): F(W_{gp})=-W_{gp}^2

  3. Smooth pasting condition (ensuring differentiability at retirement): F'(W_{gp}) = -2 W_{gp}

Solution Method

As solving method usually require transforming second-order ODEs into a first-orde ODE, we rewrite the HJB equation as: F''(W) = \min_{a > 0, c} \frac{F(W) - a + c - F'(W)(W - u(c) + h(a))}{r \gamma (a)^2 \sigma^2 / 2}

In the code below we use shooting methods for solving the model.

Code For Figure 1

Show the code
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
from scipy.integrate import solve_ivp
from scipy.optimize import minimize, fsolve

#TODO alternatives to find a_star that depends on F, Fp, c or W?

# Parameters and lambda function
r = 0.1      # discount rate
sigma = 1.0  # volatility
u = lambda c: np.sqrt(c)  # Agent's utility function
h = lambda a: 0.5*a**2 + 0.4*a  # Effort cost function
gamma = lambda a: a + 0.4  # Effort cost function
c_star = lambda Fp: 0 if Fp >= 0 else minimize(lambda c: c + Fp*u(c), 0.5, bounds=[(0.0, 2.0)], tol=1e-12).x[0]
a_star = lambda F, Fp, c, W: minimize(lambda a: - (gamma(a)*a - h(a)), 0.5, bounds=[(0.0, 1)], tol=1e-12).x[0]

def F_ode(W: float, Y: np.array) -> np.array:
    """
    Y[0] = F(W),  Y[1] = F'(W).
    Computes F''(W) from the rearranged HJB equation:
    
    F''(W) = 2 * [F(W) - a + c - F'(W)*(W - u(c) + h(a))]
             ---------------------------------------------
                    r * sigma^2 * (h'(a))^2
    """
    F, Fp = Y
    c = c_star(Fp)
    a = a_star(F, Fp, c, W)
    gamma_star = gamma(a)
    drift = W - u(c) + h(a)
    F_pp = (F - a + c - Fp*drift) / (0.5*r*(sigma**2)*(gamma_star**2))
    return [Fp, F_pp]

def shooting(Fp_guess, return_sol=False):
    """
    First: Given a candidate for Fprime(0), integrate the ODE frontward
    from W=0 to 1.2 (as in the paper), using the boundary conditions at retirement:
    F(0)=0,  F(W_gp) = -W_gp^2,    F'(W_gp) = -2*W_gp.
    
    Second: for each Fp(0),  find W_gp, by finding the first case where  F'(W_gp) = -2*W_gp.
    
    Third: find the Fp(0) such that its W_gp associated meet the condition F(W_gp) = -W_gp^2
    """
    Fp_guess = float(Fp_guess)
    w_grid = np.linspace(0, 1.2, 200)  
    # First:     
    sol = solve_ivp(F_ode,[w_grid[0],w_grid[-1]],  y0=[0,Fp_guess], method='RK45', t_eval=w_grid) # "Lower Boundary Condition" is in y0[0]=0, meaning F(0)=0.
    # Second: Smooth-Pasting Condition, select the first sol.t (=W) where F'(W_gp) = -2*W_gp
    d=sol.y[1,1:] + 2*sol.t[1:]
    idx = np.argmin(np.abs(d))
    idx = idx+1
    W_gp = sol.t[idx]
    # Third: Retirement (Upper) Boundary Condition: compute a loss function for the condition F_Wgp = -W_gp**2
    F0_val = sol.y[0,idx] + W_gp**2
    # return the solution
    if return_sol:
        return sol, W_gp
    return float(F0_val)

# root finding for Fp using an initial guess based on the original plot (fp=0.5)
root = fsolve(shooting, 0.5, xtol=1e-12)

# solve the ODE with the optimal Fp(0)
sol, W_gp = shooting(root[0], return_sol=True)

# Plot the results until W_gp + additional_y_axis
additional_y_axis = 0.1
F_W = sol.y[0]
Fp_W = sol.y[1]
W_vals = sol.t
c_vals, a_vals, drift_vals = [], [], []
for W, F, Fp in zip(W_vals, F_W, Fp_W):
    if W > W_gp + additional_y_axis:
        break
    c = c_star(Fp)
    a = minimize(
        lambda a: (F - a + c - Fp*(W - u(c) + h(a))) / (r*gamma(a)**2*sigma**2/2),  
        0.5, bounds=[(0.0, 1)]
    ).x[0]  # maybe replace by equation 10 in the paper? notice the paper has a typo in the equation.
    c_vals.append(c)
    a_vals.append(a)
    drift_vals.append(r*(W - u(c) + h(a)))

c_vals, a_vals, drift_vals = map(np.array, (c_vals, a_vals, drift_vals))

F0 = -W_vals**2
# Create a figure with two columns, left is wide, right is narrower
fig = plt.figure(figsize=(8, 6))
gs = gridspec.GridSpec(nrows=1, ncols=2, width_ratios=[2.5, 1])

# --- LEFT PLOT (Profit) ---
ax_left = fig.add_subplot(gs[0])
ax_left.plot(W_vals[W_vals <= W_gp + additional_y_axis], F_W[W_vals <= W_gp + additional_y_axis], label='F(W)')
ax_left.plot(W_vals[W_vals <= W_gp + additional_y_axis], F0[W_vals <= W_gp + additional_y_axis], label='F₀(W)')
ax_left.axvline(W_gp, color='red', ls='--', label='W_gp')
ax_left.set_xlabel('W')
ax_left.set_ylabel('Profit')
ax_left.legend()
ax_left.set_title('Principal’s Profit')

# --- RIGHT COLUMN: three stacked plots ---
gs_right = gridspec.GridSpecFromSubplotSpec(3, 1, subplot_spec=gs[1],
                                            hspace=0.4)
# Effort subplot
ax_eff = fig.add_subplot(gs_right[0])
ax_eff.plot(W_vals[W_vals <= W_gp + additional_y_axis], a_vals)
ax_eff.set_ylabel('Effort a(W)')
ax_eff.axvline(W_gp, color='red', ls='--')
ax_eff.set_ylim([0, 1])  # adjust if you want

# Consumption subplot
ax_cons = fig.add_subplot(gs_right[1])
ax_cons.plot(W_vals[W_vals <= W_gp + additional_y_axis], c_vals)
ax_cons.set_ylabel('Consumption c(W)')
ax_cons.axvline(W_gp, color='red', ls='--')
ax_cons.set_ylim([0, 1])  # adjust as needed

# Drift subplot
ax_dr = fig.add_subplot(gs_right[2])
ax_dr.plot(W_vals[W_vals <= W_gp + additional_y_axis], drift_vals)
ax_dr.set_ylabel('Drift of W')
ax_dr.set_xlabel('W')
ax_dr.axvline(W_gp, color='red', ls='--')

plt.tight_layout()
plt.show()