CVA Calculation with QuantLib and Python

Today I am going to present a way to calculate the credit value adjustment (CVA) for a netting set of plain vanilla interest rate swaps. This Monte-Carlo method is based on the code example of my previous post about the expected exposure and PFE calculation and the first steps will be exactly the same.

What is the credit value adjustment (CVA)?

The credit value adjustment is the difference between the risk-free price of a netting set and the the price which takes the possibility of the default of the counterparty into account. A netting set is a portfolio of deals with one counterparty for which you have a netting agreement. That means you are allowed to set positive against negative market values. A netting agreement will reduce your exposure and therefore the counterparty credit risk. If the time of default and value of the portfolio are independent then the CVA 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) at time t and the default probability PD(t).

The Monte Carlo approach

In my previous post we have seen a Monte Carlo method to estimate the expected exposure EE(t) at time t on an given time grid 0=t_0<\dots<T=t_n. We can reuse this estimator to 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)).

As we see in the formula we need two more ingredients to calculate the CVA, the discounted expected exposure and the the default probabilities of the counterparty.

Calculate the discounted expected exposure

To convert the expected exposure at time t in its present value expressed in time-zero dollars we only need to add a few more lines of code.

vec_discount_function = np.vectorize(t0_curve.discount)
discount_factors = vec_discount_function(time_grid)

We use the todays yield curve to calculate the discount factor for each point on our grid. With the numpy function vectorize
we generate a vectorized wrapper of the discount method of the QuantLib YieldTermStructure. This vectorized version can be applied on a array of ql.Dates or times instead of a scalar input. Basicalliy it’s equivalent to the following code snippet:

discount_factors = np.zeros(time_grid.shape)
for i in range(len(time_grid)):
	time = time_grid[i]
	discount_factors[i] = t0_curve.discount(time)

As the numpy documentation states the np.vectorize function is provided primarily for convenience, not for performance.

After the generating the market scenarios and pricing the netting set under each scenario we will calculate the discounted NPVs for each deal in the portfolio:

# Calculate the discounted npvs
discounted_cube = np.zeros(npv_cube.shape)
for i in range(npv_cube.shape[2]):
    discounted_cube[:,:,i] = npv_cube[:,:,i] * discount_factors
# Netting
discounted_npv = np.sum(discounted_cube, axis=2)

Discounted_NPV_Paths

Using the discounted npv cube we can calculate the discounted expected exposure:

# Calculate the exposure and discounted exposure
E = portfolio_npv.copy()
dE = discounted_npv.copy()
E[E&lt;0] = 0
EE = np.sum(E, axis=0)/N
dE[dE&lt;0] = 0
dEE = np.sum(dE, axis=0)/N

discounted_exposure

discounted_expected_exposure

The default probability of the counterparty

To derive the default probability one could either use market implied quotes (e.g. CDS) or use rating information (e.g. based on historical observations). We assume that the survival probability is given by sp(t)=\exp(-\int_0^t \lambda(t)) with a deterministic piecewise-constant hazard rate function \lambda(t). Given a grid of dates d_0<\dots<d_n and the corresponding backward flat hazard rate $latex\lambda_0,\dots,\lambda_n$ we can use the Quantlib HazardRateCurve to build a default curve.

# Setup Default Curve 
pd_dates =  [today + ql.Period(i, ql.Years) for i in range(11)]
hzrates = [0.02 * i for i in range(11)]
pd_curve = ql.HazardRateCurve(pd_dates,hzrates,ql.Actual365Fixed())
pd_curve.enableExtrapolation()

The QuantLib provides a real bunch of different types of DefaultTermStructures. You can either bootstrap a default curve from CDS quotes or you build a interpolated curve like we do here and combine one of the many interpolators (Linear, Backward Flat, etc.) with one of the possible traits (hazard rate, default probability, default density).

defaultCurve

With the default termstructure we can calculate the probability for a default between the times t_i and t_{i+1} for all i in our time grid.

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

Again we use the numpy function vectorize to apply a scalar function on an array. The method defaultProbability takes two times as input, t and T. It returns the probability of default between t and T.

Now we have all pieces together and the following code snippet gives us the CVA of our netting set:

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

You can download the code as an IPython (Juypter) notebook from here or just clone my repository (IPythonscripts) at GitHub.

If you want to read more about the QuantLib I would recommend to have a look on the blog and book “Implementing QuantLib” by Luigi. Another fantastic blog “Fooling around with QuantLib” by Peter has a very good and detailed post the Gsr model. Actually Peter has implemented this model in C++ and contributed it to the QuantLib.

I hope you have enjoyed reading my post and you will have fun playing around with the notebook. In one of my following posts I will extend this simulation by add a new product class to the netting set: European and bermudan swaptions.

So long.

Expected Exposure and PFE simulation with QuantLib and Python

In this post I will show how to use the Python bindings of the QuantLib library to calculate the expected exposure (EE) for a netting set of interest rate swaps in a IPython notebook. The technique I will present is very simple and works out of the box with standard QuantLib instruments and models. I will use a forward Monte Carlo Simulation to generate future market scenarios out of one-factor gaussian short rate model and evaluate the NPV of all swaps in the netting set under each scenario. The source code to this post (ExpectedExposureSimulation.ipynb) can be found in my repository IPythonScripts on GitHub or at nbviewer .

Methodology

First we define a time grid. On each date/time in our grid we want to calculate the expected exposure. For each date in our time grid we will simulate N states of the market and for each of these states we will calculate the NPV all of instruments in our portfolio / netting set. This results in N x (size of the netting set) simulated paths of NPVs. These paths can be used for EE, CVA (Credit value adjustment) or PFE (potential future exposure) calculations. swapPathsMc In the next step will we will floor each path at zero. This give the exposure of the portfolio on a path at each time. exposureSwap The expected exposure is given by the average of all paths: expExposureSwap The total number of NPV evaluations is (size of time grid) x (size of portfolio) x N. For a big portfolio and a very dense time grid it can be very time consuming task even if the single pricing is done pretty fast.

Assumption made in this example

For simplicity we restrict the portfolio to plain vanilla interest rate swaps in one currency. Further we assume that we live in a “single curve” world. We will use the same yield curve for discounting and forwarding. No spreads between the different tenor curves neither CSA discounting are taken into account. For the swap pricing we will need future states of the yield curve. In our setup we assume the the development of the yield curve follow an one factor Hull-White model. At the moment we make no assumption on how it is calibrated and assume its already calibrated. In our setting we will simulate N paths of the short rate following the Hull-White dynamics. At each time on each path the yield curve depend only on the state of our short rate process. We will use QuantLib functionalities to simulate the market states and perform the swap pricing on each path. The calculation of the expected exposure will be done in Python.

Technical Implementation

1. Setup of the market state at time zero (today)

rate = ql.SimpleQuote(0.03)
rate_handle = ql.QuoteHandle(rate)
dc = ql.Actual365Fixed()
yts = ql.FlatForward(today, rate_handle, dc)
yts.enableExtrapolation()
hyts = ql.RelinkableYieldTermStructureHandle(yts)
t0_curve = ql.YieldTermStructureHandle(yts)
euribor6m = ql.Euribor6M(hyts)

As mentioned above we live in a single curve world, we use a flat yield curve as discount and forward curve. During the Monte Carlo Simulation we will relink the Handle to the yieldTermStrucutre htys to our simulated yield curve. The original curve is stored in yts and the handle t0_curve.

2. Setup portfolio / netting set

# Setup a dummy portfolio with two Swaps
def makeSwap(start, maturity, nominal, fixedRate, index, typ=ql.VanillaSwap.Payer):
    &quot;&quot;&quot;
    creates a plain vanilla swap with fixedLegTenor 1Y
    
    parameter:
        
        start (ql.Date) : Start Date
        
        maturity (ql.Period) : SwapTenor
        
        nominal (float) : Nominal
        
        fixedRate (float) : rate paid on fixed leg
        
        index (ql.IborIndex) : Index
        
    return: tuple(ql.Swap, list&lt;Dates&gt;) Swap and all fixing dates
    
        
    &quot;&quot;&quot;
    end = ql.TARGET().advance(start, maturity)
    fixedLegTenor = ql.Period(&quot;1y&quot;)
    fixedLegBDC = ql.ModifiedFollowing
    fixedLegDC = ql.Thirty360(ql.Thirty360.BondBasis)
    spread = 0.0
    fixedSchedule = ql.Schedule(start,
                                end, 
                                fixedLegTenor, 
                                index.fixingCalendar(), 
                                fixedLegBDC,
                                fixedLegBDC, 
                                ql.DateGeneration.Backward,
                                False)
    floatSchedule = ql.Schedule(start,
                                end,
                                index.tenor(),
                                index.fixingCalendar(),
                                index.businessDayConvention(),
                                index.businessDayConvention(),
                                ql.DateGeneration.Backward,
                                False)
    swap = ql.VanillaSwap(typ, 
                          nominal,
                          fixedSchedule,
                          fixedRate,
                          fixedLegDC,
                          floatSchedule,
                          index,
                          spread,
                          index.dayCounter())
    return swap, [index.fixingDate(x) for x in floatSchedule][:-1]

The method makeSwap create a new QuantLib plain vanilla swap (see my previous post). We use this method to setup a netting set with two swaps:

portfolio = [makeSwap(today + ql.Period(&quot;2d&quot;),
                      ql.Period(&quot;5Y&quot;),
                      1e6,
                      0.03,
                      euribor6m),
             makeSwap(today + ql.Period(&quot;2d&quot;),
                      ql.Period(&quot;4Y&quot;),
                      5e5,
                      0.03,
                      euribor6m,
                      ql.VanillaSwap.Receiver),
            ]

Our netting set consists of two swaps, one receiver and one payer swap. Both swaps differ also in notional and time to maturity. Finally we create a pricing engine and link each swap in our portfolio with it.

engine = ql.DiscountingSwapEngine(hyts)

for deal, fixingDates in portfolio:
 deal.setPricingEngine(engine)
 print(deal.NPV())

In our Monte Carlo Simulation we can relink the handle hyts and use the same pricing engine. So we don’t need to create new pricing engines or relink the the deals to a new engine. We just need to call the method NPV of the instruments after relinking the yield term structure handle.

3. Monte-Carlo-Simulation of the “market”

We select a weekly time grid, including all fixing days of the portfolio. To generate the future yield curves we are using the GSR model and process of the QuantLib.

volas = [ql.QuoteHandle(ql.SimpleQuote(0.0075)),
ql.QuoteHandle(ql.SimpleQuote(0.0075))]
meanRev = [ql.QuoteHandle(ql.SimpleQuote(0.02))]
model = ql.Gsr(t0_curve, [today+100], volas, meanRev, 16.)
process = model.stateProcess()

The GSR model allows the mean reversion and the volatility to be piecewise constant. In our case here both parameter are set constant. For a more detailed view on the GSR model have a look on the C++ examples “Gaussian1dModels” in the QuantLib or here. Given a time t_0 and state x(t_0) of the process we know the conditional transition density for x(t_1) for t_1 > t_0. Therefore we don’t need to discretize the process between the evaluation dates. As a random number generator we are using the Mersenne Twister.

#%%timeit
# Generate N paths
N = 1500
x = np.zeros((N, len(time_grid)))
y = np.zeros((N, len(time_grid)))
pillars = np.array([0.0, 0.5, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10])
zero_bonds = np.zeros((N, len(time_grid), 12))
for j in range(12):
    zero_bonds[:, 0, j] = model.zerobond(pillars[j],
                                         0,
                                         0)
for n in range(0,N):
    dWs = generator.nextSequence().value()
    for i in range(1, len(time_grid)):
        t0 = time_grid[i-1]
        t1 = time_grid[i]
        x[n,i] = process.expectation(t0, 
                                     x[n,i-1], 
                                     dt[i-1]) + dWs[i-1] * process.stdDeviation(t0,
                                              x[n,i-1],
                                              dt[i-1])
        y[n,i] = (x[n,i] - process.expectation(0,0,t1)) / process.stdDeviation(0,0,t1)
        for j in range(12):
            zero_bonds[n, i, j] = model.zerobond(t1+pillars[j],
                                                 t1,
                                                 y[n, i])
                                                 

We also save the zero bonds prices on each scenario for a set of maturities (6M, 1Y,…,10Y). We use this prices as discount factors for our scenario yield curve.

4. Pricing on path & netting

On each date t and on each path p we will evaluate the netting set. First we build a new yield curve using the scenario discount factors from the step before.

date = date_grid[t]
ql.Settings.instance().setEvaluationDate(date)
ycDates = [date, 
           date + ql.Period(6, ql.Months)] 
ycDates += [date + ql.Period(i,ql.Years) for i in range(1,11)]
yc = ql.DiscountCurve(ycDates, 
                      zero_bonds[p, t, :], 
                      ql.Actual365Fixed())
yc.enableExtrapolation()
hyts.linkTo(yc)

After relinking the yield termstructure handle is the revaluation of the portfolio is straight forward. We just need to take fixing dates into account and store the fixings otherwise the pricing will fail.

for i in range(len(portfolio)):
    npv_cube[p, t, i] = portfolio[i][0].NPV()
5. Calculation EE and PFE

After populating the cube of fair values (1st dimension is simulation number, 2nd dimension is the time and 3rd dimension is the deal) we can calculate the expected exposure and the potential future exposure.

# Calculate the portfolio npv by netting all NPV
portfolio_npv = np.sum(npv_cube,axis=2)
# Calculate exposure
E = portfolio_npv.copy()
E[E&amp;amp;amp;lt;0]=0
EE = np.sum(E, axis=0)/N

Expected Expsoure Netting Set With PFE(t) we mean the 95 % quantile of the empirical distribution of the portfolio exposure at time t.

PFE_curve = np.apply_along_axis(lambda x: np.sort(x)[0.95*N],0, E)
MPFE = np.max(PFE_curve)

PFE netting set

Conclusion

With very few lines of code you can build a simulation engine for exposure profiles for a portfolio of plain vanilla swaps. These paths allows us to calculate the expected exposure or potential future exposure. Of course is the problem set in real world applications much more complex, we haven’t covered calibration (historical or market implied) or multi currencies / multi asset classes yet. And its getting even more complicated if you have multi-callable products in your portfolio.

But nevertheless I hope you have enjoyed reading this little tutorial and got an first insight into exposure simulation with QuantLib and Python. In one of my next post I will may extend this example about CVA calculation.

So long!

A brief introduction to the QuantLib in Python…

QuantLib is an open-source framework for quantitative finance  written in C++. There is an active community who develop and extend the library. QuantLib covers a wide range of financial instruments and markets like IR, FX and Equities and provide pricing engines and models, optimization algorithm, a Monte-Carlo framework, business day conventions, day count conventions, holidays calendars and even more… and grows continuously.

There are existing interfaces to different language like Python, Ruby, Java and more via swig. In the following I will describe how to build the QuantLib for Python.

In the following description I assume that you have already installed Python with IPython including all dependencies, numpy, scipy, pandas and matplotlib. A good starting point are distributions like Anaconda or WinPython (for Windows only). These distribution are shipped with many more useful libraries and provide an easy way to install more packages. Its worth to check them out, if haven’t tried them yet.

For windows you can find pre-build executables of the QuantLib at Christoph Gohlke’s website. If you want to have the most recent version you should clone the current master branch from GitHub and build the library and the python package yourself.

To build QuantLib from sources under Windows you need the correct Visual Studio depending on the Python version you are targeting. For Python 3.3 you need Visual Studio 2010 and for Python 2.7 you need Visual Studio 2008. Additional dependencies are the boost library and swig.

I will not cover how to install the dependencies but you can find binary installer for swig under windows on the swig project website and there also exists precompiled binaries for the boost library on source forge. Add the path to swig command to your PATH environment variable.

The following steps are valid for Python 2.7 and 3.3:

0. Clone the git repository either with using git:

cd c:\path\to\quantlib
git clone https://github.com/lballabio/quantlib .

or download as a zipped file from GitHub.

1. Open the QuantLib_vcXX.sln and build it in “Release” or “Release static runtime” configuration. For more details check the install documentation on the QuantLib project site. You will find the solution under

c:\path\to\quantlib\QuantLib

2. Open command line window and set required environment variables

SET QL_DIR=c:\path\to\quantlib\QuantLib
SET INCLUDE=c:\path\to\boost;%INCLUDE%
SET LIB=c:\path\to\boost\lib;%LIB%

3. Create the interface code with swig

cd c:\path\to\quantlib\QuantLib-SWIG\python
python setup.py wrap

4. Build and install the package

python setup.py build
python setup.py install

5. Test your install

python setup.py test

If you use Visual Studio 2010 and target Python 2.7 you will encounter the error message “Unable to find vcversall.bat“. The following “hack” works for me:

Set an additional environment variable at step 2:

SET VS90COMNTOOLS=%VS100COMNTOOLS%

Under unix you can also find some pre-build packages (e.g. for Ubuntu) in your package-manager or you build it from sources yourself with the usual configure, make and make install chain (see here and for Mac OS X here). After you build the library you can continue with step 3 above. Be sure you have installed swig before. Just download the source from the swig project and install via the usual unix install commands.

After we have managed to install QuantLib Python package I want to show you how you can use the QuantLib C++ classes from Python. You can download the following simple example  from my repository “IPythonScripts” on GitHub.

In the first example we will see how to setup an interest rate swap and use a discounting pricing engine to calculate the NPV and the fair rate for this swap.

First we need to import the QuantLib package, since we want also produce some plots we import numpy and matplotlib.

import numpy as np
import matplotlib.pyplot as plt
import QuantLib as ql

In the next step we will setup the yield curve (aka yieldtermstructure). For simplicity we will use an flat forward curve. I will use the same curve for discounting and calculation of the forwards. QuantLib supports multi curve setups, as I will maybe show in a later post or you can see in the examples how to setup more realistic yield curves.


# Set Evaluation Date
today = ql.Date(31,3,2015)
ql.Settings.instance().setEvaluationDate(today)
# Setup the yield termstructure
rate = ql.SimpleQuote(0.03)
rate_handle = ql.QuoteHandle(rate)
dc = ql.Actual365Fixed()
disc_curve = ql.FlatForward(today, rate_handle, dc)
disc_curve.enableExtrapolation()
hyts = ql.YieldTermStructureHandle(disc_curve)

The yieldTermStructure object provides an method which gives us the discount factor for a particular date (QuantLib.Date object) or time in years (with 0 = evaluationDate). This method is called discount.  I am using the numpy method vectorize to apply this function on arrays or list of times and then generate a plot of the discount curve.

discount = np.vectorize(hyts.discount)
tg = np.arange(0,10,1./12.)
plt.plot(tg, discount(tg))
plt.xlabel(&amp;quot;time&amp;quot;)
plt.ylabel(&amp;quot;discount factor&amp;quot;)
plt.title(&amp;quot;Flat Forward Curve&amp;quot;)

discountcurve

In the next step we will setup up a plain vanilla EURIBOR 6M Swap with maturity in 10 years. Therefore we generate an index and using the handle to our yield curve as forward curve and two schedules, one for the fixed rate leg with annual payments and one for the float leg with semi annual payments.

start = ql.TARGET().advance(today, ql.Period(&quot;2D&quot;))

end = ql.TARGET().advance(start, ql.Period(&quot;10Y&quot;))

nominal = 1e7

typ = ql.VanillaSwap.Payer

fixRate = 0.03

fixedLegTenor = ql.Period(&quot;1y&quot;)

fixedLegBDC = ql.ModifiedFollowing

fixedLegDC = ql.Thirty360(ql.Thirty360.BondBasis)

index = ql.Euribor6M(ql.YieldTermStructureHandle(disc_curve))

spread = 0.0

fixedSchedule = ql.Schedule(start, end, fixedLegTenor, index.fixingCalendar(), fixedLegBDC, fixedLegBDC, ql.DateGeneration.Backward, False)
floatSchedule = ql.Schedule(start, end, index.tenor(), index.fixingCalendar(), index.businessDayConvention(), index.businessDayConvention(), ql.DateGeneration.Backward, False)

swap = ql.VanillaSwap(typ, nominal, fixedSchedule, fixRate, fixedLegDC, floatSchedule, index, spread, index.dayCounter())

The last step before we can calculate the NPV we need a pricing engine. We are going to use the discountingSwapEngine. As the name suggests, it will discount all future payments to the evaluation date and calculate the difference between the present values of the two legs.

engine = ql.DiscountingSwapEngine(ql.YieldTermStructureHandle(disc_curve))
swap.setPricingEngine(engine)

print(swap.NPV())
print(swap.fairRate())

That was a very brief demonstration of the QuantLib. The examples which are shipped with the QuantLib give a much more detailed insight. So you may have a look on these.