Python-SciPy: Optimizing Smooth Libor Forward Curve

In my previous post, I presented how to use Python SciPy Optimization package for solving zero-coupon rate term structure from a given set of zero-coupon bond prices numerically. In this post, the same approach will be used in order to solve smooth Libor forward curve from a given set of vanilla swaps. Example of applying this scheme using Solver and Excel/VBA can be found in here.

In order to value fixed income derivatives cash flows, relevant forward rates and discount factors have to be defined from bootstrapped zero-coupon curve. In a Libor world, we use cash and FRA contracts (or futures contracts) in a short-end of the curve, while in a long-end of the curve we use swap contracts. Let us assume for a moment, that we bootstrap USD 3M zero-coupon curve, in order to get 3M forward rates on a quarterly basis. While bootstrapping approach usually gives smooth curve within short-end of the curve, it will usually fail to do this within long-end of the curve. This happens, because there is no swap contracts available in a long-end of the curve and hereby, we need to do a lot of interpolations for required swap rates. The resulting forward curve then looks like a chain of waterfalls and for the long-end of the curve, smoothness will be lost. This topic has been thoroughly covered in chapter three within this excellent book by Richard Flavell.

In a nutshell, we first define the shortest end of the curve (using spot cash rate) and then we just guess the rest of the forward rates. After this, we minimize sum of squared differences between adjacent forward rates, subject to constraints. In this case, constraints are present values of swaps, which by definition will be zero. When this optimization task has been performed, the resulting forward curve will successfully re-price the current swap market.

Bloomberg forward curve screen for the same set of data is available in here. Finally, thanks again for reading my blog.
-Mike

import numpy as np
import scipy.optimize as opt
import matplotlib.pyplot as pl

# sum of squared errors of all decision variables
# args: 0 = array of given initial rates, 1 = scaling factor
def ObjectiveFunction(x, args):
    # concatenate given initial rates and given 'guesses' for forward rates
    x = np.concatenate([args[0], x])
    return np.sum(np.power(np.diff(x), 2) * args[1])

# x = array of Libor forward rates
# args: 0 = swap rate, 1 = years to maturity, 2 = floating leg payments per year, 
# 3 = notional, 4 = array of given initial rates
def VanillaSwapPV(x, args):
    # PV (fixed payer swap) = -swapRate * Q(T) + (1 - DF(T))
    # where Q(T) is sumproduct of fixed leg discount factors and corresponding time accrual factors
    # assumption : fixed leg is always paid annually
    # since fixed leg is paid annually and maturity is always integer, 
    # time accrual factor will always be exactly 1.0
    
    # concatenate given initial rates and given 'guesses' for forward rates
    x = np.concatenate([args[4], x])
    DF = 0.0
    Q  = 0.0
    floatingLegTenor = 1 / args[2]
    nextFixedLegCouponDate = 1
    currentTimePoint = 0.0
    nCashFlows = int(args[1] * args[2])
    
    for i in range(nCashFlows):
        currentTimePoint += floatingLegTenor
        # first rate is always spot rate: calculate first spot discount factor
        if (i == 0):
            DF = 1 / (1 + x[i] * floatingLegTenor)
        # other rates are always forward rates
        if (i > 0):
            # solve current spot discount factor using previous spot discount factor
            # and current forward discount factor: DF(0,T) = DF(0,t) * DF(t,T), where (0 < t < T)
            DF *= (1 / (1 + x[i] * floatingLegTenor))
            if ((currentTimePoint - nextFixedLegCouponDate) == 0):
                Q += DF
                nextFixedLegCouponDate += 1
                
    return (-args[0] * Q + (1.0 - DF)) * args[3]

# array of given initial rates (the first given rate has to be spot cash rate)
# here, we give only one: 3M spot Libor fixing rate ('cash-to-first-future')
initialRates = np.array([0.006276])
scalingFactor = 1000000.0

# vanilla swap pricing functions as constraints
# args: swapRate, yearsToMaturity, floatingLegPaymentFrequency, notional, array of given initial rates
swaps = ({'type': 'eq', 'fun': VanillaSwapPV, 'args': [[0.0078775, 1, 4, scalingFactor, initialRates]]},
        {'type': 'eq', 'fun': VanillaSwapPV, 'args': [[0.009295, 2, 4, scalingFactor, initialRates]]},
        {'type': 'eq', 'fun': VanillaSwapPV, 'args': [[0.0103975, 3, 4, scalingFactor, initialRates]]},
        {'type': 'eq', 'fun': VanillaSwapPV, 'args': [[0.01136, 4, 4, scalingFactor, initialRates]]},
        {'type': 'eq', 'fun': VanillaSwapPV, 'args': [[0.012268, 5, 4, scalingFactor, initialRates]]},
        {'type': 'eq', 'fun': VanillaSwapPV, 'args': [[0.013183, 6, 4, scalingFactor, initialRates]]},
        {'type': 'eq', 'fun': VanillaSwapPV, 'args': [[0.01404, 7, 4, scalingFactor, initialRates]]},
        {'type': 'eq', 'fun': VanillaSwapPV, 'args': [[0.014829, 8, 4, scalingFactor, initialRates]]},
        {'type': 'eq', 'fun': VanillaSwapPV, 'args': [[0.01554, 9, 4, scalingFactor, initialRates]]},
        {'type': 'eq', 'fun': VanillaSwapPV, 'args': [[0.016197, 10, 4, scalingFactor, initialRates]]},
        {'type': 'eq', 'fun': VanillaSwapPV, 'args': [[0.016815, 11, 4, scalingFactor, initialRates]]},
        {'type': 'eq', 'fun': VanillaSwapPV, 'args': [[0.017367, 12, 4, scalingFactor, initialRates]]},
        {'type': 'eq', 'fun': VanillaSwapPV, 'args': [[0.018645, 15, 4, scalingFactor, initialRates]]},
        {'type': 'eq', 'fun': VanillaSwapPV, 'args': [[0.01996, 20, 4, scalingFactor, initialRates]]},
        {'type': 'eq', 'fun': VanillaSwapPV, 'args': [[0.020615, 25, 4, scalingFactor, initialRates]]},
        {'type': 'eq', 'fun': VanillaSwapPV, 'args': [[0.02097, 30, 4, scalingFactor, initialRates]]},
        {'type': 'eq', 'fun': VanillaSwapPV, 'args': [[0.021155, 40, 4, scalingFactor, initialRates]]},
        {'type': 'eq', 'fun': VanillaSwapPV, 'args': [[0.021017, 50, 4, scalingFactor, initialRates]]})

# initial guesses for Libor forward rates up to 50 years
# number of rates: 50 * 4 - 1 (years to maturity * floating leg payments per year - first given cash spot rate)
initialGuesses = np.full(199, 0.006)
model = opt.minimize(ObjectiveFunction, initialGuesses, args = ([initialRates, scalingFactor]), method = 'SLSQP', options = {'maxiter': 500}, constraints = swaps)

# print selected model results
print('Success: ' + str(model.success))
print('Message: ' + str(model.message))
print('Number of iterations: ' + str(model.nit))
print('Objective function (sum of squared errors): ' + str(model.fun))
#print('Changing variables (Libor forward rates): ' + str(model.x * 100))
pl.plot(model.x)
pl.show()

Python-SciPy: Solving Zero-Coupon Term Structure Using Optimization

This post is presenting how to use Python SciPy Optimization package for solving out zero-coupon rate term structure from a given set of zero-coupon bond prices. Needless to say, we do not need any numerical method to do this, since we have exact analytical formulas for backing out zero-coupon rates from zero-coupon bond prices. However, approach presented in this relatively simple example can be applied into more interesting schemes, such as solving out smooth Libor forward curve from a given set of vanilla swaps. Example of such scheme using Solver and Excel/VBA can be found in here.

Let us assume zero-coupon rate term structure as shown below. Let us also assume, that the only information we can observe is the zero-coupon bond prices. Our task is to solve zero-coupon rates term structure.

In a nutshell, we first set initial guesses for ten zero-coupon rates. After this we minimize sum of squared differences between adjacent zero-coupon rates (in objective function), but subject to constraints. In this case, constraints are differences between observed zero-coupon bond market prices and calculated bond prices, which need to be pushed to zero. When this optimization task has been completed, the resulting zero-coupon term structure will successfully re-price the current zero-coupon bond market, while curve smoothness is maximized. Example program results are shown below.

Thanks for reading this blog.
-Mike

import numpy as np
import scipy.optimize as opt
import matplotlib.pyplot as pl

# sum of squared errors of decision variables
def ObjectiveFunction(x, args):
    return np.sum(np.power(np.diff(x), 2) * args[0])

# zero coupon bond pricing function
def ZeroCouponBond(x, args):
    # return difference between calculated and market price
    return ((1 / (1 + x[args[3]])**args[2]) - args[0]) * args[1]

# zero coupon bond pricing functions as constraints
# args: market price, scaling factor, maturity, index number for rate array
zeroPrices = ({'type': 'eq', 'fun': ZeroCouponBond, 'args': [[0.998801438274071, 1000000.0, 1.0, 0]]},
            {'type': 'eq', 'fun': ZeroCouponBond, 'args': [[0.996210802629012, 1000000.0, 2.0, 1]]},
            {'type': 'eq', 'fun': ZeroCouponBond, 'args': [[0.991943543964159, 1000000.0, 3.0, 2]]},
            {'type': 'eq', 'fun': ZeroCouponBond, 'args': [[0.981028206597786, 1000000.0, 4.0, 3]]},
            {'type': 'eq', 'fun': ZeroCouponBond, 'args': [[0.962851266220459, 1000000.0, 5.0, 4]]},
            {'type': 'eq', 'fun': ZeroCouponBond, 'args': [[0.946534719794057, 1000000.0, 6.0, 5]]},
            {'type': 'eq', 'fun': ZeroCouponBond, 'args': [[0.924997530805076, 1000000.0, 7.0, 6]]},
            {'type': 'eq', 'fun': ZeroCouponBond, 'args': [[0.912584111300984, 1000000.0, 8.0, 7]]},
            {'type': 'eq', 'fun': ZeroCouponBond, 'args': [[0.892632531026722, 1000000.0, 9.0, 8]]},
            {'type': 'eq', 'fun': ZeroCouponBond, 'args': [[0.877098137542374, 1000000.0, 10.0, 9]]})

# initial guesses for ten zero-coupon rates
initialGuess = np.full(10, 0.005)
model = opt.minimize(ObjectiveFunction, initialGuess, args = ([1000000.0]), method = 'SLSQP', constraints = zeroPrices)

# print selected model results
print('Success: ' + str(model.success))
print('Message: ' + str(model.message))
print('Number of iterations: ' + str(model.nit))
print('Objective function (sum of squared errors): ' + str(model.fun))
print('Changing variables (zero-coupon rates): ' + str(model.x * 100))
pl.plot(model.x)
pl.show()

Python: Path Generator for Correlated Processes

One reader was interested to know, how to simulate correlated asset paths by using just Python libraries, without using QuantLib. This blog post is presenting the result of woodshedding this stuff. A couple of notes:

• GeneratePaths method can be used to simulate paths for a single process or multiple processes, based on a given (or non-existing) correlation matrix.
• Method returns Numpy array having dimensions based on the given number of processes, number of paths and number of time steps.
• Discretized stochastic processes for generator method to be simulated are defined as lambda methods. This approach makes it relatively easy to implement (almost) any desired single-factor model.

The program is enough commented and should be self-explainable. Thanks for reading this blog.
-Mike

import numpy as np
import matplotlib.pyplot as pl

# returns ndarray with the following dimensions: nProcesses, nPaths, nSteps
def GeneratePaths(spot, process, maturity, nSteps, nPaths, correlation = None):
    dt = maturity / nSteps
    
    # case: given correlation matrix, create paths for multiple correlated processes
    if (isinstance(correlation, np.ndarray)):
        nProcesses = process.shape[0]
        result = np.zeros(shape = (nProcesses, nPaths, nSteps))
        
        # loop through number of paths
        for i in range(nPaths):
            # create one set of correlated random variates for n processes
            choleskyMatrix = np.linalg.cholesky(correlation)
            e = np.random.normal(size = (nProcesses, nSteps))            
            paths = np.dot(choleskyMatrix, e)
            # loop through number of steps
            for j in range(nSteps):
                # loop through number of processes
                for k in range(nProcesses):
                    # first path value is always current spot price
                    if(j == 0):
                        result[k, i, j] = paths[k, j] = spot[k]
                    else:
                        # use SDE lambdas (inputs: previous spot, dt, current random variate)
                        result[k, i, j] = paths[k, j] = process[k](paths[k, j - 1], dt, paths[k, j])

    # case: no given correlation matrix, create paths for a single process
    else:
        result = np.zeros(shape = (1, nPaths, nSteps))
        # loop through number of paths
        for i in range(nPaths):
            # create one set of random variates for one process
            path = np.random.normal(size = nSteps)
            # first path value is always current spot price
            result[0, i, 0] = path[0] = spot
            # loop through number of steps
            for j in range(nSteps):
                if(j > 0):
                    # use SDE lambda (inputs: previous spot, dt, current random variate)
                    result[0, i, j] = path[j] = process(path[j - 1], dt, path[j])
    return result

# Geometric Brownian Motion parameters
r = 0.03
v = 0.25

# define lambda for process (inputs: spot, dt, random variate)
BrownianMotion = lambda s, dt, e: s + r * s * dt + v * s * np.sqrt(dt) * e   

# general simulation-related parameters
maturity = 1.0
nPaths = 10
nSteps = 250

# case: one process
SingleAssetPaths = GeneratePaths(100.0, BrownianMotion, maturity, nSteps, nPaths)
for i in range(nPaths):
    pl.plot(SingleAssetPaths[0, i, :])
pl.show()

# case: two correlated processes
matrix = np.array([[1.0, 0.999999], [0.999999, 1.0]])
spots = np.array([100.0, 100.0])
processes = np.array([BrownianMotion, BrownianMotion])
MultiAssetPaths = GeneratePaths(spots, processes, maturity, nSteps, nPaths, matrix)
f, subPlots = pl.subplots(processes.shape[0], sharex = True)
for i in range(processes.shape[0]): 
    for j in range(nPaths):
        subPlots[i].plot(MultiAssetPaths[i, j, :])
pl.show()

C#/Python: Creating Python Wrapper for C# Class by Using Python for .NET

Interoperability is just amazing concept. Sometimes, instead re-inventing the wheel again for a new language, it might be easier to recycle the old stuff. Creating Python wrapper for C++ program was presented in previous post. This post will present simple steps for creating Python wrapper for a simple C#.NET class.

Python for .NET

Python for .NET is a package, that gives Python programmers nearly seamless integration with .NET Common Language Runtime. With this package, you can script .NET applications or build entire applications in Python, using .NET services and components written in any language that targets the CLR (C#, VB.NET, F#, C++/CLI). Download package from here.

Create C# class

In Visual Studio, create a new C# class library project under the namespace CsLibrary. Create C# class Functions into this class project. The content of this class is presented below.

using System;

namespace CsLibrary
{
    // C# class
    public class Functions
    {

        public double Add(double a, double b) {
            return a + b;
        }

        public double Subtract(double a, double b) {
            return a - b;
        }

        public double Multiply(double a, double b) {
            return a * b;
        }

        public double Divide(double a, double b) {
            return a / b;
        }

    }
}

Create wrapper for C# class

In order to access this class from Python, we need to have wrapper class, which inherits from DynamicObject. Add the following new class to existing solution. Effectively, this class is just object adapter, which delegates function calls to our previous C# class.

using System;
using System.Dynamic;

namespace CsLibrary {

    public class PyWrapper : DynamicObject {

        // wrapped C# class
        private Functions functions;

        // ctor
        public PyWrapper() {
            functions = new Functions();
        }        
        
        public double Add(double a, double b) {
            return functions.Add(a, b);
        }

        public double Subtract(double a, double b) {
            return functions.Subtract(a, b);
        }

        public double Multiply(double a, double b) {
            return functions.Multiply(a, b);
        }

        public double Divide(double a, double b) {
            return functions.Divide(a, b);
        }
    }
}

Use wrapper class in Python

The following example program is using C# class wrapper (CsLibrary.dll) in Python.

import clr
clr.AddReference(r'..\CsLibrary\bin\Release\CsLibrary.dll')
from CsLibrary import PyWrapper
wrapper = PyWrapper()

x = 1.23
y = 2.34

print(wrapper.Add(x, y))
print(wrapper.Subtract(x, y))
print(wrapper.Multiply(x, y))
print(wrapper.Divide(x, y))

# 3.57
# -1.1099999999999999
# 2.8781999999999996
# 0.5256410256410257

That's it. The files can be downloaded from my GitHub page. Thanks for reading my blog.
-Mike

C++/Python: Creating Python Wrapper for C++ Class by Using SWIG

If the need sometimes arises, existing C++ libraries can be interfaced relatively easy to be used in Python by using SWIG wrapper. SWIG stands for Simplified Wrapper and Interface Generator and it is an open-source software tool used to connect computer programs or libraries written in C or C++ with scripting languages, such as Python. As one extremely motivating example, QuantLib C++ libraries have been made available for Python by using SWIG wrappers. Complete documentation for SWIG 3.0 can be found in here. The actual SWIG package can be downloaded in here. By taking a look at the documentation, one may understand quickly, that in this post only the tip of the iceberg will get scratched.

Also, it should be noted at this point, that this post is assuming all operations taking place in Linux environment. Corresponding instructions for Windows users can be found in SWIG documentation. This blog post is presenting, how to interface custom C++ random generator class to be used in Python. The original (template) class implementation can be found in here.

Header and implementation files for simple C++ random number generator class are presented below.

Header file (RG.h)

#include <algorithm>
#include <functional>
#include <random>

class Generator {
public:
 Generator(unsigned long seed);
 void operator()(std::vector<double>& v);
private:
 std::function<double(double)> randomGenerator;
 std::normal_distribution<double> distribution;
 std::mt19937 uniformGenerator;
};

Implementation file (RG.cpp)

#include "RG.h"

Generator::Generator(unsigned long seed) {
  // construct lambda method for processing standard normal random number
  randomGenerator = [this](double x)-> double {
   x = distribution(uniformGenerator);
   return x;
  };
  uniformGenerator.seed(seed);
}

// fill client-given vector with random variates from standard normal distribtuion
void Generator::operator()(std::vector<double>& v) {
  std::transform(v.begin(), v.end(), v.begin(), randomGenerator);
}

The class functionality is simple and straightforward. First, create a class implementation by giving desired seed value in constructor. In constructor, C++ function object will be created. This function object will ultimately create our normally distributed random variates. Client will request a new set of random variates by giving a reference for STL vector by using parenthesis operator, which has been overloaded here for syntactic reasons only.

SWIG interface file (RG.i)

A simple SWIG interface for our class can be built by simply wrapping the header file as follows.

%module RG
%{
#include "RG.h"
%}

%include "std_vector.i"
%template(DoubleVector) std::vector<double>;

%include "RG.h"

Python does not know anything about std::vector. The std_vector.i library provides support for C++ STL vector class. Then, all we need to do is to instantiate different versions of vector for all the specific types we would like to use in Python. For our example class, we want to use vector consisting of double data types. Using this library involves the use of the %template directive. In Python code, C++ vector will be used by using alias DoubleVector.

SWIG build

All in all we have now three files in our folder: RG.h, RG.cpp and RG.i. With terminal opened in the folder consisting the previous three files, the following three commands will create Python wrapper for our C++ class.

$ swig -c++ -python RG.i
$ g++ -fpic -c -I/home/mikejuniperhill/anaconda3/include/python3.7m RG.cpp RG_wrap.cxx
$ g++ -shared RG.o RG_wrap.o -o _RG.so

In the second command, a correct parameter for -I argument can be found by using the following command in terminal.

$ python3-config --cflags

After these three steps, our C++ program can be accessed and used in Python program. The following simple example program is requesting a vector of standard normal random numbers from our wrapped C++ class and uses those for calculating Monte Carlo PV for one European call option.

# import RG module
import SwigTester.RG as RG
import numpy as np

# set seed, create generator object
seed = 0
generator = RG.Generator(seed)

# set vector size, initialize vector
steps = 1000
vector = RG.DoubleVector(steps)

# request generator to fill vector with standard normal variates 
generator(vector)
# copy content to numpy array
e = np.array(vector)

# use variates to value one european call option
# set parameters
s0 = 100.0
x = 100.0
t = 1.25
v = 0.25
r = 0.01

# use vectorization in numpy, pre-calculate components
drift = (r - 0.5 * v * v) * t
diffusion = v * np.sqrt(t)
df = np.exp(-r * t)

# get option value as discounted average of all terminal payoffs
c0 = np.mean(np.maximum(s0 * np.exp(drift + diffusion * e) - x, 0.0)) * df
print(c0)

Note, that I imported SwigTester.RG, because I created Python wrapper module directly into this specific folder.

As a SWIG newcomer, one may expect to face all kinds of cryptic installation and other infuriating compatibility issues. As always, Google is the best friend along this cruel and unusual phase. But, have faith - and you'll be rewarded with all wonders and joys of SWIG wrapper.

All relevant files can be found directly from my GitHub repository. Thanks for reading my blog.
-Mike

Python: Market Scenario Files Generator for Third-party Analytics Software

Third-party analytics software usually requires specific set of market data for performing its calculations. In this post, I am publishing one of my utility Python programs for creating different types of stress scenario markets, based on given base market and set of prepared XML configurations. The complete program can be found in my GitHub repository.

Configurations

The following screenshot shows configurations for this program. SourceFilePath attribute captures the source market data CSV file and TargetFolderPath captures the folder, into which all market scenario files will be created. Finally, ScenarioConfigurationsPath captures the folder, which contains all XML scenario configuration files. This configuration XML file should be stored in a chosen directory.

<Configurations>
    <!-- attributes for scenario generator settings -->
    <ScenarioConfigurationsPath>\\Temp\ScenarioConfigurations\</ScenarioConfigurationsPath>
    <SourceFilePath>\\Temp\baseMarket.csv</SourceFilePath>
    <TargetFolderPath>\\Temp\Scenarios\</TargetFolderPath>
</Configurations>

Market data

The following screenshot shows given base market data. Due to brevity reasons, only EUR swap curve has been used here as an example. All market data points are defined here as key-value pairs (ticker, value).

IR.EUR-EURIBOR.CASH-1BD.MID,-0.00365
IR.EUR-EURIBOR.CASH-1W.MID,-0.00373
IR.EUR-EURIBOR.CASH-1M.MID,-0.00363
IR.EUR-EURIBOR.CASH-2M.MID,-0.00336
IR.EUR-EURIBOR.CASH-3M.MID,-0.00309
IR.EUR-EURIBOR.CASH-6M.MID,-0.00237
IR.EUR-EURIBOR-6M.FRA-1M-7M.MID,-0.00231
IR.EUR-EURIBOR-6M.FRA-2M-8M.MID,-0.00227
IR.EUR-EURIBOR-6M.FRA-3M-9M.MID,-0.002255
IR.EUR-EURIBOR-6M.FRA-4M-10M.MID,-0.00218
IR.EUR-EURIBOR-6M.FRA-5M-11M.MID,-0.00214
IR.EUR-EURIBOR-6M.FRA-6M-12M.MID,-0.002075
IR.EUR-EURIBOR-6M.FRA-7M-13M.MID,-0.00198
IR.EUR-EURIBOR-6M.FRA-8M-14M.MID,-0.00186
IR.EUR-EURIBOR-6M.FRA-9M-15M.MID,-0.001775
IR.EUR-EURIBOR-6M.FRA-10M-16M.MID,-0.00169
IR.EUR-EURIBOR-6M.FRA-11M-17M.MID,-0.00159
IR.EUR-EURIBOR-6M.FRA-12M-18M.MID,-0.00141
IR.EUR-EURIBOR-6M.SWAP-2Y.MID,-0.001603
IR.EUR-EURIBOR-6M.SWAP-3Y.MID,-0.000505
IR.EUR-EURIBOR-6M.SWAP-4Y.MID,0.00067
IR.EUR-EURIBOR-6M.SWAP-5Y.MID,0.00199
IR.EUR-EURIBOR-6M.SWAP-6Y.MID,0.003315
IR.EUR-EURIBOR-6M.SWAP-7Y.MID,0.0046
IR.EUR-EURIBOR-6M.SWAP-8Y.MID,0.00584
IR.EUR-EURIBOR-6M.SWAP-9Y.MID,0.00698
IR.EUR-EURIBOR-6M.SWAP-10Y.MID,0.00798
IR.EUR-EURIBOR-6M.SWAP-11Y.MID,0.00895
IR.EUR-EURIBOR-6M.SWAP-12Y.MID,0.009775
IR.EUR-EURIBOR-6M.SWAP-15Y.MID,0.01165
IR.EUR-EURIBOR-6M.SWAP-20Y.MID,0.013245
IR.EUR-EURIBOR-6M.SWAP-25Y.MID,0.013725
IR.EUR-EURIBOR-6M.SWAP-30Y.MID,0.01378

We can clearly see, that the system used for constructing market data tickers leads to scheme, in which every market data point will have one and only one unique ticker. This will then guarantee, that we can drill down and stress individual market data points with regex expressions, if so desired. This data should be copied into CSV file (directory has been defined in previous configuration file).

Scenario configurations

The following screenshot shows XML configurations for one market scenario. One such scenario can have several different scenario items (Say, stress these rates up, stress those rates down, apply these changes to all FX rates against EUR and set hard-coded values for all CDS curves). From these configurations, ID and description are self-explainable. Attribute regExpression captures all regex expressions (scenario items), which will be searched from risk factor tickers. As soon as regex match is found, the program will use corresponding operationType attribute to identify desired stress operation (addition, multiplication or hard-coded value). Finally, the amount of change which will be applied in risk factor value is defined within stressValue attribute. This XML configuration should be stored (directory has been defined in program configuration file).

<!-- operation types : 0 = ADDITION, 1 = MULTIPLICATION, 2 = HARD-CODED VALUE -->
<Scenario>
  <ID>CURVE.STRESS</ID>
  <description>custom stress scenario for EUR swap curve</description>
  <regExpression>^IR.EUR-EURIBOR.CASH,^IR.EUR-EURIBOR-6M.FRA,^IR.EUR-EURIBOR-6M.SWAP</regExpression>
  <operationType>1,0,2</operationType>
  <stressValue>1.25,0.015,0.05</stressValue>
</Scenario>

Finally, the following screenshot shows resulting market data, when all configured scenario items have been applied. This is the content of output CSV file, created by this Python program.

IR.EUR-EURIBOR.CASH-1BD.MID,-0.0045625
IR.EUR-EURIBOR.CASH-1W.MID,-0.0046625
IR.EUR-EURIBOR.CASH-1M.MID,-0.0045375
IR.EUR-EURIBOR.CASH-2M.MID,-0.004200000000000001
IR.EUR-EURIBOR.CASH-3M.MID,-0.0038624999999999996
IR.EUR-EURIBOR.CASH-6M.MID,-0.0029625000000000003
IR.EUR-EURIBOR-6M.FRA-1M-7M.MID,0.01269
IR.EUR-EURIBOR-6M.FRA-2M-8M.MID,0.01273
IR.EUR-EURIBOR-6M.FRA-3M-9M.MID,0.012745
IR.EUR-EURIBOR-6M.FRA-4M-10M.MID,0.01282
IR.EUR-EURIBOR-6M.FRA-5M-11M.MID,0.01286
IR.EUR-EURIBOR-6M.FRA-6M-12M.MID,0.012924999999999999
IR.EUR-EURIBOR-6M.FRA-7M-13M.MID,0.01302
IR.EUR-EURIBOR-6M.FRA-8M-14M.MID,0.013139999999999999
IR.EUR-EURIBOR-6M.FRA-9M-15M.MID,0.013224999999999999
IR.EUR-EURIBOR-6M.FRA-10M-16M.MID,0.013309999999999999
IR.EUR-EURIBOR-6M.FRA-11M-17M.MID,0.01341
IR.EUR-EURIBOR-6M.FRA-12M-18M.MID,0.01359
IR.EUR-EURIBOR-6M.SWAP-2Y.MID,0.05
IR.EUR-EURIBOR-6M.SWAP-3Y.MID,0.05
IR.EUR-EURIBOR-6M.SWAP-4Y.MID,0.05
IR.EUR-EURIBOR-6M.SWAP-5Y.MID,0.05
IR.EUR-EURIBOR-6M.SWAP-6Y.MID,0.05
IR.EUR-EURIBOR-6M.SWAP-7Y.MID,0.05
IR.EUR-EURIBOR-6M.SWAP-8Y.MID,0.05
IR.EUR-EURIBOR-6M.SWAP-9Y.MID,0.05
IR.EUR-EURIBOR-6M.SWAP-10Y.MID,0.05
IR.EUR-EURIBOR-6M.SWAP-11Y.MID,0.05
IR.EUR-EURIBOR-6M.SWAP-12Y.MID,0.05
IR.EUR-EURIBOR-6M.SWAP-15Y.MID,0.05
IR.EUR-EURIBOR-6M.SWAP-20Y.MID,0.05
IR.EUR-EURIBOR-6M.SWAP-25Y.MID,0.05
IR.EUR-EURIBOR-6M.SWAP-30Y.MID,0.05

Handy way to create and test regex expressions is to use any online tool available. As an example, the first scenario item (^IR.EUR-EURIBOR.CASH) has been applied to a given base market data. The last screenshot below shows all regex matches.

Have a great start for the year 2019 and thanks a lot again for reading my blog.
-Mike

QuantLib-Python: Multiprocessing Method Wrapper

In this post, I published a program for simulating term structure up to 30 years with daily time step, using Hull-White one-factor model. The resulting curve was able to replicate the given initial yield curve without any notable differences. The only serious issue here was the cost of running the program in time units.

In order to improve this specific issue, I got familiar with some multi-threading possibilities in Python. Now, there are some deep issues related to the way Python threads are actually implemented and especially the issues with thread locking known as GIL. Related stuff has been completely chewed in here. In order to avoid facing GIL-related issues, another way is to use Python multiprocessing, which allows the programmer to fully leverage multiple processors on a given machine. Moreover, separate processes are completely separate and one process cannot affect another's variables.

Comparison statistics

Again, I simulated term structure up to 30 years with daily time step, using 10000 paths. I did some quick profiling, in order to find the bottlenecks in the original program (sequential). By looking the column task share, we can see, that there are two tasks, which are consuming the most part of the complete running time: path generations (~54%) and path integrations (~46%). After this, I isolated these two parts and processed these by using multiprocessing scheme.

By using multiprocessing (two configured processes), I managed to decrease the complete running time from 163 seconds to 107 seconds. In general, for all those parts which were enjoying the benefits of multiprocessing, improvement ratio (multiprocessing per sequential) is around 0.65.

CPU puzzle

In order to understand the reason for this improvement, let us take a look at processor architecture in my laptop. First, let us check the "Grand Promise" made by System Monitor.

Based on this view, I actually expected to have four CPU for processing. I was then really surprised, that adding third and fourth process was not decreasing running time any further, than having just two processes. After some usual Stackoverflow gymnastics, I finally got the definition to calculate the real number of CPU available in my laptop.

CPU available is "Core(s) per socket * Socket(s)", which is 2 * 1 in my laptop. So, all in all I have only two CPU available for processing, not four as was shown in that System Monitor. This means, that having more than two CPU available in a laptop, one should expect even better improvement ratio than reported in my statistics here.

Wrapper method

In order to avoid code duplication and to come up with something a bit more generic, I started to dream about the possibility to create a mechanism for applying multiprocessing for any given method, if so desired. Such solution is possible by using Python lambda methods.

# method for executing given lambdas in parallel
def MultiprocessingWrapper(targetFunctionList):
    processList = []
    aggregatedResults = []
    queue = mp.Manager().Queue()

    # execute lambda from a given list based on a given index
    # storing results into queue
    def Worker(index):
        result = targetFunctionList[index]()
        queue.put(result)
    
    # start processes, call worker method with index number
    for i in range(len(targetFunctionList)):
        process = mp.Process(target = Worker, args = (i,))
        processList.append(process)
        process.start()
    
    # join processes, extract queue into results list
    for process in processList:
        process.join()
        aggregatedResults.append(queue.get())
        
    # return list of results for a client
    return aggregatedResults

Previous method is receiving a list of lambda methods as its argument. The idea is, that there would be always one lambda method for each process to be started. Why? We might face a situation, in which different set of parameters would be required for each lambda (say, all lambdas would be using the same uniform random generator, but with a different values for seeding the generator). In wrapper method, process is then created to start each configured lambda method. Results calculated by given lambda will be stored into queue. Finally, all results will be imported from queue into result list and returned for a client.

As concrete example, the following program segment is first creating two lambda methods for starting GeneratePaths method, but with different seed value for both processes. After this, wrapper method is then fed with the list of lambdas and processed paths will be received into list of results.

    # task of generating paths is highly time critical
    # use multiprocessing for path generation
    # create lambdas for multiprocessing wrapper
    # target signature: def GeneratePaths(seed, process, timeGrid, n)
    nPaths = 10000
    nProcesses = 2
    seeds = [1834, 66023]
    nPathsPerProcess = int(nPaths / nProcesses)
    target_1 = lambda:GeneratePaths(seeds[0], HW1F, grid.GetTimeGrid(), nPathsPerProcess)
    target_2 = lambda:GeneratePaths(seeds[1], HW1F, grid.GetTimeGrid(), nPathsPerProcess)
    targetFunctionList = [target_1, target_2]
    results = MultiprocessingWrapper(targetFunctionList)

The complete program can be found in my GitHub repository. Finally, a couple of discussion threads, which may be useful in order to understand some QuantLib-related issues, are given in here and here. Thanks for reading my blog and again, Merry Christmas for everyone.
-Mike
  

QuantLib-Python: Exposure Simulation

This Python program is using QuantLib library tools for simulating exposures for one selected Bloomberg vanilla benchmark swap transaction. Based on simulated exposures, the program will then calculate Expected Positive Exposure (EPE) and Expected Negative Exposure (ENE), as well as corresponding CVA and DVA statistics.

Then, some (unfortunate) limitations: the program can only handle one transaction at a time, so simulating exposures for netting sets having several transactions is not possible. Also, the program can simulate only one risk factor at a time, so simulating exposures for transactions exposed to more than one risk factor is not possible. However, with some careful re-designing, these properties could also be implemented by using QuantLib library tools.

The complete program can be found in my GitHub repository. Thanks for reading this blog. Merry Christmas for everyone.
-Mike


Simulated exposures

Data

A few notes on data.

• Swap transaction is 5Y receiver vs. 3M USD Libor + spread. At inception, swap PV has been solved to be zero. Details can be found in the screenshot below.

• Interest rate data for spot term structure (discount factors) has been retrieved from Bloomberg Swap Manager as of 12.12.2018.

• Default term structures for the both parties (counterparty, self) are created from flat CDS term structures (100 bps), as seen on Bloomberg Swap Manager CVA tab.

• Short rate simulations are processed by using Hull-White one-factor model, which uses parameters calibrated to a given set of flat 20% swaption volatilities, as seen on Bloomberg Swap Manager CVA tab.

Results

Bloomberg Swap Manager results: CVA = 6854 and DVA = 1557. Program results (for one run): CVA = 6727 and DVA = 1314, using weekly time steps and 1000 paths. However, "close enough" results can be achieved with considerably smaller amount of paths and less dense time grid.

Screens

Bloomberg swap transaction



Bloomberg CVA



Bloomberg DVA



Bloomberg EPE



Program EPE

Bloomberg ENE



Program ENE