Exposure simulation / PFE and CVA for multi-callable swaps / Bermudan swaptions… Part 1 of 3

In my previous posts we have seen a Monte-Carlo method to generate market scenarios and calculate the expected exposure, potential future exposure and credit value adjustment for a netting set of plain vanilla swaps. In the next three posts we will add multi-callable swaps (Bermudan swaptions) to the netting set.

Roadmap to multi callable products

In the first part we will see how to use a Monte-Carlo simulation for a single-callable swap (European Swaption) pricing. We don’t worry about the model calibration in this posts. This post is intend to fix the used notation and provide a simple example about basic Monte-Carlo pricing. The techniques presented here will be used and extend in the following two posts.

In the second part we develop a regression based approach (aka backward induction method or American Monte-Carlo or Longstaff-Schwartz method) to calculate the npv of a Bermudan swaption. Further we will calibrate the Gaussian short rate model to fit to a set of market prices of European swaptions.

At this point we will be able to calculate the npv of a multi-callable swap but the Monte-Carlo pricing is not suitable for an exposure simulation. Since a nested pricing simulation in a exposure simulation is very time consuming.

Therefore we will modify the American Monte-Carlo method in the third part and make it useable for our exposure simulation. The approach in the third post will follow the method presented by Drs. Alexandre Antonov, Serguei Issakov and Serguei Mechkov in their research paper ‘Backward Induction for Future Values’ and the methodology presented in the book ‘Modelling, Pricing, and Hedging Counterparty Credit Exposure: A Technical Guide’ by Giovanni Cesari et al.

But for now let’s start with something easy: European swaptions.

European swaption pricing

Since my first post we have been living in a single curve world and our model parameter have been being exogenous. To make things even more easy we have been using a flat yield curve. For now we don’t leave this comfortable world and apply the same setting for this example.

How to create a swaption with QuantLib?

An European payer/receiver swaption with physical delivery is an option that allows the option holder at option expiry to enter a payer/receiver swap. The rate paid/received on the fixed leg equals the strike of the swaption.

Given a plain vanilla swap, one can create an European swaption in the QuantLib with very few lines of code. All we need is the expiry date and the settlement type (cash settled or physical delivery).

def makeSwaption(swap, callDates, settlement):
    if len(callDates) == 1:
        exercise = ql.EuropeanExercise(callDates[0])
    else:
        exercise = ql.BermudanExercise(callDates)
    return ql.Swaption(swap, exercise, settlement)

settlementDate = today + ql.Period("1Y")

swaps = [makeSwap(settlementDate,
                  ql.Period("5Y"),
                  1e6,
                  0.03047096,
                  euribor6m)
        ]

calldates = [euribor6m.fixingDate(settlementDate)]
 
swaptions = [makeSwaption(swap, 
                          calldates, 
                          ql.Settlement.Physical) 
             for swap, fd in swaps]

Monte-Carlo pricing

At option expiry T_E the npv of the swaption is V(T_E) = \max(Swap(T_E), 0.0) with Swap(T_E) donating the value of the underlying swap at expiry.

In the Gaussian short rate model under the T-forward measure the zerobond with maturity in T years is our numeraire N(t). Under the usual conditions and using the T-forward measure we can calculate the npv at time 0 by

V(0) = N(0) E[\frac{V(T_E)}{N(T_E)} | F_0].

In our model the numeraire itself is not deterministic so we have to simulate it too.

The Monte-Carlo pricing will consist of three steps

– generate M paths of the short rate process and
– evaluate the swap npv V_i and calculate the numeraire price N_i at option expiry T_E for each path i=0\dots,M-1
– and finally approximate the expected value by \frac{1}{M} \sum_{i=0}^{M-1} \max(\frac{V_i}{N_i},0.0).

Implementation

Instead of using the QuantLib swap pricer we will do the path pricing in Python. Therefore we need to extract the needed information from the instrument.

We convert all dates into times (in years from today). We use the day count convention Act/365.

mcDC = yts.dayCounter()

def timeFromReferenceFactory(daycounter, ref):
    """
    returns a function, that calculate the time in years
    from a the reference date *ref* to date *dat* 
    with respect to the given DayCountConvention *daycounter*
    
    Parameter:
        dayCounter (ql.DayCounter)
        ref (ql.Date)
        
    Return:
    
        f(np.array(ql.Date)) -> np.array(float)
    """
    def impl(dat):
        return daycounter.yearFraction(ref, dat)
    return np.vectorize(impl)

timeFromReference = timeFromReferenceFactory(mcDC, today)

In the first step we extract all fixed leg cashflows and payment dates to numpy arrays.

That are all information we need to calculate the fixed leg npv on a path. We calculate the discount factors for each payment time and multiply the cashflow array with the discount factors array element-wise. The sum of this result gives us the fixed leg npv.

fixed_leg = swap.leg(0)
n = len(fixed_leg)
fixed_times = np.zeros(n)
fixed_amounts = np.zeros(n)
for i in range(n):
    cf = fixed_leg[i]
    fixed_times[i] = timeFromReference(cf.date())
    fixed_amounts[i] = cf.amount()

For the floating leg npv we extract all payment, accrual period start and end dates. We assume that the index start and end dates coincide with the accruals start and end dates and that all periods are regular. With this information we can estimate all floating cashflows by estimating the index fixing through

fixing(t_s) = (\frac{df(t_s)}{df(t_e)}-1) \frac{1}{dcf_{idx}(t_s,t_e)},

with the discount factor df(t) at time t, the year fraction dcf_{idx} between accrual start time t_s and accrual end time t_e using the index day count convention.

float_leg = swap.leg(1)
n = len(float_leg)
float_times = np.zeros(n)
float_dcf = np.zeros(n)
accrual_start_time = np.zeros(n)
accrual_end_time = np.zeros(n)
nominals = np.zeros(n)
for i in range(n):
    # convert base classiborstart_idx Cashflow to
    # FloatingRateCoupon
    cf = ql.as_floating_rate_coupon(float_leg[i])
    iborIndex = cf.index()
    value_date = cf.referencePeriodStart()
    index_mat = cf.referencePeriodEnd()
    # year fraction
    float_dcf[i] = cf.accrualPeriod()
    # calculate the start and end time
    accrual_start_time[i] = timeFromReference(value_date)
    accrual_end_time[i] = timeFromReference(index_mat)
    # payment time
    float_times[i] = timeFromReference(cf.date())
    # nominals 
    nominals[i] = cf.nominal()

We could extend this about gearings and index spreads, but we set the gearing to be one and the spread to be zero.

To calculate the swap npv we need the discount factors for all future payment times (fixed and floating leg), accrual period start and end dates. We store all times together in one array. To get the discount factors we apply the method zeroBond of the GSR model on this array element-wise.

# Store all times for which we need a discount factor in one array
df_times = np.concatenate([fixed_times,
				       ibor_start_time,
				       ibor_end_time,
				       float_times])
df_times = np.unique(df_times)

# Store indices of fix leg payment times in 
# the df_times array
fix_idx = np.in1d(df_times, fixed_times, True)
fix_idx = fix_idx.nonzero()
# Indices of the floating leg payment times 
# in the df_times array
float_idx = np.in1d(df_times, float_times, True)
float_idx = float_idx.nonzero()
# Indices of the accrual start and end time
# in the df_times array
accrual_start_idx = np.in1d(df_times, ibor_start_time, True)
accrual_start_idx = accrual_start_idx.nonzero()
accrual_end_idx = np.in1d(df_times, ibor_end_time, True)
accrual_end_idx = accrual_end_idx.nonzero()

Our pricing algorithm for the underlying swap is:

# Calculate all discount factors
discount = np.vectorize(lambda T: model.zerobond(T, t, x_t))
dfs = discount(df_times)
# Calculate fixed leg npv
fix_leg_npv = np.sum(fixed_amounts * dfs[fix_idx])
# Estimate the index fixings
index_fixings = (dfs[accrual_start_idx] / dfs[accrual_end_idx] - 1) 
index_fixings /= float_dcf
# Calculate the floating leg npv
float_leg_npv = np.sum(nominals * index_fixings * float_dcf * dfs[float_idx])
npv = float_leg_npv - fix_leg_npv

Our time grid for the simulation consists of two points, today and option expiry.

The path generation is very similar like the one in the previous posts, but this time we not only simulate the underlying process but also the numeraires, and we calculate all needed discount factors on a path.

M = 100000
m = len(time_grid)
x = np.zeros((M, m))
y = np.zeros((M, m))
numeraires = np.zeros((M,m))
dfs = np.zeros((M, m, len(df_times)))

for n in range(0,M):
    numeraires[n, 0] = model.numeraire(0, 0)
    
for n in range(0,M):
    dWs = generator.nextSequence().value()
    for i in range(1, len(time_grid)):
        t0 = time_grid[i-1]
        t1 = time_grid[i]
        e = process.expectation(t0, 
                                x[n,i-1], 
                                dt[i-1])
        std = process.stdDeviation(t0,
                                   x[n,i-1],
                                   dt[i-1])
        x[n,i] = e + dWs[i-1] * std 
        e_0_0 = process.expectation(0,0,t1)
        std_0_0 = process.stdDeviation(0,0,t1)
        y[n,i] = (x[n,i] - e_0_0) / std_0_0
        df = np.vectorize(lambda T : model.zerobond(T, t1, y[n,i]))
        numeraires[n ,i] = model.numeraire(t1, y[n, i])
        dfs[n,i] = df(df_times)                            

Given the matrix of numeraires and discount factors we can calculate the npv on the path very fast using numpy arrays.

index_fixings = dfs[:,-1, accrual_start_idx][:,0,:] / dfs[:, -1, accrual_end_idx][:,0,:] - 1
index_fixings /= float_dcf
floatLeg_npv = np.sum(index_fixings * float_dcf * dfs[:,-1, float_idx][:,0,:] * nominals, 
                     axis = 1) 
fixedLeg_npv = np.sum(fixed_amounts * dfs[:, -1, fix_idx][:,0,:], axis=1)
npv = (floatLeg_npv - fixedLeg_npv)
# Apply payoff function 
npv[npv < 0] = 0
# Deflate NPV
npv = npv / numeraires[:,-1] 
npv = np.mean(npv) * numeraires[0,0]

Some remarks

To extract the information from the swap we use the method leg. This method is not a part of the QuantLib 1.5, but you could clone my QuantLib fork on GitHub (branch: SwigSwapExtension) and build the Swig binding yourself. I also send a pull request to Luigi. Maybe it will be part of the official QuantLib at a later time.

In the real world there are quotes for European swaptions in terms of implied volatility available and one would like use a model that is consistent with the market quotes. This is done by model calibration (choice the model parameter so that the model give the same premium for the quoted swaptions). Of cause one could use the Monte-Carlo pricing to calibrate the model, but this would be very time consuming process. The Gaussian short rate model provide some faster and very convenient routines for that. In the next part we will see how to calibrate the model and use the calibrated model to price Bermudan swaptions.

As usual you can download the notebook on nbviewer or GitHub.

Stay tuned for the next part coming soon…

Advertisements

3 thoughts on “Exposure simulation / PFE and CVA for multi-callable swaps / Bermudan swaptions… Part 1 of 3

  1. Hi,

    Could you please elaborate a little bit more on the three last computations you are doing.

    # Apply payoff function
    npv[npv < 0] = 0
    # Deflate NPV
    npv = npv / numeraires[:,-1]
    npv = np.mean(npv) * numeraires[0,0]

    Why do you need to deflate the NPV? And what does the syntax numeraires[:,-1] do? I am programming something similar in C++, so I would really appreciate a few comments about this.

    Like

    1. Hi,

      under the T-Forward measure the deflated payoff (in units of the numeraire) is a martingal. By dividing the payoff by the numeraire value at expiry I convert the payoff into units of the numeraire. After approximating the expected deflated payoff with the mean i convert it back in time-0 Dollars by muliplying it with the time-0 value of the numeraire.

      numeraires is a 2 dimensional Array (row x columns) and stores for each simulation (rows) the numeraire values on the timegrid (0, expiry) (columns). numeraires[:,-1] Returns for all rows the last column. The multiply and divide operator for numpy arrays are element-wise.

      Like

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s