Pricing Margin Call Stock Loans
Numerical Implementation

Instructor: Dr. Le Nhat Tan

In this R session, we attempt to price a margin call stock loan using a numerical procedure.
Note that the pricing formula for this asset depends on the price of its non-recourse counterpart.

1 Non-Recourse Stock Loan

We provide the stock loan parameters:

  1. \(S,\) the current stock price;
  2. \(E,\) the initial loan amount;
  3. \(\eta,\) the loan interest rate;
  4. \(r,\) the risk-free interest rate;
  5. \(\delta,\) the stock dividend rate;
  6. \(\sigma,\) the stock volatility;
  7. \(T,\) the stock loan time to expiry.

We first apply the numerical procedure at \(n=25\) discrete points.

S = 1.7; E = 0.7; eta = 0.1; r = 0.06; delta = 0.03; sigma = 0.4; T = 15; n = 25
Delta_t = T / n
endpoints = seq(0, T, length.out = n + 1)
B = rep(0, n + 1)
V = rep(0, n + 1)

We provide functions being used in the numerical procedure: \[a(\tau)=Ee^{\eta(T-\tau)},d_{\pm}(x,y,z)=\frac{\ln x-\ln z+(r-\delta\pm\sigma^2/2)y}{\sigma\sqrt{y}},\]

a = function(tau) {
    return(E * exp(eta * (T - tau)))
}

d_plus = function(x, y, z) {
    if (y == 0) {
        if (x < z) {return(-Inf)}
        if (x == z) {return(0)}
        return(Inf)
    }
    return((log(x) - log(z) + (r - delta + sigma ^ 2 / 2) * y) /
              (sigma * sqrt(y)))
}

d_minus = function(x, y, z) {
    if (y == 0) {
        if (x < z) {return(-Inf)}
        if (x == z) {return(0)}
        return(Inf)
    }
    return((log(x) - log(z) + (r - delta - sigma ^ 2 / 2) * y) /
              (sigma * sqrt(y)))
}

\[V_E(S,\tau)=Se^{-\delta\tau}\cdot\Phi(d_+(S,\tau,a(0)))-a(0)e^{-r\tau}\cdot\Phi(d_-(S,\tau,a(0))),\]

V_E = function(x, y) {
    return(x * exp(-delta * y) * pnorm(d_plus(x, y, a(0)))
           - a(0) * exp(-r * y) * pnorm(d_minus(x, y, a(0))))
}

\[J_1(x,y,z,w)=x\delta e^{-\delta(y-z)}\Phi(d_+(x,y-z,w))-a(y)(r-\eta)e^{-(r-\eta)(y-z)}\Phi(d_-(x,y-z,w)),\]

J_1 = function(x, y, z, w) {
    return(x * delta * exp(-delta * (y - z)) * pnorm(d_plus(x, y - z, w))
           - a(y) * (r - eta) * exp(-(r - eta) * (y - z)) * pnorm(d_minus(x, y - z, w)))
}

\[J_2(x,y,w)=\left\{\begin{matrix}x\delta-(r-\eta)a(y), & x>w\\0.5x\delta-0.5(r-\eta)a(y), & x=w\\0, & x<w\end{matrix}\right.,\]

J_2 = function(x, y, w) {
    if (x < w) {return(0)}
    result = x * delta - (r - eta) * a(y)
    if (x > w) {return(result)}
    return(result / 2)
}

\[J_{3,i}(x)=a(x_i)-x+V_E(x,x_i)+\frac{T}{2n}\left(J_1(x,x_i,0,B_1)+J_2(x,x_i,x)\right)+\Delta t\sum_{j=2}^{i-1}J_1(x,x_i,x_j,B_j),\]

J_3 = function(x, i) {
    x_i = endpoints[i]
    result = a(x_i) - x + V_E(x, x_i) + Delta_t * (J_1(x, x_i, 0, B[1]) + J_2(x, x_i, x)) / 2
    if (i == 2) {return(result)}
    for (j in 2:(i - 1)) {
        x_j = endpoints[j]
        result = result + Delta_t * J_1(x, x_i, x_j, B[j])
    }
    return(result)
}

\[\frac{\partial}{\partial x}V_E(x,y)=e^{-\delta y}\Phi(d_+(x,y,a(0)))+\frac{e^{-\delta y}}{\sigma\sqrt{y}}\phi(d_+(x,y,a(0)))-\frac{a(0)e^{-ry}}{\sigma x\sqrt{y}}\phi(d_-(x,y,a(0))),\]

V_E_derivative = function(x, y) {
    return(exp(-delta * y) * pnorm(d_plus(x, y, a(0)))
           + exp(-delta * y) * dnorm(d_plus(x, y, a(0))) / (sigma * sqrt(y))
           - a(0) * exp(-r * y) * dnorm(d_minus(x, y, a(0))) / (sigma * x * sqrt(y)))
}

\[\frac{\partial}{\partial x}J_1(x,y,z,w)=\delta e^{-\delta(y-z)}\Phi(d_+(x,y-z,w))+\frac{\delta e^{-\delta(y-z)}}{\sigma\sqrt{y-z}}\phi(d_+(x,y-z,w))-\frac{a(y)(r-\eta)e^{-(r-\eta)(y-z)}}{\sigma x\sqrt{y-z}}\phi(d_-(x,y-z,w)),\]

J_1_derivative = function(x, y, z, w) {
    return(delta * exp(-delta * (y - z)) * pnorm(d_plus(x, y - z, w)) +
           delta * exp(-delta * (y - z)) * dnorm(d_plus(x, y - z, w)) / (sigma * sqrt(y - z)) -
           a(y) * (r - eta) * exp(-(r - eta) * (y - z)) * dnorm(d_minus(x, y - z, w)) / (sigma * x * sqrt(y - z)))
}

\[J_{3,i}'(x)=-1+\frac{\partial}{\partial x}V_E(x,x_i)+\frac{T}{2n}\left(\frac{\partial}{\partial x}J_1(x,x_i,0,B_1)+\frac{\delta}{2}\right)+\Delta t\sum_{j=2}^{i-1}\frac{\partial}{\partial x}J_1(x,x_i,x_j,B_j).\]

J_3_derivative = function(x, i) {
    x_i = endpoints[i]
    result = -1 + V_E_derivative(x, x_i) + Delta_t * (J_1_derivative(x, x_i, 0, B[1]) + delta / 2) / 2
    if (i == 2) {return(result)}
    for (j in 2:(i - 1)) {
        x_j = endpoints[j]
        result = result + Delta_t * J_1_derivative(x, x_i, x_j, B[j])
    }
    return(result)
}

The Newton-Raphson algorithm gives the \(k^{\textrm{th}}\) order approximation \(B_i^{(k)}\) of \(B_i\) as \[B_i^{(k)}=B_i^{(k-1)}-\frac{J_{3,i}\left(B_i^{(k-1)}\right)}{J_{3,i}'\left(B_i^{(k-1)}\right)}\] with the initial guess \(B_i^{(0)}=B_{i-1}.\) By the Trapezoidal Rule, \[V_i=V_E(S,x_i)+\frac{T}{2n}\left(J_1(S,x_i,0,B_1)+J_2(S,x_i,B_i)\right)+\Delta t\sum_{j=2}^{i-1}J_1(S,x_i,x_j,B_j).\]

# function to calculate V_i
V_value = function(i) {
    x_i = endpoints[i]
    if (i == 1) {return(V_E(S, x_i))}
    result = V_E(S, x_i) + Delta_t * (J_1(S, x_i, 0, B[1]) + J_2(S, x_i, B[i]))
    if (i == 2) {return(result)}
    for (j in 2:(i - 1)) {
        x_j = endpoints[j]
        result = result + Delta_t * J_1(S, x_i, x_j, B[j])
    }
    return(result)
}

# optimal exit price & option price at maturity
B[1] = max(a(0), (r - eta) * a(0) / delta)
V[1] = V_value(1)

# Newton-Raphson algorithm
i = 2
while (i <= n + 1) {
    B[i] = B[i - 1]
    continue = TRUE
    while (continue) {
        approx = B[i] - J_3(B[i], i) / J_3_derivative(B[i], i)
        if (abs(approx - B[i]) / approx < 10 ^ (-6)) {continue = FALSE}
        B[i] = approx
    }
    V[i] = V_value(i)
    i = i + 1
}

An alternative is to use the Secant algorithm, with the \(k^{\textrm{th}}\) order approximation \[B_i^{(k)}=\frac{B_i^{(k-2)}\cdot J_{3,i}\left(B_i^{(k-1)}\right)-B_i^{(k-1)}\cdot J_{3,i}\left(B_i^{(k-2)}\right)}{J_{3,i}\left(B_i^{(k-1)}\right)-J_{3,i}\left(B_i^{(k-2)}\right)}\] with initial guesses \(B_i^{(0)}=B_{i-2},B_i^{(1)}=B_{i-1}.\) Here we choose \(B_2^{(0)}=E.\)

Note that the Secant algorithm helps avoiding calculations on the term \(J'_{3,i}(x).\)

# Secant algorithm
i = 2
while (i <= n + 1) {
    if (i == 2) {previous = E} else {previous = B[i - 2]}
    current = B[i - 1]
    continue = TRUE
    while (continue) {
        B[i] = (previous * J_3(current, i) - current * J_3(previous, i)) / (J_3(current, i) - J_3(previous, i))
        if (abs(current - B[i]) / current < 10 ^ (-6)) {continue = FALSE}
        previous = current
        current = B[i]
    }
    V[i] = V_value(i)
    i = i + 1
}

Visualizations are provided below.

n = 100 # increase time step for curve smoothness
Delta_t = T / n
endpoints = seq(0, T, length.out = n + 1)
B = rep(0, n + 1)
V = rep(0, n + 1)
B[1] = max(a(0), (r - eta) * a(0) / delta)
V[1] = V_value(1)
euro = rep(0, n + 1)
payoff = S - a(endpoints)
payoff = (payoff > 0) * payoff
i = 2
while (i <= n + 1) {
    if (i == 2) {previous = E} else {previous = B[i - 2]}
    current = B[i - 1]
    continue = TRUE
    while (continue) {
        B[i] = (previous * J_3(current, i) - current * J_3(previous, i)) / (J_3(current, i) - J_3(previous, i))
        if (abs(current - B[i]) / current < 10 ^ (-6)) {continue = FALSE}
        previous = current
        current = B[i]
    }
    V[i] = V_value(i)
    euro[i] = V_E(S, endpoints[i])
    i = i + 1
}

# Plot
plot(endpoints, B, type = 'l',
     xlab = expression(paste('Time to Maturity (', tau, ')')), ylab = '',
     main = 'Optimal Exit Boundary vs. Accumulated Loan Principal',
     col = 'violet', lwd = 2, ylim = c(0.5, 5), cex.lab = 1.25)
lines(endpoints, a(endpoints), col = 'blue', lwd = 2, lty = 2)
grid()
legend(x = c(9.3, 15), y = c(4.3, 5),
       legend = c('Optimal Exit Boundary', 'Accumulated Loan Principal'),
       col = c('violet', 'blue'), lwd = 2, lty = c(1, 2))

plot(endpoints, euro, type = 'l',
     xlab = expression(paste('Time to Maturity (', tau, ')')), ylab = '',
     main = 'American vs. European Stock Loan',
     col = 'blue', lwd = 2, ylim = c(-0.1, 1.1), lty = 2, cex.lab = 1.25)
lines(endpoints, V, col = 'violet', lwd = 2)
lines(endpoints, payoff, col = '#007F00', lwd = 2, lty = 4)
grid()
legend(x = c(0, 4.9), y = c(0.75, 1),
       legend = c('American Stock Loan', 'European Stock Loan', 'Payoff'),
       col = c('violet', 'blue', '#007F00'), lwd = 2, lty = c(1, 2, 4))

2 Margin Call Stock Loan

The margin call stock loan has one extra parameter: \(\Delta,\) the payback proportion if the stock price drops down below the accumulated loan.

# Parameters
S = 1.7; E = 1; eta = 0.1; r = 0.06; delta = 0.03; sigma = 0.4; T = 5; n = 50
payback = 0.1

If at any time, the stock price falls down to the accumulated loan \(a(\tau),\) the borrower must pays back a percentage \(\Delta\) of the loan and the contract then becomes a non-recourse one with the corresponding loan value \((1-\Delta)a(\tau).\) We first use the numerical procedure in Section 1 to calculate \[R(\tau)=V_{st}(a(\tau),\tau,(1-\Delta)a(\tau))-\Delta\cdot a(\tau).\] Here \(V_{st}(a(\tau),\tau,(1-\Delta)a(\tau))\) is the price of a stock loan at time to maturity \(\tau\) with the stock price \(a(\tau)\) and accumulated loan \((1-\Delta)a(\tau).\)

Delta_t = T / n
endpoints = seq(0, T, length.out = n + 1)
B = rep(0, n + 1)
V_M = rep(0, n + 1)
prices = a(endpoints)
paybacks = payback * prices
E = E * (1 - payback)

a = function(tau) {
    return(E * exp(eta * (T - tau)))
}

B[1] = max(a(0), (r - eta) * a(0) / delta)

V_value = function(i) {
    S = prices[i]
    x_i = endpoints[i]
    if (i == 1) {return(V_E(S, x_i))}
    result = V_E(S, x_i) + Delta_t * (J_1(S, x_i, 0, B[1]) + J_2(S, x_i, B[i]))
    if (i == 2) {return(result)}
    for (j in 2:(i - 1)) {
        x_j = endpoints[j]
        result = result + Delta_t * J_1(S, x_i, x_j, B[j])
    }
    return(result)
}

V_M[1] = V_value(1)

i = 2
while (i <= n + 1) {
    if (i == 2) {previous = E} else {previous = B[i - 2]}
    current = B[i - 1]
    continue = TRUE
    while (continue) {
        B[i] = (previous * J_3(current, i) - current * J_3(previous, i)) / (J_3(current, i) - J_3(previous, i))
        if (abs(current - B[i]) / current < 10 ^ (-6)) {continue = FALSE}
        previous = current
        current = B[i]
    }
    V_M[i] = V_value(i)
    i = i + 1
}

R = V_M - paybacks

We visualise \(\Delta\cdot a(\tau)\) and \(R(\tau)\) at each \(\tau\in[0, T]\) below.

plot(endpoints, paybacks, type = 'l',
     xlab = expression(paste('Time to Maturity (', tau, ')')),
     ylab = '', main = '', cex.lab = 1.25,
     col = 'violet', lwd = 2, ylim = c(-0.01, 0.25))
lines(endpoints, R, col = 'blue', lwd = 2, lty = 2)
legend(x = c(3.5, 5), y = c(0.2, 0.242),
       legend = c('Payback Proportion', 'Rebate'),
       col = c('violet', 'blue'), lwd = 2, lty = c(1, 2))
grid()

We provide a function for interpolating values on the rebate curve.

R_Linear_Inter = function(x) {
    if (x < 0) {return(R[2] + (R[1] - R[2]) * (endpoints[2] - x) / Delta_t)}
    if (x > T) {return(R[n] + (R[n + 1] - R[n]) * (x - endpoints[n]) / Delta_t)}
    i = 0
    while (endpoints[i + 1] < x & i < n) {i = i + 1}
    return(R[i] + (R[i + 1] - R[i]) * (x - endpoints[i]) / Delta_t)
}

plot(endpoints, R,  xlab = expression(paste('Time to Maturity (', tau, ')')),
     ylab = '', main = '', cex.lab = 1.25,
     col = 'violet', lwd = 2, xlim = c(0, 6), ylim = c(-0.05, 0.25))
x0 = c(-0.1212, 0.1212, 1.2323, 2.3434, 3.4545, 4.5656, 5.6767)
points(x0, lapply(x0, FUN = R_Linear_Inter), col = 'blue', pch = 19)
legend(x = c(4.3, 6), y = c(0.2, 0.25),
       legend = c('Rebate', 'Interpolated Value'),
       col = c('violet', 'blue'), lwd = 2)
grid()

We provide functions being used in the numerical procedure: \[M_1(x,y,z)=xe^{-\delta y}\Phi(d_+(x,y,z))-a(0)e^{-ry}\Phi(d_-(x,y,z)),\]

E = 1

B = rep(0, n + 1)

a = function(tau) {
    return(E * exp(eta * (T - tau)))
}

B[1] = max(a(0), (r - eta) * a(0) / delta)

M_1 = function(x, y, z) {
    return(x * exp(-delta * y) * pnorm(d_plus(x, y, z))
           - a(0) * exp(-r * y) * pnorm(d_minus(x, y, z)))
}

\[\gamma=\frac{2r}{\sigma^2},q=\frac{2\delta}{\sigma^2},k=\gamma-q-\frac{2\eta}{\sigma^2}-1,M(x,y,z)=M_1(x,y,z)-\left(\frac{a(y)}{x}\right)^kM_1\left(\frac{a(y)^2}{x},y,z\right),\]

k = 2 * (r - delta - eta) / (sigma ^ 2) - 1

M = function(x, y, z) {
    return(M_1(x, y, z) - M_1(a(y) ^ 2 / x, y, z) * (a(y) / x) ^ k)
}

\[Q_1(x,y,z,w)=x\delta e^{-\delta(y-z)}\Phi(d_+(x,y-z,w))-a(y)(r-\eta)e^{-(r-\eta)(y-z)}\Phi(d_-(x,y-z,w)),\]

Q_1 = function(x, y, z, w) {
    return(x * delta * exp(-delta * (y - z)) * pnorm(d_plus(x, y - z, w))
           - a(y) * (r - eta) * exp(-(r - eta) * (y - z)) * pnorm(d_minus(x, y - z, w)))
}

\[Q_2(x,y,z,w)=Q_1(x,y,z,w)-\left(\frac{a(y)}{x}\right)^kQ_1\left(\frac{a(y)^2}{x},y,z,w\right),\]

Q_2 = function(x, y, z, w) {
    return(Q_1(x, y, z, w) - Q_1(a(y) ^ 2 / x, y, z, w) * (a(y) / x) ^ k)
}

\[Q_3(x,y,w)=\left\{\begin{matrix}x\delta-(r-\eta)a(y), & x>w\\0.5x\delta-0.5(r-\eta)a(y), & x=w\\0, & x<w\end{matrix}\right.\]

Q_3 = function(x, y, w) {
    return(J_2(x, y, w))
}

\[\alpha=-\frac{k}{2},\beta=-\alpha^2-\gamma,K_2(x,y)=\frac{\ln x-\ln a(y)}{\sigma\sqrt{2y}}\]

alpha = -k / 2
beta = -alpha ^ 2 - 2 * r / (sigma ^ 2)
K_2 = function(x, y) {
    return((log(x) - log(a(y))) / (sigma * sqrt(2 * y)))
}

\[K_3(x,y,v)=\left(\frac{x}{a(y)}\right)^{\alpha}\cdot R\left(y-\frac{y\cdot K_2(x,y)^2}{(v+K_2(x,y))^2}\right)\cdot\exp\left(\frac{(\ln x-\ln a(y))^2\beta}{4(v+K_2(x,y))^2}-(v+K_2(x,y))^2+v\right)\]

K_3 = function(x, y, v) {
    result = (x / a(y)) ^ alpha
    result = result * R_Linear_Inter(y - y * K_2(x, y) ^ 2 / (v + K_2(x, y)) ^ 2)
    result = result * exp(beta * (log(x) - log(a(y))) ^ 2 / 4 / (v + K_2(x, y)) ^ 2 - (v + K_2(x, y)) ^ 2 + v)
    return(result)
}

\[K_1(x,y)=\frac{2}{\sqrt{\pi}}\int_0^{\infty}e^{-v}K_3(x,y,v)dv\approx\frac{2}{\sqrt{\pi}}\sum_{i=1}^{n+1}\omega_iK_3(x,y,l_i)\]

GLQ = gaussLaguerre(n + 1)
K_1 = function(x, y) {
    result = 0
    for (i in 1:(n + 1)) {
        result = result + GLQ$w[i] * K_3(x, y, GLQ$x[i])
    }
    return(result * 2 / sqrt(pi))
}

\[Q_{4,i}(x)=a(x_i)-x+M(x,x_i,a(0))+K_1(x,x_i)+\Delta t\sum_{j=2}^{i-1}Q_2(x,x_i,x_j,B_j)+\frac{T}{2n}\left(Q_3(x,x_i,x)+Q_2(x,x_i,0,B_1)\right)\]

Q_4 = function(x, i) {
    x_i = endpoints[i]
    result = a(x_i) - x + M(x, x_i, a(0)) + K_1(x, x_i) + Delta_t * (Q_3(x, x_i, x) + Q_2(x, x_i, 0, B[1])) / 2
    if (i == 2) {return(result)}
    for (j in 2:(i - 1)) {
        x_j = endpoints[j]
        result = result + Delta_t * Q_2(x, x_i, x_j, B[j])
    }
    return(result)
}

The Secant algorithm gives the \(k^{\textrm{th}}\) order approximation \(B_i^{(k)}\) of \(B_i\) as \[B_i^{(k)}=\frac{B_i^{(k-2)}\cdot Q_{4,i}\left(B_i^{(k-1)}\right)-B_i^{(k-1)}\cdot Q_{4,i}\left(B_i^{(k-2)}\right)}{Q_{4,i}\left(B_i^{(k-1)}\right)-Q_{4,i}\left(B_i^{(k-2)}\right)}\] with initial guesses \(B_i^{(0)}=B_{i-2},B_i^{(1)}=B_{i-1}.\) Here we choose \(B_2^{(0)}=E.\) The margin call stock loan value is approximated as \[V_i=M(S,x_i,a(0))+\frac{T}{2n}(K_1(S,x_i)+Q_3(S,x_i,B_i)+Q_2(S,x_i,0,B_1))+\Delta t\sum_{j=2}^{i-1}(K(S,x_i,x_j)+Q_2(S,x_i,x_j,B_j)).\]

# function to calculate V_i
V_M_value = function(i) {
    x_i = endpoints[i]
    if (i == 1) {return(V_E(S, x_i))}
    result = M(S, x_i, a(0)) + Delta_t * (K_1(S, x_i) + Q_3(S, x_i, B[i]) + Q_2(S, x_i, 0, B[1])) / 2
    if (i == 2) {return(result)}
    for (j in 2:(i - 1)) {
        x_j = endpoints[j]
        result = result + Delta_t * (K_1(S, x_i) + Q_2(S, x_i, x_j, B[j]))
    }
    return(result)
}

V_M[1] = V_M_value(1)

# Secant algorithm
i = 2
while (i <= n + 1) {
    if (i == 2) {
        previous = B[i - 1]
        current = previous + 0.001
    } else {
        previous = B[i - 1]
        current = 2 * B[i - 1] - B[i - 2]
    }
    continue = TRUE
    while (continue) {
        B[i] = (previous * Q_4(current, i) - current * Q_4(previous, i)) / (Q_4(current, i) - Q_4(previous, i))
        if (abs(current - B[i]) / current < 10 ^ (-6)) {continue = FALSE}
        previous = current
        current = B[i]
    }
    V_M[i] = V_M_value(i)
    i = i + 1
}

plot(endpoints, B, xlab = 'Time to Maturity (\U1D70F)', ylab = '', type = 'l',
     main = 'Optimal Exit Price', col = 'violet', lwd = 2)

3 Class Structure

non_recourse = function(S = 1.7, E = 1, eta = 0.1, r = 0.06,
                        delta = 0.03, sigma = 0.4, T = 5, n = 25) {
    Delta_t = T / n
    endpoints = seq(0, T, length.out = n + 1)
    B = rep(0, n + 1)
    V = rep(0, n + 1)
    
    a = function(tau) {
        return(E * exp(eta * (T - tau)))
    }
    
    B[1] = max(a(0), (r - eta) * a(0) / delta)
    
    d_plus = function(x, y, z) {
        if (y == 0) {
            if (x < z) {return(-Inf)}
            if (x == z) {return(0)}
            return(Inf)
        }
        return((log(x) - log(z) + (r - delta + sigma ^ 2 / 2) * y) /
               (sigma * sqrt(y)))
    }

    d_minus = function(x, y, z) {
        if (y == 0) {
            if (x < z) {return(-Inf)}
            if (x == z) {return(0)}
            return(Inf)
        }
        return((log(x) - log(z) + (r - delta - sigma ^ 2 / 2) * y) /
               (sigma * sqrt(y)))
    }
    
    V_E = function(x, y) {
        return(x * exp(-delta * y) * pnorm(d_plus(x, y, a(0)))
               - a(0) * exp(-r * y) * pnorm(d_minus(x, y, a(0))))
    }
    
    J_1 = function(x, y, z, w) {
        return(x * delta * exp(-delta * (y - z)) * pnorm(d_plus(x, y - z, w))
               - a(y) * (r - eta) * exp(-(r - eta) * (y - z)) * pnorm(d_minus(x, y - z, w)))
    }
    
    J_2 = function(x, y, w) {
        if (x < w) {return(0)}
        result = x * delta - (r - eta) * a(y)
        if (x > w) {return(result)}
        return(result / 2)
    }
    
    J_3 = function(x, i) {
        x_i = endpoints[i]
        result = a(x_i) - x + V_E(x, x_i) + Delta_t * (J_1(x, x_i, 0, B[1]) + J_2(x, x_i, x)) / 2
        if (i == 2) {return(result)}
        for (j in 2:(i - 1)) {
            x_j = endpoints[j]
            result = result + Delta_t * J_1(x, x_i, x_j, B[j])
        }
        return(result)
    }
    
    V_value = function(i) {
        x_i = endpoints[i]
        if (i == 1) {return(V_E(S, x_i))}
        result = V_E(S, x_i) + Delta_t * (J_1(S, x_i, 0, B[1]) + J_2(S, x_i, B[i]))
        if (i == 2) {return(result)}
        for (j in 2:(i - 1)) {
            x_j = endpoints[j]
            result = result + Delta_t * J_1(S, x_i, x_j, B[j])
        }
        return(result)
    }
    
    V[1] = V_value(1)
    
    i = 2
    while (i <= n + 1) {
        if (i == 2) {previous = E} else {previous = B[i - 2]}
        current = B[i - 1]
        continue = TRUE
        while (continue) {
            B[i] = (previous * J_3(current, i) - current * J_3(previous, i)) / (J_3(current, i) - J_3(previous, i))
            if (abs(current - B[i]) / current < 10 ^ (-6)) {continue = FALSE}
            previous = current
            current = B[i]
        }
        V[i] = V_value(i)
        i = i + 1
    }
    
    return(list(optimal_exit = B, price = V))
}
n = 100
endpoints = seq(0, T, length.out = n + 1)

nr1 = non_recourse(n = n, T = T)
nr2 = non_recourse(n = n, T = T, eta = 0.15)

plot(endpoints, nr1$optimal_exit, type = 'l',
     xlab = expression(paste('Time to Maturity (', tau, ')')), ylab = '',
     main = 'Optimal Exit Boundary', cex.lab = 1.25,
     col = 'violet', lwd = 2, ylim = c(1.6, 2.8))
lines(endpoints, nr2$optimal_exit, col = 'blue', lwd = 2, lty = 2)
grid()
legend(x = c(4, 5), y = c(2.6, 2.8),
       legend = c(expression(paste(eta, ' = 0.1')),
                  expression(paste(eta, ' = 0.15'))),
       col = c('violet', 'blue'), lwd = 2, lty = c(1, 2))

nr2 = non_recourse(n = n, T = T, sigma = 0.2)
nr3 = non_recourse(n = n, T = T, sigma = 0.45)

plot(endpoints, nr1$optimal_exit, type = 'l',
     xlab = expression(paste('Time to Maturity (', tau, ')')), ylab = '',
     main = 'Optimal Exit Boundary', cex.lab = 1.25,
     col = 'violet', lwd = 2, ylim = c(1, 3))
lines(endpoints, nr2$optimal_exit, col = 'blue', lwd = 2, lty = 2)
lines(endpoints, nr3$optimal_exit, col = '#007F00', lwd = 2, lty = 4)
grid()
legend(x = c(0, 0.9), y = c(1, 1.4),
       legend = c(expression(paste(sigma, ' = 0.2')),
                  expression(paste(sigma, ' = 0.4')),
                  expression(paste(sigma, ' = 0.45'))),
       col = c('blue', 'violet', '#007F00'), lwd = 2, lty = c(2, 1, 4))

nr2 = non_recourse(n = n, T = T, r = 0.08)

plot(endpoints, nr1$optimal_exit, type = 'l',
     xlab = expression(paste('Time to Maturity (', tau, ')')), ylab = '',
     main = 'Optimal Exit Boundary', cex.lab = 1.25,
     col = 'violet', lwd = 2, ylim = c(1.6, 2.8))
lines(endpoints, nr2$optimal_exit, col = 'blue', lwd = 2, lty = 2)
grid()
legend(x = c(4.1, 5), y = c(2.6, 2.8),
       legend = c('r = 0.06', 'r = 0.08'),
       col = c('violet', 'blue'), lwd = 2, lty = c(1, 2))

nr2 = non_recourse(n = n, T = T, delta = 0.02)

plot(endpoints, nr1$optimal_exit, type = 'l',
     xlab = expression(paste('Time to Maturity (', tau, ')')), ylab = '',
     main = 'Optimal Exit Boundary', cex.lab = 1.25,
     col = 'violet', lwd = 2, ylim = c(1.6, 2.8))
lines(endpoints, nr2$optimal_exit, col = 'blue', lwd = 2, lty = 2)
grid()
legend(x = c(4.1, 5), y = c(2.6, 2.8),
       legend = c(expression(paste(delta, ' = 0.02')),
                  expression(paste(delta, ' = 0.03'))),
       col = c('blue', 'violet'), lwd = 2, lty = c(2, 1))

margin_call = function(S = 1.7, E = 1, eta = 0.1, r = 0.06, delta = 0.03,
                       sigma = 0.4, T = 5, n = 37, payback = 0.1) {
    Delta_t = T / n
    endpoints = seq(0, T, length.out = n + 1)
    B = rep(0, n + 1)
    V = rep(0, n + 1)
    
    a = function(tau) {
        return(E * exp(eta * (T - tau)))
    }
    
    prices = a(endpoints)
    paybacks = payback * prices
    E1 = E * (1 - payback)

    a = function(tau) {
        return(E1 * exp(eta * (T - tau)))
    }

    B[1] = max(a(0), (r - eta) * a(0) / delta)
    
    d_plus = function(x, y, z) {
        if (y == 0) {
            if (x < z) {return(-Inf)}
            if (x == z) {return(0)}
            return(Inf)
        }
        return((log(x) - log(z) + (r - delta + sigma ^ 2 / 2) * y) /
               (sigma * sqrt(y)))
    }

    d_minus = function(x, y, z) {
        if (y == 0) {
            if (x < z) {return(-Inf)}
            if (x == z) {return(0)}
            return(Inf)
        }
        return((log(x) - log(z) + (r - delta - sigma ^ 2 / 2) * y) /
               (sigma * sqrt(y)))
    }
    
    V_E = function(x, y) {
        return(x * exp(-delta * y) * pnorm(d_plus(x, y, a(0)))
               - a(0) * exp(-r * y) * pnorm(d_minus(x, y, a(0))))
    }
    
    J_1 = function(x, y, z, w) {
        return(x * delta * exp(-delta * (y - z)) * pnorm(d_plus(x, y - z, w))
               - a(y) * (r - eta) * exp(-(r - eta) * (y - z)) * pnorm(d_minus(x, y - z, w)))
    }
    
    J_2 = function(x, y, w) {
        if (x < w) {return(0)}
        result = x * delta - (r - eta) * a(y)
        if (x > w) {return(result)}
        return(result / 2)
    }
    
    J_3 = function(x, i) {
        x_i = endpoints[i]
        result = a(x_i) - x + V_E(x, x_i) + Delta_t * (J_1(x, x_i, 0, B[1]) + J_2(x, x_i, x)) / 2
        if (i == 2) {return(result)}
        for (j in 2:(i - 1)) {
            x_j = endpoints[j]
            result = result + Delta_t * J_1(x, x_i, x_j, B[j])
        }
        return(result)
    }

    V_value = function(i) {
        S = prices[i]
        x_i = endpoints[i]
        if (i == 1) {return(V_E(S, x_i))}
        result = V_E(S, x_i) + Delta_t * (J_1(S, x_i, 0, B[1]) + J_2(S, x_i, B[i]))
        if (i == 2) {return(result)}
        for (j in 2:(i - 1)) {
            x_j = endpoints[j]
            result = result + Delta_t * J_1(S, x_i, x_j, B[j])
        }
        return(result)
    }

    V[1] = V_value(1)

    i = 2
    while (i <= n + 1) {
        if (i == 2) {previous = E1} else {previous = B[i - 2]}
        current = B[i - 1]
        continue = TRUE
        while (continue) {
            B[i] = (previous * J_3(current, i) - current * J_3(previous, i)) / (J_3(current, i) - J_3(previous, i))
            if (abs(current - B[i]) / current < 10 ^ (-6)) {continue = FALSE}
            previous = current
            current = B[i]
        }
        V[i] = V_value(i)
        i = i + 1
    }

    R = V - paybacks
    
    R_Linear_Inter = function(x) {
        if (x < 0) {return(R[2] + (R[1] - R[2]) * (endpoints[2] - x) / Delta_t)}
        if (x > T) {return(R[n] + (R[n + 1] - R[n]) * (x - endpoints[n]) / Delta_t)}
        i = 0
        while (endpoints[i + 1] < x & i < n) {i = i + 1}
        return(R[i] + (R[i + 1] - R[i]) * (x - endpoints[i]) / Delta_t)
    }
    
    a = function(tau) {
        return(E * exp(eta * (T - tau)))
    }
    
    B[1] = max(a(0), (r - eta) * a(0) / delta)

    M_1 = function(x, y, z) {
        return(x * exp(-delta * y) * pnorm(d_plus(x, y, z))
               - a(0) * exp(-r * y) * pnorm(d_minus(x, y, z)))
    }
    
    k = 2 * (r - delta - eta) / (sigma ^ 2) - 1

    M = function(x, y, z) {
        return(M_1(x, y, z) - M_1(a(y) ^ 2 / x, y, z) * (a(y) / x) ^ k)
    }
    
    Q_1 = function(x, y, z, w) {
        return(x * delta * exp(-delta * (y - z)) * pnorm(d_plus(x, y - z, w))
               - a(y) * (r - eta) * exp(-(r - eta) * (y - z)) * pnorm(d_minus(x, y - z, w)))
    }
    
    Q_2 = function(x, y, z, w) {
        return(Q_1(x, y, z, w) - Q_1(a(y) ^ 2 / x, y, z, w) * (a(y) / x) ^ k)
    }
    
    Q_3 = function(x, y, w) {
        return(J_2(x, y, w))
    }
    
    alpha = -k / 2
    beta = -alpha ^ 2 - 2 * r / (sigma ^ 2)
    
    K_2 = function(x, y) {
        return((log(x) - log(a(y))) / (sigma * sqrt(2 * y)))
    }
    
    K_3 = function(x, y, v) {
        result = (x / a(y)) ^ alpha
        result = result * R_Linear_Inter(y - y * K_2(x, y) ^ 2 / (v + K_2(x, y)) ^ 2)
        result = result * exp(beta * (log(x) - log(a(y))) ^ 2 / 4 / (v + K_2(x, y)) ^ 2 - (v + K_2(x, y)) ^ 2 + v)
        return(result)
    }
    
    GLQ = gaussLaguerre(n + 1)
    
    K_1 = function(x, y) {
        result = 0
        for (i in 1:(n + 1)) {
            result = result + GLQ$w[i] * K_3(x, y, GLQ$x[i])
        }
        return(result * 2 / sqrt(pi))
    }
    
    Q_4 = function(x, i) {
        x_i = endpoints[i]
        result = a(x_i) - x + M(x, x_i, a(0)) + K_1(x, x_i) + Delta_t * (Q_3(x, x_i, x) + Q_2(x, x_i, 0, B[1])) / 2
        if (i == 2) {return(result)}
        for (j in 2:(i - 1)) {
            x_j = endpoints[j]
            result = result + Delta_t * Q_2(x, x_i, x_j, B[j])
        }
        return(result)
    }
        
    V_value = function(i) {
        x_i = endpoints[i]
        if (i == 1) {return(V_E(S, x_i))}
        result = M(S, x_i, a(0)) + Delta_t * (K_1(S, x_i) + Q_3(S, x_i, B[i]) + Q_2(S, x_i, 0, B[1])) / 2
        if (i == 2) {return(result)}
        for (j in 2:(i - 1)) {
            x_j = endpoints[j]
            result = result + Delta_t * (K_1(S, x_i) + Q_2(S, x_i, x_j, B[j]))
        }
        return(result)
    }

    V[1] = V_value(1)
    
    i = 2
    while (i <= n + 1) {
        previous = B[i - 1]
        if (i == 2) {current = previous + 0.001} else {current = 2 * B[i - 1] - B[i - 2]}
        continue = TRUE
        while (continue) {
            B[i] = (previous * Q_4(current, i) - current * Q_4(previous, i)) / (Q_4(current, i) - Q_4(previous, i))
            if (abs(current - B[i]) / current < 10 ^ (-6)) {continue = FALSE}
            previous = current
            current = B[i]
        }
        V[i] = V_value(i)
        i = i + 1
    }
    
    return(list(optimal_exit = B, price = V))
}
n = 80
endpoints = seq(0, T, length.out = n + 1)

mc1 = margin_call(n = n, T = T)
mc2 = margin_call(n = n, T = T, payback = 0)
mc3 = margin_call(n = n, T = T, payback = 0.05)

plot(endpoints, mc1$optimal_exit, type = 'l',
     xlab = expression(paste('Time to Maturity (', tau, ')')), ylab = '',
     main = 'Optimal Exit Boundary', cex.lab = 1.25,
     col = 'violet', lwd = 2, ylim = c(1.6, 2.6))
lines(endpoints, mc2$optimal_exit, col = 'blue', lwd = 2, lty = 2)
lines(endpoints, mc3$optimal_exit, col = '#007F00', lwd = 2, lty = 4)
grid()
legend(x = c(3.7, 5), y = c(2.4, 2.6),
       legend = c('Non-Recourse',
                  expression(paste(Delta, ' = 0.05')),
                  expression(paste(Delta, ' = 0.1'))),
       col = c('blue', 'violet', '#007F00'), lwd = 2, lty = c(2, 1, 4))

n = 40
endpoints = seq(0, T, length.out = n + 1)

mc1 = margin_call(n = n, T = T)
mc2 = margin_call(n = n, T = T, eta = 0.15)

plot(endpoints, mc1$optimal_exit, type = 'l',
     xlab = expression(paste('Time to Maturity (', tau, ')')), ylab = '',
     main = 'Optimal Exit Boundary', cex.lab = 1.25,
     col = 'violet', lwd = 2, ylim = c(1.6, 2.6))
lines(endpoints, mc2$optimal_exit, col = 'blue', lwd = 2, lty = 2)
grid()
legend(x = c(4.05, 5), y = c(2.46, 2.6),
       legend = c(expression(paste(eta, ' = 0.1')),
                  expression(paste(eta, ' = 0.15'))),
       col = c('violet', 'blue'), lwd = 2, lty = c(1, 2))

mc2 = margin_call(n = n, T = T, sigma = 0.2)
mc3 = margin_call(n = n, T = T, sigma = 0.45)

plot(endpoints, mc1$optimal_exit, type = 'l',
     xlab = expression(paste('Time to Maturity (', tau, ')')), ylab = '',
     main = 'Optimal Exit Boundary', cex.lab = 1.25,
     col = 'violet', lwd = 2, ylim = c(1.2, 2.6))
lines(endpoints, mc2$optimal_exit, col = 'blue', lwd = 2, lty = 2)
lines(endpoints, mc3$optimal_exit, col = '#007F00', lwd = 2, lty = 4)
grid()
legend(x = c(0, 0.9), y = c(1.2, 1.5),
       legend = c(expression(paste(sigma, ' = 0.2')),
                  expression(paste(sigma, ' = 0.4')),
                  expression(paste(sigma, ' = 0.45'))),
       col = c('blue', 'violet', '#007F00'), lwd = 2, lty = c(2, 1, 4))

mc2 = margin_call(n = n, T = T, r = 0.08)

plot(endpoints, mc1$optimal_exit, type = 'l',
     xlab = expression(paste('Time to Maturity (', tau, ')')), ylab = '',
     main = 'Optimal Exit Boundary', cex.lab = 1.25,
     col = 'violet', lwd = 2, ylim = c(1.6, 2.5))
lines(endpoints, mc2$optimal_exit, col = 'blue', lwd = 2, lty = 2)
grid()
legend(x = c(4.1, 5), y = c(2.33, 2.47),
       legend = c('r = 0.06', 'r = 0.08'),
       col = c('violet', 'blue'), lwd = 2, lty = c(1, 2))

mc2 = margin_call(n = n, T = T, delta = 0.02)

plot(endpoints, mc1$optimal_exit, type = 'l',
     xlab = expression(paste('Time to Maturity (', tau, ')')), ylab = '',
     main = 'Optimal Exit Boundary', cex.lab = 1.25,
     col = 'violet', lwd = 2, ylim = c(1.6, 2.4))
lines(endpoints, mc2$optimal_exit, col = 'blue', lwd = 2, lty = 2)
grid()
legend(x = c(4.1, 5), y = c(2.27, 2.4),
       legend = c(expression(paste(delta, ' = 0.02')),
                  expression(paste(delta, ' = 0.03'))),
       col = c('blue', 'violet'), lwd = 2, lty = c(2, 1))