American Monte Carlo

In this post we are going to simulate exposures of a bermudan swaption (with physical settlement) through an backward induction (aka American Monte-Carlo) method. These exposure simulations are often needed at the counterparty credit risk management. We will use this exposures to calculate the Credit Value Adjustment (CVA). One could easily extend this notebook to calculate the PFE or other xVAs. The implementation is based on the paper “Backward Induction for Future Values” by Dr. Antonov et al. .

We assume that the value of a derivate and the default probability of the counterpart are independent (no wrong-way risk). For our example our netting set consists of only one Bermudan swaption (but it can be easily extended). For simplicity we assume a flat yield termstructure.

In our example we consider a 5Y Euribor 6M payer swap at 3 percent. The bermudan swaption allows us to enter this swap on each fixed rate payment date.

Just as a short reminder the unilateral CVA of a derivate is given by:

$CVA = (1-R) \int_0^Tdf(t)EE(t)dPD(t),$

with recovery rate R, portfolio maturity T (the latest maturity of all deals in the netting set), discount factor df(t) at time t, the expected exposure EE(t) of the netting set at time t and the default probability PD(t).

We can approximate the integral by the discrete sum

$CVA = (1-R) \sum_{i=1}^n df(t_i)EE(t_i)(PD(t_{i-1})-PD(t_i)).$

In one of my previous posts we calculated the CVA of a plain vanilla swap. We performed the following steps:

Generate a timegrid T
Generate N paths of the underlying market values which influences the value of the derivate
For each point in time calculate the positive expected exposure
Approximate the integral

As in the plain vanilla case we will use a short rate process and assume its already calibrated.

The first two steps are exactly the same as in the plain vanilla case.

But how can we calculate the expected exposure of the swaption?

In the T-forward measure the expected exposure is given by:

$EE(t) = P(0,T) E^T[\frac{1}{P(t,T)} \max(V(t), 0)]$

with future value of the bermudan swaption $V(t)$ at time $t$ . The npv of a physical settled swaption can be negative at time t if the option has been exercised before time t and the underlying swap has a negative npv at time t.

For each simulated path we need to calculate the npv of the swaption $V(t_i,x_i)$ conditional the state $x_i$ at time $t_i$.

In the previous post we saw a method how to calculate the npv of a bermudan swaption. But for the calculation of the npv we used a Monte Carlo Simulation. This would result again in a nested Monte-Carlo-Simulation which is, for obvious reasons, not desirable.

If we have a look on the future value, we see that is very much the same as the continuation value in the bermudan swaption pricing problem.

But instead of calculate the continuation values only one exercises date we calculate it now on each time in our grid. We use a the same regression based approximation.

We approximate the continuation value through some function of current state of the stochastic process x_i:

$V(t_i, x_i) \approx f(x_i)$ .

A common choice for the function is of the form:

$f(x) = a_0 + a_1 g_1(x)+ a_2 g_2(x) + \dots a_n g_n(x)$

with $g_i (x)$ a polynom of the degree $i$ .

The coefficients of this function f are estimated by the ordinary least square error method.

We can almost reuse the code from the previous post, only a few extensions are needed.

We need to add more times to our time grid (we add a few more points in time after the maturity to have some nicer plots):

date_grid = [today + ql.Period(i, ql.Months) for i in range(0,66)] + calldates + fixing_dates

When we exercise the option we will enter the underlying swap. Therefore we need to update the swaption npvs to the underlying swap npv for all points in time after the exercise time on each grid:

if t in callTimes:
cont_value = np.maximum(cont_value_hat, exercise_values)
swaption_npvs[cont_value_hat < exercise_values, i:] = swap_npvs[cont_value_hat < exercise_values, i:].copy()

For our example we can observe the following simulated exposures:

If we compare it with the exposures of the underlying swap, we can see that the some exposure paths coincide. This is the case when the swaption has been exercised. In some of this cases we observe also negative npvs after exercising the option.

For the CVA calculation we need the positive expected exposure, which can be easily calculated:

swap_npvs[swap_npvs<0] = 0
swaption_npvs[swaption_npvs<0] = 0
EE_swaption = np.mean(swaption_npvs, axis=0)
EE_swap = np.mean(swap_npvs, axis=0)

Given the default probability of the counterparty as a default termstructure we can now calculate the CVA for the bermudan swaption:

# Calculation of the default probs
defaultProb_vec = np.vectorize(pd_curve.defaultProbability)
dPD = defaultProb_vec(time_grid[:-1], time_grid[1:])

# Calculation of the CVA
recovery = 0.4
CVA = (1-recovery) * np.sum(EE_swaption[1:] * dPD)

As usual the complete notebook is available here.

Hey Everyone,

today I want to continue my last post and show you today how to calculate the NPV of a bermudan swaption through Monte-Carlo simulation. The techniques I will show you in this post can be easily extended to simulate exposure paths of a multi callable swap. This exposure simulations can be used for a xVA calculation like the Credit Value Adjustment (CVA) or Potential Future Exposure (PFE).

A physical settled Bermudan swaption is a path dependent derivate. The option has multiple exercise dates and the swaption holder has the right to enter the underlying swap at any of this exercise dates.

On each of the exercise dates the holder have to decide whether it is optimal to exercise the option now or continue to hold the option and may exercise it on a later date.

Lets consider a bermudan swaption with two exercise dates.

At the latest exercise date T the payoff is the well known european swaption payoff

$V(T) = \max(S(T), 0)$ ,

where $S(T)$ is the npv of the underlying swap at time t.

At the first exercise date $t_i \le T$ the npv of the swaption given by

$V(t_i) = \max(S(t_i), N(t_i) E [ \frac{V(T)} {N(T)} | F_{t_i}]$ .

and so to the NPV at time 0

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

One way to solve problem is performing a Monte-Carlo-Simulation. But a naive Monte Carlo approach would require a nested Monte-Carlo Simulation on each path to calculate the continuation value
$C(t_i ) = N(t_i) E[\frac{V(T)}{N(T)} | F_{t_i} ]$ at time $t_i$ .

Lets say we use 100.000 samples in our simulation, so a bermudan swaption with two exercise dates would require 100.000 x 100.000 samples. Which is very time consuming and grows exponential with the number of exercise dates.

Instead of calculate the continuation value also with a Monte-Carlo simulation we will use a approximation. We use the approximation and algorithm developed by Longstaff and Schwarz.

We approximate the continuation value through some function of some state variable (e.g swap rate, short rate, etc…) $x_i$ :

$C(t_i) \approx f(x_i)$ .

A common choice for the function is of the form:

$f(x) = a_0 + a_1 g_1(x)+ a_2 g_2(x) + \dots a_n g_n(x)$

with $g_i (x)$ a polynom of the degree $i$ .

Lets choice $g_i (x)= x^{i-1}$ and $n=5$ for our example.

The coefficients of this function f are estimated by the ordinary least square error method. Thats why some people call this method also the OLS Monte-Carlo Method.

So how does the the algorithm look like?

We would simulate n paths of our stochastic process (in our case the short rate process) and go backward through time (backward induction).

We start at the last exercise date $T$ and calculate for each path the deflated terminal value of the swaption $V_T$ (exactly like we did in the european case).

Then we go back into time to the next exercise date and approximate the needed continuation value by our ordinary least square estimate.

Our choice for the state variable is the state variable of the short rate process itself.

We choice our coefficient so that the square error of $f(x_i) - V(T)$ is minimized, where $x_i$ is a vector of all simulated states over all paths and $V(T)$ the corresponding vector of the calculated npvs at time T from the previous step.

Now we use this coefficients to calculate the continuation value on each path by inserting the state of the current path into the formula.

On path j the deflated value of the swaption will be

$V_{i,j} = \frac{1}{N(t-i, x_{i,j})} max(S(t_i, x_{i,j}), N(t_i, x_{i,j}) C(x_{i,j})$

where $x_{i,j}$ is the the state of the process at time i on path j.

The npv at time 0 is then

$V(0) = N(0) \frac{1}{n} \sum_j=0^n V_{i,j}$ .

Some remarks:

One could also perform a 2nd MC Simulation (forward in time). First estimate the needed coefficients with a smaller numbers of paths. After train the model with this small set generate a larger sets of paths and go straight forward in time and use the approximation of the continuation value of the first run. On each path you only need to go forward until the holder will exercise the option.
The exercise time by this approach is only approximatively optimal, since its uses a approximation of the continuation value. Therefore our price will be an (asymptotic) lower bound of the real price.
The choice of appropriate functions g_i and the degree n is not obvious and can be tricky (see below).

Now lets have a look how this algorithm could be implemented in Python and Quantlib.

We use the notebook from my previous post as our starting point. We use the same yield curves, model (Gaussian short rate model) and the same underlying swap. The underlying 5y swap starts in 1Y.

We setup a bermudan swaption with 2 exercise dates, the first exercise date is on the start date of the swap and the 2nd date on the start date of the 2nd fixed leg accrual period.

Lets add both dates in our list of exercise dates:

calldates = [ settlementDate, 
              euribor6m.valueDate(ql.Date(5,4,2017))
            ]

We need to evaluate the value of the underlying swap at both exercise dates on each simulated path. To do that we introduce two new auxiliary functions. The first gives us all future (relative to an arbitrary evaluation time t) payment times and amounts from the fixed rate leg of the swap:

def getFixedLeg(swap, t):
    """
    returns all future payment times and amount of the fixed leg of the underlying swap

    Parameter:
        swap (ql.Swap)
        t (float) 

    Return:
        (np.array, np.array) (times, amounts)

    """
    fixed_leg = swap.leg(0)
    n = len(fixed_leg)
    fixed_times=[]
    fixed_amounts=[]
    npv = 0
    for i in range(n):
        cf = fixed_leg[i]
        t_i = timeFromReference(cf.date())
        if t_i > t:
            fixed_times.append(t_i)
            fixed_amounts.append(cf.amount())
    return np.array(fixed_times), np.array(fixed_amounts)

The 2nd will give us all future payment times, accrual start and end times, notionals, gearings and day count fractions of the floating leg. We need all this information to estimate the fixing of the float leg.

def getFloatingLeg(swap, t):

    float_leg = swap.leg(1)
    n = len(float_leg)
    float_times = []
    float_dcf = []
    accrual_start_time = []
    accrual_end_time = []
    nominals = []
    for i in range(n):
        # convert base classiborstart_idx Cashflow to
        # FloatingRateCoupon
        cf = ql.as_floating_rate_coupon(float_leg[i])
        value_date = cf.referencePeriodStart()
        t_fix_i = timeFromReference(value_date)
        t_i = timeFromReference(cf.date()) 
        if t_fix_i >= t:
            iborIndex = cf.index()

            index_mat = cf.referencePeriodEnd()
            # year fraction
            float_dcf.append(cf.accrualPeriod())
            # calculate the start and end time
            accrual_start_time.append(t_fix_i)
            accrual_end_time.append(timeFromReference(index_mat))
            # payment time
            float_times.append(t_i)
            # nominals 
            nominals.append(cf.nominal())
    return np.array(float_times), np.array(float_dcf), np.array(accrual_start_time), np.array(accrual_end_time), np.array(nominals)

With these two function we can evaluate the the underlying swap given the time and state of the process:

def swapPathNPV(swap, t):
    fixed_times, fixed_amounts = getFixedLeg(swap, t)
    float_times, float_dcf, accrual_start_time, accrual_end_time, nominals = getFloatingLeg(swap, t)
    df_times = np.concatenate([fixed_times, 
                           accrual_start_time, 
                           accrual_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, accrual_start_time, True)
    accrual_start_idx = accrual_start_idx.nonzero()
    accrual_end_idx = np.in1d(df_times, accrual_end_time, True)
    accrual_end_idx = accrual_end_idx.nonzero()
    # Calculate NPV
    def calc(x_t):
        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
        return npv
    return calc

This functions returns us a pricing function, which takes the current state of the process to price calculate the NPV of the swap at time t. The technique we use here is called closure function. The inner function is aware of the underlying swap and the ‘fixed’ time t and can access all variables defined in the outer function.

Remark on the pricing function
This is a very simple function. Its not possible to calculate the correct NPV of the swap for a time t between to fixing times of the floating leg. The current floating period will be ignored. So use that function only to evaluate the swap NPV only on fixing dates. We will extend this function to be capable to evaluate the NPV on any time t in the next post.

Use of the function at time t=0

The generation of the time grid hasn’t changed from the last time. It just consists of three points.

Also the generation of our sample path is the same as last time. After we generate our path we calculate the deflated payoffs of our swap at time T:

pricer = np.vectorize(swapPathNPV(swap, time_grid[-1]))
cont_value = pricer(y[:,-1]) / numeraires[:,-1]
cont_value[cont_value < 0] = 0

First we generate a vectorised function of the pricing function and use it on the array of our sample paths y and then apply the maximum function on the result.

In the next step we go one step back in time and calculate the deflated exercise value of the swaption at that time:

pricer = np.vectorize(swapPathNPV(swap, time_grid[-2]))
exercise_values = pricer(y[:,-2]) / numeraires[:,-2]
exercise_values[exercise_values < 0] = 0

Now we estimate the coefficients of continuation value function. We use the library statsmodels and fit an OLS model to the data.

states = y[:, -2]
Y = np.column_stack((states, states**2, states**3, states**4))
Y = sm.add_constant(Y)
ols = sm.OLS(cont_value, Y)
ols_result = ols.fit()

With this coefficients we can now calculate the continuation value on each path, given the state:

cont_value_hat = np.sum(ols_result.params * Y, axis=1)

The deflated value of the swaption at the first exercise is the maximum out of exercise value and continuation value:

npv_amc = np.maximum(cont_value_hat, exercise_values)

The npv at time 0 is the mean of the simulated deflated npvs at the first exercise date times the value of the numeraire at time 0:

npv_amc = np.mean(npv_amc) * numeraires[0,0]

To check the quality of our regression function $f$ we can have a look on a scatter plot:

As we can see the regression function doesn’t seems to fit that good to the left tail. So we could either increase the degree of our function, try other polynomial function, change to another state variable or try piecewise regression functions.

As usual you can download the source code from my github account or find it on nbViewer.

In the next post we are going to use this regression based approach to generate exposure paths for a multi callable swap.

I hope you enjoy the post. Till next time.

Jupyter notebooks – a Swiss Army Knife for Quants

A blog about quantitative finance, data science in fraud detection, machine and deep learning by Matthias Groncki

Exposure Simulation Part III / CVA for Bermudan Swaptions

Exposure Simulation / CVA and PFE for multi-callable swaps Part II