Swap Curve construction using Global Optimization

In my previous post, I discussed option Implied Volatility and Binomial Model. In this post, we switch gears and discuss Swap Curve construction. Building a swap curve is so fundamental to interest rate derivatives pricing that it is probably one of most closely guarded proprietary information on the trading desk.

Pricing an interest rate swap involves bootstrapping a blended curve of different instruments based on their maturities and market liquidity. Usually cash deposits, Eurodollar futures are used at the short end and market swap rates are used at the long end. At the long end, however, we have only a subset of market swap rates available and bootstrapping requires all the missing rates to be interpolated from known rates. This makes interpolation methodology a critical part of curve building. Also, since forward rates are a gradient of the discount rates, any misalignment in the latter is magnified in the former.

There is an alternative approach to bootstrapping called global optimization approach, where the short end of the curve is bootstrapped as usual but at the longer end we “guess” the forward rates, compute the par rates of the market swaps and minimize the error between the actual par rates and the computed par rates. We also add a smoothness constraint to the minimization procedure so that overall gradient of the curve is minimized. This approach is illustrated in the excellent book Swaps and Other derivatives 2nd Ed by Richard Flavell

I will use QuantLib to generate swap schedules and to deal with business day conventions. QuantLib can of course generate a fully built swap curve but I will use Scipy’s optimize package for curve building. My objective was to match Spreadsheet 3.9 “Building a blended curve” from the above book. Unfortunately, the spreadsheet does not actually show the equations for Excel’s Solver used for optimization but shows the final result, so that leaves considerable ambiguity in understanding which I hope I will be able to clear.

Note on QuantLib and Python

There are numerous resources online on how to build QuantLib from sources and then build the Python extensions, I would like to point you to the precompiled package for QuantLib-Python maintained by Christoph Gohlke. If you are on windows, you can just install the whl package and get started.

First some common formulae we will be using:

$$Discount Factor : DF_t = \frac{1}{(1 + r_t * d_t)}$$ where $d_t$ is year fraction and $r_t$ is annual rate

$$Forward Rate : F_{t/T} = \frac{[(DF_t/DF_T)- 1]}{(T-t)}$$ where $DF_t$ is disocunt factor to t and $DF_T$ is disocunt factor to T (both from today)

$$Par  Rate  of  Swap: \frac{(1-D_T)}{\sum_{n=1}^n(\Delta_n * D_n)}$$ where $D_T$ is maturity date discount factor, $\Delta_n$ are the time fractions between 2 reset dates and $D_n$ are the various reset date discount factors.

import QuantLib as ql
import numpy as np
from scipy.optimize import least_squares
from scipy.interpolate import interp1d
import math
import timeit
import matplotlib.pyplot as plt

# we are trying to match Flavell's spreadsheet 3.9
today = ql.Date(4, ql.February,2008)
calendar = ql.JointCalendar(ql.UnitedStates(), ql.UnitedKingdom())
ql.Settings.instance().evaluationDate = today
settle_date = calendar.advance(today, 2, ql.Days)
#define a Swap class to keep state of various dates and year fractions

class MySwap(object):
    def __init__(self, floating_dates, fixed_dates, fixed_rate):
        super(MySwap, self).__init__()
        self.fixed_rate = fixed_rate
        self.floating_dates = floating_dates
        self.fixed_dates = fixed_dates
        self.floating_leg_from_dates = floating_dates[:-1]
        self.floating_leg_to_dates = floating_dates[1:]
        self.fixed_leg_from_dates = fixed_dates[:-1]
        self.fixed_leg_to_dates = fixed_dates[1:]
        self.yf_floating_leg = [ql.Thirty360(ql.Thirty360.BondBasis).yearFraction(d1,d2) for d1,d2
                                in zip(self.floating_leg_from_dates,self.floating_leg_to_dates)]
        self.yf_fixed_leg = [ql.Thirty360(ql.Thirty360.BondBasis).yearFraction(d1, d2) for d1, d2
                                in zip(self.fixed_leg_from_dates, self.fixed_leg_to_dates)]
    
    # for sorting 2 swaps by maturity date   
    def __gt__(self, swap2):
        return self.floating_leg_to_dates[-1] > swap2.floating_leg_to_dates[-1]
#class for Eurodollar futures

class MyFuture(object):
    def __init__(self, begin_date, mkt_yield):
        super(MyFuture, self).__init__()
        self.begin_date = begin_date
        self.mkt_yield=  mkt_yield
#class for Cash deposits

class MyDeposit(object):
    def __init__(self, begin_date, end_date, mkt_yield):
        super(MyDeposit, self).__init__()
        self.begin_date = begin_date
        self.end_date = end_date
        self.mkt_yield=  mkt_yield
        self.yf = ql.Thirty360(ql.Thirty360.BondBasis).yearFraction(begin_date, end_date)
        self.disc_factor = 1/(1+(self.yf*mkt_yield))
# A utility function to recursively calculate a series of discount factors from year fractions
# and forward rates

def get_disc_factors(yf, rates):
    df = [0] * len(yf)
    df[0] = 1.0 / (1.0 + (yf[0] * rates[0]))
    for i in range(1, len(yf)):
        df[i] = 1.0 / (1.0 + (yf[i] * rates[i])) * df[i - 1]
    return df
# this is the objective function which will be minimized using least squares method
# we calculate the par swap rate from the "guessed" forward rates 
# and minimize the difference between this par swap rate and the market quoted par swap rate
# the first 4 forward rates have been derived from deposit and Eurodollar futures and
# are known

#@profile
def scipy_func_to_minimize(fwdGuesses):

    x = np.concatenate([known_fwds,fwdGuesses])
    df_float_leg_longest_swap = get_disc_factors(swaps[-1].yf_floating_leg, x)
    
    def par_rate(yf,df):
        up = 1 - df[-1]
        down = np.sum([y*d for y,d in zip(yf,df)])
        return up/down

    def swap_rate(swap):

        df = np.interp([d.serialNumber() for d in swap.fixed_leg_to_dates],
                       [dl.serialNumber() for dl in swaps[-1].floating_leg_to_dates],
                       df_float_leg_longest_swap)
        return par_rate(swap.yf_fixed_leg, df)

    fi = [swap_rate(s) - s.fixed_rate for s in swaps]
    sq = smoothness_func(x)
    fi.append(sq)
    return fi
# This function ensures a smooth forward curve by minimizing the gradient of the 
# forward rates
# This function is added as an additional constraint to the above 
# objective function

def smoothness_func(fwdGuesses):
    sq = 0.0
    for i, f in enumerate(fwdGuesses[:-1]):
        sq += math.pow((fwdGuesses[i] - fwdGuesses[i + 1]), 2)
    return sq
# A utility function to define a vanilla IRS, specifically USD-3M-Libor 
# fixed-float swap

def makeVanilla1Y3MSwap(maturity, fixedRate, index):

    end = calendar.advance(settle_date, maturity)
    fixedLegTenor = ql.Period("1y")
    fixedLegBDC = ql.ModifiedFollowing
    fixedLegDC = ql.Thirty360(ql.Thirty360.BondBasis)
    spread = 0.0
    fixedSchedule = ql.Schedule(settle_date,
                                end,
                                fixedLegTenor,
                                ql.TARGET(),
                                fixedLegBDC,
                                fixedLegBDC,
                                ql.DateGeneration.Forward,
                                False)
    floatSchedule = ql.Schedule(settle_date,
                                end,
                                index.tenor(),
                                index.fixingCalendar(),
                                index.businessDayConvention(),
                                index.businessDayConvention(),
                                ql.DateGeneration.Forward,
                                False)
    swap = ql.VanillaSwap(ql.VanillaSwap.Payer,
                          1.0,
                          fixedSchedule,
                          fixedRate,
                          fixedLegDC,
                          floatSchedule,
                          index,
                          spread,
                          index.dayCounter())
    return swap
# Market data as of 2/4/2008

libor3m = ql.USDLibor(ql.Period(3, ql.Months))
deposits = [(calendar.advance(settle_date, ql.Period(1, ql.Weeks)),0.032175),
            (calendar.advance(settle_date, ql.Period(1, ql.Months)),0.031813),
            (calendar.advance(settle_date, ql.Period(3, ql.Months)),0.03145)]

depos = [MyDeposit(settle_date, t[0], t[1]) for t in deposits]

# ED Futures start dates : Mar-08,Jun-08,Sep-08,Dec-08,Mar-09
fQuotes =   [97.000,97.410,97.520,97.495,97.395]
futures = []

date_holder = settle_date
for fz in fQuotes:
    fut_begin_date = ql.IMM_nextDate(date_holder)
    futures.append(MyFuture(fut_begin_date, (100.-fz)/100.))
    date_holder = fut_begin_date

ql_swaps = [makeVanilla1Y3MSwap(ql.Period(2, ql.Years),2.795e-2,libor3m),
            makeVanilla1Y3MSwap(ql.Period(3, ql.Years),3.035e-2,libor3m),
            makeVanilla1Y3MSwap(ql.Period(4, ql.Years),3.275e-2,libor3m),
            makeVanilla1Y3MSwap(ql.Period(5, ql.Years),3.5050e-2,libor3m),
            makeVanilla1Y3MSwap(ql.Period(6, ql.Years),3.715e-2,libor3m),
            makeVanilla1Y3MSwap(ql.Period(7, ql.Years),3.885e-2,libor3m),
            makeVanilla1Y3MSwap(ql.Period(8, ql.Years),4.025e-2,libor3m),
            makeVanilla1Y3MSwap(ql.Period(9, ql.Years),4.155e-2,libor3m),
            makeVanilla1Y3MSwap(ql.Period(10, ql.Years),4.265e-2,libor3m),
            makeVanilla1Y3MSwap(ql.Period(12, ql.Years),4.435e-2,libor3m),
            makeVanilla1Y3MSwap(ql.Period(15, ql.Years),4.615e-2,libor3m),
            makeVanilla1Y3MSwap(ql.Period(20, ql.Years),4.755e-2,libor3m),
            makeVanilla1Y3MSwap(ql.Period(25, ql.Years),4.805e-2,libor3m),
            makeVanilla1Y3MSwap(ql.Period(30, ql.Years),4.815e-2,libor3m)]

swaps = [MySwap([d for d in x.floatingSchedule()], [d for d in x.fixedSchedule()], x.fixedRate()) for x in ql_swaps]
swaps.sort()
#Starting to build fwd curve, add deposits first

fwd_curve = [(dep[0], dep[1]) for dep in deposits]

# get next 3 forward rates from futures rates by interpolating on the 
# provided rates 

fut_curve_dates = [calendar.advance(settle_date, ql.Period(3*i, ql.Months)) for i in range(1, 5)]
fut_interpolation = interp1d([f.begin_date.serialNumber() for f in futures], [f.mkt_yield for f in futures])
fwd_curve.extend([(fut_curve_dates[idx], fut_interpolation(fut_curve_dates[idx-1].serialNumber())[()]) 
                  for idx, f in enumerate(fut_curve_dates) if idx > 0])
# At this point we know first 4 forward rates 
# and we guess the next 116 forwards to be same as the 4th forward rate
known_fwds = []
for d in fut_curve_dates:
    known_fwds.extend([f[1] for f in fwd_curve if f[0] == d])
fwd_guesses = [known_fwds[-1]] * (len(swaps[-1].floating_leg_to_dates)-len(known_fwds))
# Now we proceed with optimization and get the optimized forwards 
# that provide the arbitrage-free and smooth forward curve

start_time = timeit.default_timer()
x = least_squares(scipy_func_to_minimize, fwd_guesses, verbose=1).x
print('### Scipy optimization took:%s seconds' % str(timeit.default_timer() - start_time))
optimized_fwds = np.concatenate([known_fwds, x])

`gtol` termination condition is satisfied.
Function evaluations: 20, initial cost: 1.7812e-03, final cost 1.5585e-09, first-order optimality 1.04e-08.
### Scipy optimization took:2.5262434 seconds

#Start building a discount curve as a list of tuples (date,disc_factor), 
# first discount factor = 1.0

dates_df = [(settle_date, 1.0)]
for idx, fwd in enumerate(optimized_fwds):
    temp_tuple = (swaps[-1].floating_leg_to_dates[idx], 
                  dates_df[idx][1]/(1 + (optimized_fwds[idx] * swaps[-1].yf_floating_leg[idx])))
    dates_df.append(temp_tuple)
#Start building a zero curve

zc_rates = []
cumulative_floating_yf = np.cumsum(swaps[-1].yf_floating_leg)
cumulative_floating_yf = np.insert(cumulative_floating_yf, 0, 0.)
entire_date_range = swaps[-1].floating_leg_from_dates + [swaps[-1].floating_leg_to_dates[-1]]
zc_rates = [ (dt, -np.log(dc[1])/yf) for idx, (dc, yf, dt) in enumerate(zip(dates_df,cumulative_floating_yf, entire_date_range)) if idx > 0]
zc_temp_tup = (swaps[-1].floating_leg_from_dates[0], zc_rates[0][1])
zc_rates.insert(0, zc_temp_tup)
#We already have forward rates, just put them in a list of tuples

fwd_rates = [(d,fwd) for d, fwd in zip(swaps[-1].floating_leg_from_dates, optimized_fwds)]
# Plot all 3 curves

_ = plt.figure(figsize=(10, 8.3))

_ = plt.subplot(221)
x_axis = cumulative_floating_yf[1:]
y_axis = list(f[1] for f in fwd_rates)
_ = plt.plot(x_axis, y_axis)
plt.grid()
_ = plt.xlabel('Time')
_ = plt.ylabel('Forward Rates')
_ = plt.title('Forward Curve')

_ = plt.subplot(222)
x_axis = cumulative_floating_yf
y_axis = list(d[1] for d in dates_df)
_ = plt.plot(x_axis, y_axis)
plt.grid()
_ = plt.xlabel('Time')
_ = plt.ylabel('Discount Factors')
_ = plt.title('Discount Curve')

_ = plt.subplot(223)
x_axis = cumulative_floating_yf
y_axis = list(z[1] for z in zc_rates)
_ = plt.plot(x_axis, y_axis)
plt.grid()
_ = plt.xlabel('Time')
_ = plt.ylabel('Zero Rates')
_ = plt.title('Zero Curve')
plt.show()

Option Volatility and Binomial Model

In my previous post Options and Volatility Smile , we used Black-Scholes formula to derive Implied Volatility from given option strike and tenor. Most of the options traded on exchanges are American (with a few index options being European) and can be exercised at any time prior to expiration. Whether it is optimal to exercise an option early depends on whether the stock pays dividend or the level of interest rates and is a very complex subject. What I want to focus on is using Binomial model to price an American option. I summarize below Binomial model theory from the excellent book Derivative Markets 3rd Ed. by Robert McDonald

There are many flavors of the Binomial model but they all have following steps in common:

  • Simulate future prices of underlying stock at various points in time until expiry
  • Calculate the payoff of the option at expiry
  • Discount the payoff back to today to get the option price

1. Simulate future prices of the underlying

Assume the price of a stock is $S_0$ and over a small time period $\Delta t$, the price could either be $S_0u$ or $S_0d$ where u is a factor by which price rises and d is a factor by which price drops. The stock is assumed to follow a random walk, also assume that $p$ is the probability of the stock price rising and $(1-p)$ is the probability of it falling. There are many ways to approach the values of $u$,$d$ and $p$ and the various Binomial models differ in the ways these three parameters are calculated.

In Cox-Ross-Rubinstein (CRR) model, $u = \frac{1}{d}$ is assumed. Since we have 3 unknowns, 2 more equations are needed and those come from risk neutral pricing assumption. Over a small $\Delta t$ the expected return of the stock is

$$pu + (1-p)d = e^{r \Delta t}$$

and the expected variance of the returns is

$$pu^2 + (1-p)d^2 – (e^{r \Delta t})^2 = \sigma ^2 \Delta t$$

Solving for $u$, $d$ and $p$, we get

$$p = \frac{e^{r\Delta t} – d}{u-d}$$

$$u = e^{\sigma \sqrt{\Delta t}}$$

$$d = e^{-\sigma\sqrt{\Delta t}}$$

The CRR model generates a following tree as we simulate multi step stock price movements, this a recombining tree centered around $S_0$.


2. Calculating payoffs at expiry

In this step, we calculate the payoff at each node that corresponds to expiry.

For a put option, $payoff = max(K – S_N, 0)$

For a call option, $payoff = max(S_N – K, 0)$

where $N$ is node at expiry with a stock price $S_N$ and $K$ is the strike.

3. Discounting the payoffs

In this step, we discount the payoffs at expiry back to today using backward induction where we start at expiry node and step backwards through time calculating option value at each node of the tree.

For American put, $V_n = max(K – S_n, e^{-r \Delta t} (p V_u + (1-p) V_d))$

For American call, $V_n = max(S_n – K, e^{-r \Delta t} (p V_u + (1-p) V_d))$

$V_n$ is the option value at node n

$S_n$ is the stock price at node n

$r$ is risk free interest rate

$\Delta t$ is time step

$V_u$ is the option value from the upper node at n+1

$V_d$ is the option value from the lower node at n+1

All the variants of Binomial model, including CRR, converge to Black-Scholes in the limit $\Delta t \to 0$ but the rate of convergence is different. The variant of Binomial model that I would like to use is called Leisen-Reimer Model which converges much faster. Please see the original paper for formulas and a C++ implementation at Volopta.com  which I have ported to Python in the next section.

The python code is going to look very similar to Options and Volatility Smile post except that we will swap out Black-Scholes framework with Leisen-Reimer model. We will also use the same option chain data AAPL_BBG_vols.csv

A note on Python code

I usually do not write code like below, I am purposely avoiding using any classes so that the focus remains on the objective which is to understand the technique.

import math
import numpy as np
import matplotlib as mpl
import pandas as pd
import matplotlib.pyplot as plt
from scipy.integrate import quad
from scipy.optimize import brentq
from scipy.interpolate import interp2d
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm

def BinomialLeisenReimer(S0, K, r, T, N, sigma, div):

    is_call = False

    odd_N = N if (N % 2 == 1) else (N + 1)

    d1 = (math.log(S0 / K) + ((r - div) + (sigma ** 2) / 2.) * T) / (sigma * math.sqrt(T))

    d2 = (math.log(S0 / K) + ((r - div) - (sigma ** 2) / 2.) * T) / (sigma * math.sqrt(T))

    dt = T / odd_N
    r_n = math.exp((r - div) * dt)
    method = 2.
    Term1 = math.pow((d1 / (N + 1.0 / 3.0 - (1 - method) * 0.1 / (N + 1))), 2) * (N + 1.0 / 6.0)
    pp = 0.5 + math.copysign(1, d1) * 0.5 * math.sqrt(1 - math.exp(-Term1))

    Term2 = math.pow((d2 / (N + 1.0 / 3.0 - (1 - method) * 0.1 / (N + 1))), 2) * (N + 1.0 / 6.0)
    p = 0.5 + math.copysign(1, d2) * 0.5 * math.sqrt(1 - math.exp(-Term2))

    u = r_n * pp / p
    d = (r_n - p * u) / (1 - p)

    qu = p
    qd = 1 - p

    #Initialize price tree

    STs = [np.array([S0])]
    # Simulate the possible stock prices path
    for i in range(N):
        prev_branches = STs[-1]
        st = np.concatenate((prev_branches * u, [prev_branches[-1] * d]))
        STs.append(st)  # Add nodes at each time step

    #Initialize payoff tree
    payoffs = np.maximum(0, (STs[N]-K) if is_call else (K-STs[N]))

    def check_early_exercise(payoffs, node):
        early_ex_payoff = (STs[node] - K) if is_call else (K - STs[node])
        return np.maximum(payoffs, early_ex_payoff)

    #Begin tree traversal
    for i in reversed(range(N)):
        # The payoffs from NOT exercising the option
        payoffs = (payoffs[:-1] * qu + payoffs[1:] * qd) * df
        # Payoffs from exercising, for American options
        early_ex_payoff = (STs[i] - K) if is_call else (K - STs[i])
        payoffs = np.maximum(payoffs, early_ex_payoff)

    return payoffs[0]

headers = ['Date','Strike','Call_Bid','Call_Ask','Call','Maturity','Put_Bid','Put_Ask','Put']
dtypes = {'Date': 'str', 'Strike': 'float', 'Call_Bid': 'float', 'Call_Ask': 'float',
          'Call':'float','Maturity':'str','Put_Bid':'float','Put_Ask':'float','Put':'float'}
parse_dates = ['Date', 'Maturity']
data = pd.read_csv('AAPL_BBG_vols.csv',skiprows=1,header=None, names=headers, dtype=dtypes, parse_dates=parse_dates)
data['LR_Imp_Vol'] = 0.0
data['TTM'] = 0.0
r = 0.0248 #risk-free rate
S0 = 153.97 # spot price as of 5/11/2017
div = 0.0182
for i in data.index:
    t = data['Date'][i]
    T = data['Maturity'][i]
    K = data['Strike'][i]
    Put = data['Put'][i]
    sigma_guess = 0.2
    time_to_maturity = (T - t).days/365.0
    data.loc[i, 'TTM'] = time_to_maturity
    def func_LR(sigma):
        return BinomialLeisenReimer(S0, K, r, time_to_maturity, ((T - t).days) - 1, sigma, div) - Put

    lr_imp_vol = brentq(func_LR, 0.03, 1.0)
    data.loc[i, 'LR_Imp_Vol'] = lr_imp_vol

ttm = data['TTM'].tolist()
strikes = data['Strike'].tolist()
imp_vol = data['LR_Imp_Vol'].tolist()
f = interp2d(strikes,ttm,imp_vol, kind='cubic')

plot_strikes = np.linspace(data['Strike'].min(), data['Strike'].max(),25)
plot_ttm = np.linspace(0, data['TTM'].max(), 25)
fig = plt.figure()
ax = fig.gca(projection='3d')
X, Y = np.meshgrid(plot_strikes, plot_ttm)
Z = np.array([f(x,y) for xr, yr in zip(X, Y) for x, y in zip(xr,yr) ]).reshape(len(X), len(X[0]))
surf = ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=cm.coolwarm,linewidth=0.1)
ax.set_xlabel('Strikes')
ax.set_ylabel('Time-to-Maturity')
ax.set_zlabel('Implied Volatility')
fig.colorbar(surf, shrink=0.5, aspect=5)