CDO Tranche Pricing

The par spread $S$ of a CDS tranche is given by $$S = \frac{\text{Protection Leg PV}}{\text{RPV01}},$$

where $$\begin{aligned} \text{RPV01} &= \sum_{i} \Delta_i Z(t_i) \frac{Q(t_{i-1},K_1,K_2) + Q(t_i, K_1, K_2)}{2} \\ \text{Protection Leg PV} &= \int_{0}^{T} Z(t) (-dQ(t,K_1,K_2)) \end{aligned}$$

Here,

  • $K_1, K_2, T$ is the tranche starting level, final level, and maturity, respectively.
  • $Z(t)$ is the discount factor curve.
  • $Q(t,K_1,K_2)$ is the tranche survival curve: $Q(t,K_1,K_2)= 1 - \frac{f(t, K_2) - f(t, K_1)}{K_2 - K_1}$
  • $t_i$ is the time of the $i^{th}$ premium payment. And $\Delta_i$ is the time-fraction of the $i^{th}$ interval, using the Act365 DCF.

And $f(t,K) = E[ \min (L(t),K)]$ is model-dependent.

With the exception for the formula of $Q$, the CDO tranche pricing formula is identical to the CDS pricing formula with $R=0$. So we first define a class for calculating CDS's.


In [2]:
from datetime import timedelta
class CDS(object):
    """A generic CDS class. Can also be used for CDO tranches.
    
    Attributes:
        today (datetime.date): The present date.
        maturity (datetime.date): The maturity of the CDO.
        remaining_payments (list): A list of datetime.date objects indicating 
            the remaining payment dates for the premium leg.
        R (float): The recovery fraction (0 to 1).
        Z (callable function): A function Z(t) that outputs the discount factor
            Z(t) for a given time t (datetime.date object). Input and output are
            both positive floats.
        Q (callable function): A function Q(t) that takes in a single
            datetime.date perameter and returns a float representing the tranche
            survival probability. This function should be well-defined for all 
            dates between `today` and `maturity`.
    """
    def __init__(self, today, maturity, remaining_payments, R, Z, Q):
        self.today = today
        self.maturity = maturity
        self.remaining_payments = remaining_payments
        self.R = R
        self.Z = Z
        self.Q = Q
    def rpv01(self):
        """Returns the value of the tranche premium leg for unit notional.
        """
        days = [self.today] + self.remaining_payments
        nodes = [(day - self.today).days / 365 for day in days]
        qvals = [self.Q(day) for day in days]
        total = 0
        for i in range(1, len(days)):
            delta = nodes[i] - nodes[i - 1]
            total += delta * self.Z(days[i]) * (qvals[i] + qvals[i - 1])
        return total / 2
    def protectionLegPV(self, N = 200):
        """Returns the value of the protection leg for unit notional.
        
        Arguments:
            N (int, optional): The number of nodes for calculating the integral.
        """
        delta = (self.maturity - self.today).days / N
        days = [today + timedelta(days = delta) * n for n in range(N + 1)]
        qvals = [self.Q(day) for day in days]
        values = [Z(days[i]) * (qvals[i - 1] - qvals[i])
                  for i in range(1, len(days))]
        return (1 - self.R) * sum(values)
    def parSpread(self, N = 200):
        """ Returns the par spread.
        
        Arguments:
            N (int, optional): The number of nodes for calculating the 
                protection leg integral.
        """
        return self.protectionLegPV(N) / self.rpv01()

Large Homogenous Portfolio (LHP) Model

The LHP tranche pricing model is given by the following equations: $$\begin{aligned} E[\min (L(T), K)] &= (1-R) \Phi_{2}(C(T), -A(K), - \beta) + K \Phi(A(K)) \\ A(K) &= \frac{1}{\beta} \left( C(T) - \sqrt{1 - \beta^{2}} \Phi^{-1} \left( \frac{K}{1-R} \right) \right) \\ C(T) &= \Phi^{-1} (1 - Q(T)) \end{aligned}$$

Here:

  • $\Phi_{2}$ is the cdf of the bivariate normal distribution. The third perameter being correlation.
  • $\beta^{2}$ the correlation
  • $R$ is the recovery rate.
  • $Q(T)$ is the survival curve

In [13]:
from scipy.stats import norm, mvn #normal and bivariate normal
from numpy       import sqrt

def Q_lhp(t, K1, K2, R, beta, Q):
    """Calculates the Tranche survival curve for the LHP model.
    
    Args:
        T (datetime.date): Should be given in "days" (no hours, minutes, etc.)
        K1 (float): The starting value of the tranche. Its value should be 
            between 0 & 1.
        K2 (float): The final value of the tranche. 
        beta (float): Correlation perameter. Its value should be between 0 and 1.
        R (float): Recovery rate. Usually between 0 and 1.
        Q (callable function): A function Q that takes in a dateime.date and
            outputs a float from 0 to 1. It is assumed that Q(0) = 1 and Q is 
            decreasing. Represents survival curve of each credit.
    """
    if Q(t) == 1:
        return 1 # prevents infinity
    def emin(K):
    # Calculates E[min(L(T), K)] in LHP model
        C = norm.ppf(1 - Q(t))
        A = (1 / beta) * (C - sqrt(1 - beta * beta) * norm.ppf(K / (1 - R)))
        return (1 - R) * mvn.mvndst(upper = [C, -1 * A],
                                    lower = [0,0],
                                    infin = [0,0], # set lower bounds = -infty
                                    correl =  -1 * beta)[1] + K * norm.cdf(A)
    return 1 - (emin(K2) - emin(K1)) / (K2 - K1)

Gaussian Heterogenous Model

The tranche survival curve is given by

$$\begin{aligned} Q_{\text{Gauss}} (t, K_1, K_2) & = 1 - \frac{1}{K_2 - K_1} \int_{- \infty}^{+ \infty} \phi (z) f(z) dz \\ f(z) &= \sigma \phi \left(\frac{\mu - K_1}{\sigma} \right) + (\mu - K_1) \Phi \left( \frac{\mu - K_1}{\sigma} \right) - (\mu - K_2) \Phi \left( \frac{\mu - K_2}{\sigma} \right) - \sigma \phi \left( \frac{\mu - K2}{\sigma} \right) \\ \mu &= \frac{1}{N} \sum_{i=1}^{N} p_{i}(t,z) F_i (1 - R_i) \\ \sigma^{2} &= \frac{1}{N^2} \sum_{i=1}^{N} p_i (t,z) (1 - p_i (t,z)) F_i^2 (1-R_i)^{2} \\ p_{i}(t,z) &= \Phi \left( \frac{C_i (t) - \beta_i z}{\sqrt{1 - \beta_i^2}}\right) \\ C_{i} (t) &= \Phi^{-1} (1 -Q_i (t)) \end{aligned}$$

Here:

  • $K_1, K_2$ are the boundaries of the tranche (they are between 0 and 1).
  • $N$ is the number of credits
  • $F_i$ is the fractional face value of the $i$th credit.
  • $R_i$ is the recovery rate of the $i$th credit.
  • $\beta_{i}$ is the $i$th correlation perameter
  • $Q_i$ is the survival curve of the $i$th credit.

Unlike the LHP model, this is heterogenous in that each credit has distinct values of $ R_i, \beta_i, Q_i$. Moreover, the values of $F_i$ becomes relevant.


In [14]:
from scipy.integrate import quad
def Q_gauss(t, K1, K2, Fs, Rs, betas, Qs):
    """ Calculate the tranche survival probability in the Gaussian heterogenous model.
    
    Arguments:
        t (float): Time. Positive.
        K1 (float): Starting tranche value. Between 0 and 1.
        K2 (float): Ending tranche value. Between 0 and 1.
        Fs (list): List of fractional face values for each credit. Each entry must be
            between 0 and 1.
        Rs (list): List of recovery rates for each credit. Each entry must be between
            0 and 1.
        betas (list): Correlation perameters for each credit. Each entry must be between
            0 and 1.
        Qs (list): Survival curves for each credit. Each entry must be a callable function
            that takes in a datetime.date argument and returns a number from 0 to 1.
    """
    Cs = [norm.ppf(1 - q(t)) for q in Qs]
    N = len(Fs) 
    def f(K,z):
        ps = [norm.cdf((C - beta * z) / sqrt(1 - beta ** 2)) for C, beta in zip(Cs, betas)]
        mu = 1 / N * sum([p * F * (1 - R) for p, F, R in zip(ps, Fs, Rs)])
        sigma_squared = 1 / N / N * sum([p * (1 - p) * F ** 2 * (1 - R) ** 2
                                         for p, F, R in zip(ps, Fs, Rs)])
        sigma = sqrt(sigma_squared)
        return -1 * sigma * norm.pdf((mu - K) / sigma) - (mu - K) * norm.cdf((mu - K) / sigma)
    emin = lambda K: quad(lambda z: norm.pdf(z) * f(K, z), -10, 10)[0]
    return 1 - (emin(K2) - emin(K1)) / (K2 - K1)

Adjusted Binomial Model

The adjusted binomial model is given as

$$\begin{aligned} E[min(L(T),K)] &= \sum_{k=0}^{N} h(k) \min (Lk, K) \\ h(k) &= \int_{- \infty}^{\infty} \phi (z) g(k) dz \\ g(k) &= \alpha f(k) \text{ for } k \neq l,u \\ g(l) &= f(l) + (1 - \alpha)(u - m) \\ g(u) &= f(u) - (1 - \alpha)(l - m) \\ l &= \text{floor} (m) \\ m &= p N \\ u &= l + 1 \\ p &= \frac{1}{N} \sum_{i=1}^N \frac{(1 - R_i)F_i}{L} p_i \\ \alpha &= \frac{v_E N + (u - m)^2 -((l -m)^2 - (u -m)^2)(u -m)} {v_A N + (u - m)^2 -((l -m)^2 - (u -m)^2)(u -m)} \\ v_E &= \frac{1}{N^2} \sum_{i=1}^{N} \left( \frac{(1-R_i)F_i}{L}\right)^2 p_i(1 - p_i) \\ v_A &= \frac{p(1-p)}{N} \\ f(k) &= \frac{N!}{(N - k)! k!} p^k (1 - p)^{N - k} \\ L &= \frac{1}{N} \sum_{i=1}^{N} (1 - R_i) F_i \end{aligned}$$

In [15]:
def Q_adjbinom(t, K1, K2, Fs, Rs, betas, Qs):
    """ Calculates the tranche survival probability under the adjusted
        binomial model.
        
        Arguments:
            t (datetime.date): Time.
            K1 (float): Starting tranche value (0 to 1).
            K2 (float): Final tranche value (0 to 1).
            Fs (list): List of fractional face values (floats) for each credit.
            Rs (list): List of recovery rates (floats) for each credit.
            betas (list): List of correlation perameters
            Qs (list): List of survival probabilities. These are callable functions that
                takes in a single datetime.date argument and returns a float.
        Returns:
            float: The value of the tranche survival curve.
    """
    if Qs[0](t) == 1:
        return 1.0 # initial value -- avoids weird nan return
    N = len(Fs)
    Cs = [norm.ppf(1 - Q(t)) for Q in Qs]
    L = sum([(1 - R) * F for R, F in zip(Rs, Fs)]) / N
    def choose(n, k): # Calculates binomial coeffecient: n choose k.
        if k == 0 or k == n:
            return 1
        return choose(n - 1, k - 1) + choose(n - 1, k)
    def g(k,z):
        ps = [norm.cdf((C - beta * z) / sqrt(1 - beta * beta)) for C, beta in zip(Cs, betas)]
        p_avg = sum([(1 - R) * F / L * p for R, F, p in zip(Rs, Fs, ps)]) / N
        f = lambda k: choose(N, k) * p_avg ** k * (1 - p_avg) ** (N - k)
        vA = p_avg * (1 - p_avg) / N
        vE = 1 / N / N * sum([((1 - R) * F / L) ** 2 * p * (1 - p) for R, F, p in zip(Rs, Fs, ps)])
        m = p_avg * N
        l = int(m)
        u = l + 1
        o = (u - m) ** 2 + ((l - m) ** 2 - (u - m) ** 2) * (u - m)
        alpha = (vE * N + o) / (vA * N + o)
        if k == l:
            return f(l) + (1 - alpha) * (u - m)
        if k == u:
            return f(u) - (1 - alpha) * (l - m)
        return alpha * f(k)
    I = lambda k: quad(lambda z: norm.pdf(z) * g(k, z), -10, 10)[0]
    emin = lambda K: sum([I(k) * min(L * k, K) for k in range(0, N + 1)])
    return 1 - (emin(K2) - emin(K1)) / (K2 - K1)

Exact Value (Recursive, Inhomogenous)


In [16]:
def gcd(a, b):
    if a * b == 0:
        return max(a, b)
    if a < b:
        return gcd(a, b % a)
    return gcd(a % b, b)
def Q_exact(t, K1, K2, Fs, Rs, betas, Qs):
    Cs = [norm.ppf(1 - Q(t)) for Q in Qs]
    N = len(Fs)
    g = round(3600 * Fs[0] * (1 - Rs[0]))
    for j in range(1, N):
        g = gcd(g, round(3600 * Fs[j] * (1 - Rs[j])))
    g = g / 3600
    ns = [round(F * (1 - R) / g) for F, R in zip(Fs, Rs)]
    def f(j, k, z):
        if (j, k) == (0, 0):
            return 1.0
        if j == 0:
            return 0.0
        ps = [norm.cdf((C - beta * z) / sqrt(1 - beta ** 2)) for C, beta in zip(Cs, betas)]
        if k < ns[j - 1]:
            return f(j - 1, k, z) * (1 - ps[j - 1])
        return f(j - 1, k, z) * (1 - ps[j -1]) + f(j - 1, k - ns[j - 1], z) * ps[j - 1]
    I = [quad(lambda z: norm.pdf(z) * f(N, k, z), -12, 12)[0] for k in range(sum(ns))]
    emin = lambda K: sum([I[k] * min(k * g, K) for k in range(sum(ns))])
    return 1 - (emin(K2) - emin(K1)) / (K2 - K1)

Example (N = 2 credits)


In [17]:
from numpy import exp
from datetime import date, timedelta
N = 2
K1 = 0.03
K2 = 0.07
Fs = [0.5, 0.5]
Rs = [0.4, 0.6]
def expdecay(today, rate):
    return lambda t: exp(-1 * (t - today).days / 365 * rate)
today = date(2012,1, 1)
Qs = [expdecay(today, 0.0100), expdecay(today, 0.0150)]
betas = [0.40, 0.40]
tvalues = [today + timedelta(days = 30) * n for n in range(37)] #3 years
# LHP perameters average the other perameters
R = 0.5
Q = expdecay(today, 0.0150)
beta = 0.4

In [18]:
from time import time
x = time()
lhpcurve = [Q_lhp(t, K1, K2, R, beta, Q) for t in tvalues]
y = time()
print("lhp: {} ms".format(round(1000 * (y - x))))


lhp: 119 ms

In [19]:
x = time()
gaussiancurve = [Q_gauss(t, K1, K2, Fs, Rs, betas, Qs) for t in tvalues]
y=time()
print("gauss: {} ms".format(round(1000 * (y - x))))


gauss: 21154 ms

In [20]:
x= time()
adjustedbinomialcurve = [Q_adjbinom(t, K1, K2, Fs, Rs, betas, Qs) for t in tvalues]
y = time()
print("adjbinom: {} ms".format(round(1000 * (y - x))))


adjbinom: 39434 ms

In [21]:
x = time()
exactcurve = [Q_exact(t, K1, K2, Fs, Rs, betas, Qs) for t in tvalues]
y = time()
print("exact: {} ms".format(round(1000 * (y - x))))


exact: 47635 ms

In [22]:
%matplotlib inline
import matplotlib.pyplot as plt
plt.plot(tvalues, lhpcurve, 'ro', #red
         tvalues, adjustedbinomialcurve, 'bo', #blue
         tvalues, gaussiancurve, 'go', #green
         tvalues, exactcurve, 'ko'
         )
plt.show


Out[22]:
<function matplotlib.pyplot.show>

In [ ]: