# 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. 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. The expected exposure is given by the average of all paths: 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.

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()
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;
fixedLegTenor = ql.Period(&quot;1y&quot;)
fixedLegBDC = ql.ModifiedFollowing
fixedLegDC = ql.Thirty360(ql.Thirty360.BondBasis)
fixedSchedule = ql.Schedule(start,
end,
fixedLegTenor,
index.fixingCalendar(),
fixedLegBDC,
fixedLegBDC,
ql.DateGeneration.Backward,
False)
floatSchedule = ql.Schedule(start,
end,
index.tenor(),
index.fixingCalendar(),
ql.DateGeneration.Backward,
False)
swap = ql.VanillaSwap(typ,
nominal,
fixedSchedule,
fixedRate,
fixedLegDC,
floatSchedule,
index,
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,
]
```

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()
```

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].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
``` 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)
```

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!

## 8 thoughts on “Expected Exposure and PFE simulation with QuantLib and Python”

1. Luigi Ballabio says:

Hi Matthias,
a small issue: I think WordPress (or the plugin you use for code listings) messed up your quotes. You have stuff in there like

ql.Period(&quot;5Y&quot;)

Apart from that, interesting stuff. I just referred someone on the mailing list to your post…

Like

1. Matthias Groncki says:

Hi Luigi,
thanks for your comment, I will fix the issue with the quotes.

Like

2. Peter Caspers says:

I had this, too, WordPress changed these characters again and again, even if the preview looked fine. The only way out for me was to go back to the “classic editor” which doesn’t mess up things.

Like

3. Bernd says:

Dear Matthias,

I worked through your above example and got some problems in the line where you calculate the 95% confidence bands of future exposure paths. Of course I run the calculations also with the notebook downloaded from your page (exact copy of your code). But with the same error 😦

Do you have an idea what the problem might be?

Thank’s
Bernd

# Calculate the PFE curve (95% quantile)
PFE_curve = np.apply_along_axis(lambda x: np.sort(x)[0.95*N],0, E)

in (x)
1 # Calculate the PFE curve (95% quantile)
—-> 2 PFE_curve = np.apply_along_axis(lambda x: np.sort(x)[0.95*N],0, E)
3
4 plt.figure(figsize=(7,5), dpi=300)
5 plt.plot(time_grid,PFE_curve)

IndexError: only integers, slices (`:`), ellipsis (`…`), numpy.newaxis (`None`) and integer or boolean arrays are valid indices

Liked by 1 person

1. Matthias Groncki says:

Try to cast N*0.95 to an int. Like int(N*0.95)

Liked by 1 person

1. Bernd says:

That solved the problem. Thanks!

Liked by 1 person