MoCaX: Using Chebyshev Tensors for Computational Bottlenecks in Risk Calculations

Even this Machine Learning has been the hottest hype already for years, I have personally still had very vague understanding of how this thing could be applied in risk management domain. This post is presenting one such application: Chebyshev Tensors for getting rid of computational bottlenecks in risk calculations.

All risk management people are familiar with the following generic steps.

The sequence is describing steps for the most types of computationally-intensive risk calculations. In here, pricing (revaluation of financial transactions) is the step, which makes these calculations computationally intensive, since there will usually be substantial amount of revaluation calls sent to all kinds of pricing functions. For this reason, processing times will usually explode and this will costs us time.

In the past, the banks have usually tried to resolve these computational bottlenecks by acquiring more computational muscle and relying on multi-processing schemes. However, for a large enough player with an active derivatives book containing exotic products, even these attempts might not be enough to resolve the issue. Since the burden of producing all kinds of regulatory measures seems to be just growing, there is simply more calculations to be processed, than we have hours within one day.

Chebyshev Tensors

Another way to resolve such computational issues is to find clever algorithmic approximations and that's where Chebyshev Tensors are amazing. In a nutshell, by using Chebyshev Tensors we can create extremely accurate replicas of our standard pricing functions and use them instead of our original time-consuming function in calculations..

Interpolation methods through polynomials can deliver exponential convergence to the original function, when the interpolating points are Chebyshev points. Also, derivatives of the polynomial converge exponentially to the derivatives of the original function. Finally, Chebyshev Tensors can be extended to higher dimensions.

Example: 1-dimensional MoCaX approximation on valuing European call options

First, let us value a portfolio of European call options, where the spot price varies between 50 and 150. We divide this spot price range into one million steps. Other parameters in the model are constants. We get total of one million calculation tasks and one million revaluation calls to our pricing function. Note, that processing time is about 330 seconds (5.5 minutes) to complete all calculations. 

Next, we calculate the approximations for exactly the same set of option values, but we give only 25 analytically-valued option values for Mocax object and the rest of the values will be interpolated by Mocax object. Processing time is about 5.2 seconds to complete all calculation. Maximum error is around 1.5e-12, which is already quite close to machine precision. Error function (analytical values minus approximated values) is shown below. 

So, with just 25 Chebyshev points in our grid, we get remarkably close to the original pricing function values just with the fraction of original processing time. I dare to describe this as nothing less than impressive. This simple example just quickly scratches the surface, since one can imagine the application areas in risk management calculations are endless. Some of these application areas are well described in the book by Ignacio Ruiz and Mariano Zeron.

Let us go through the code. First, we need to import MoCaX library.

import mocaxpy

After this, we create MoCax object. In constructor, the first argument could be a callable function, but in here we just pass None, since the actual values will be delivered to our model on a later stage. The second argument is the number of dimensions. Since in our pricing function we only change the spot price and all the other variables are being fixed, our model is one-dimensional. The third argument is constructed MocaxDomain object, which sets the domain for our model. In the fourth argument, we can set our error threshold, but in here we just pass None. The fifth argument is constructed MocaxNs object, which sets the accuracy for our model.

n_dimension = 1
n_chebyshev_points = 25
domain = mocaxpy.MocaxDomain([[lowest_spot, highest_spot]])
accuracy = mocaxpy.MocaxNs([n_chebyshev_points])
mocax = mocaxpy.Mocax(None, n_dimension, domain, None, accuracy)

At it's current state, our constructed (but still incomplete) Mocax object is only aware of dimensionality, domain and grid of Chebyshev points. The next step is to set the actual function values for each point in our Chebyshev grid. 

chebyshev_points = mocax.get_evaluation_points()
option_values_at_evaluation_points = [call_option(point[0], 100.0, 1.0, 0.25, 0.01) for point in chebyshev_points]
mocax.set_original_function_values(option_values_at_evaluation_points)
call_options_approximations = [mocax.eval([spot], derivativeId=0) for spot in spots]

First, we request back the information on Chebyshev grid points (total of 25 points) from Mocax object. We need to calculate function values for all these points, since we did not give any callable pricing function in our constructor. Next, we ask revaluations from the original pricing function for all Chebyshev grid points (total of 25 revaluations). Then, we handle these calculated values back to Mocax object. Finally, we request call option values from Mocax object. As a response, we will get one million values back from Mocax object: 25 of those values were being explicitly calculated by using original pricing function and the rest are being interpolated by Mocax object.

The complete example program for this section can be found in my GitHub page.

Example: 2-dimensional MoCaX approximation on valuing European call options

Again, let us value a portfolio of European call options, but this time by changing two variables: spot price and time to maturity. Spot price varies between 50 and 150 and time to maturity varies between 0.25 and 2.5. By dividing the both domains into one thousand steps, we get grid containing one million points and for full revaluation we get total of one million calls to our pricing function. As a result, we get the familiar 3D surface.

Processing time is about 330 seconds (5.5 minutes) to complete all calculations. 

Next, we calculate the approximations for option values by using Mocax object and 25 Chebyshev points for the both domains. Revaluation calls are being sent for 625 options and the rest will be interpolated by Mocax object. Processing time is about 13.7 seconds to complete all calculation. So, with just 625 Chebyshev points in our grid, we get remarkably close to the original pricing function values just with the fraction of original processing time. Error surface is shown below (maximum error is around 7.7e-08).

 

The complete example program for this section can be found in my GitHub page.

Further material

The main source for all information is MoCaX Intelligence company page and their collection of high-quality YouTube videos. 

The book "Machine Learning for Risk Calculations: A Practitioner's View" written by Ignacio Ruiz and Mariano Zeron is also available for sale.

QuantsHub YouTube channel is also containing 1-hour video on MoCaX, as presented by Ignacio Ruiz. Concerning the material presented in this post, the most essential videos to get started on using MoCaX as presented by Mariano Zeron can be found in here and in here.

The actual MoCaX library (C++ and Python) can be downloaded from software page. The downloaded package is professionally equipped with MoCaX API documentation (doxygen), as well as with program examples and PDF documentations explaining installation and the main usage of MoCaX API. 

Thanks for reading this blog.

-Mike

QuantLib-Python: Monte Carlo Valuation for Power Reverse Dual-Currency Note (PRDC)

PRDC note

Essentially, PRDC note can be thought as taking a leveraged position on FX forward curve. Floating coupon rate is a function of FX rates, usually defined as follows.

Additional FX term (at t=0) will scale FX rate used in coupon and is usually set to be equal to FX forward at inception. Mathematically, this changes the slope of coupon payoff and therefore, adding desired leverage to the product.

Model for FX stochastics

From pricing perspective, we need to get a view on FX forward curve. This can be simulated by using Garman-KohlHagen stochastic process for an exchange rate, which is available in QuantLib-Python library. In essence, this process is Geometric Brownian Motion, with an added term for foreign interest rate. We will use this model, even it does not necessarily fully reflect the statistical properties of exchange rate.

Applying Ito's rule to Geometric Brownian Motion, gives us the following solution for an exchange rate S observed at time t.

Python program

In order to implement completely new instrument and corresponding new pricing engine in QuantLib, one should first implement these in QuantLib C++ library by following carefully the designed Instrument-Engine scheme and then, interface these newly created C++ types back to Python by using swig. However in here, we will skip the previous approach completely. Instead, we just take "all the useful bits and pieces" available in QuantLib-Python library and use these to implement our own valuation code.

QuantLib-Python object building documentation by David Duarte is available in here.

The complete code is available in here.

Libraries

The following two libraries are being used.

import numpy as np import QuantLib as ql

FX path generator

Path generator implements spot rate simulation, based on a given 1D stochastic process. Generated path is exactly in line with the dates given in coupon schedule. Note, that the exact number of paths to be generated has not been defined in this class. Whoever will be the client for this class, can request as many paths as desired. 

class PathGenerator:
    def __init__(self, valuation_date, coupon_schedule, day_counter, process):        
        self.valuation_date = valuation_date
        self.coupon_schedule = coupon_schedule
        self.day_counter = day_counter
        self.process = process
        #
        self.create_time_grids()
        
    def create_time_grids(self):
        all_coupon_dates = np.array(self.coupon_schedule)
        remaining_coupon_dates = all_coupon_dates[all_coupon_dates > self.valuation_date]
        self.time_grid = np.array([self.day_counter.yearFraction(self.valuation_date, date) 
            for date in remaining_coupon_dates])
        self.time_grid = np.concatenate((np.array([0.0]), self.time_grid))
        self.grid_steps = np.diff(self.time_grid)
        self.n_steps = self.grid_steps.shape[0]
        
    def next_path(self):
        e = np.random.normal(0.0, 1.0, self.n_steps)
        spot = self.process.x0()
        dw = e * self.grid_steps
        path = np.zeros(self.n_steps, dtype=float)
        
        for i in range(self.n_steps):
            dt = self.grid_steps[i]
            t = self.time_grid[i]
            spot = self.process.evolve(t, spot, dt, dw[i])
            path[i] = spot

        return path

Monte Carlo pricer

Monte Carlo pricer class will use path generator for simulating requested number of FX paths. A call for method npv will trigger coupon rates simulation. Simulated rates will then be corrected by taking into account of intro coupon. Finally, a series of cash flows will be created, based on calculated coupon rates and schedule dates.

class MonteCarloPricerPRDC: def __init__(self, valuation_date, coupon_schedule, day_counter, notional, discount_curve_handle, payoff_function, fx_path_generator, n_paths, intro_coupon_schedule, intro_coupon_rate): self.valuation_date = valuation_date self.coupon_schedule = coupon_schedule self.day_counter = day_counter self.notional = notional self.discount_curve_handle = discount_curve_handle self.payoff_function = payoff_function self.fx_path_generator = fx_path_generator self.n_paths = n_paths self.intro_coupon_schedule = intro_coupon_schedule self.intro_coupon_rate = intro_coupon_rate # self.create_coupon_dates() def create_coupon_dates(self): self.all_coupon_dates = np.array(self.coupon_schedule) self.past_coupon_dates = self.all_coupon_dates[self.all_coupon_dates < self.valuation_date] n_past_coupon_dates = self.past_coupon_dates.shape[0] - 1 self.past_coupon_rates = np.full(n_past_coupon_dates, 0.0, dtype=float) self.remaining_coupon_dates = self.all_coupon_dates[self.all_coupon_dates > self.valuation_date] self.time_grid = np.array([self.day_counter.yearFraction(self.valuation_date, date) for date in self.remaining_coupon_dates]) self.grid_steps = np.concatenate((np.array([self.time_grid[0]]), np.diff(self.time_grid))) self.n_steps = self.grid_steps.shape[0] if(self.intro_coupon_schedule==None): self.has_intro_coupon = False else: self.intro_coupon_dates = np.array(self.intro_coupon_schedule) self.remaining_intro_coupon_dates = self.intro_coupon_dates[self.intro_coupon_dates > self.valuation_date] self.n_remaining_intro_coupon_dates = self.remaining_intro_coupon_dates.shape[0] if(self.n_remaining_intro_coupon_dates > 0): self.has_intro_coupon = True else: self.has_intro_coupon = False def simulate_coupon_rates(self): self.simulated_coupon_rates = np.zeros(self.n_steps, dtype=float) for i in range(self.n_paths): path = self.fx_path_generator.next_path() for j in range(self.n_steps): self.simulated_coupon_rates[j] += self.payoff_function(path[j]) self.simulated_coupon_rates = self.simulated_coupon_rates / self.n_paths if(self.has_intro_coupon): self.append_intro_coupon_rates() self.coupon_rates = np.concatenate((self.past_coupon_rates, self.simulated_coupon_rates)) self.n_coupon_cash_flows = self.coupon_rates.shape[0] def append_intro_coupon_rates(self): for i in range(self.n_remaining_intro_coupon_dates): self.simulated_coupon_rates[i] = self.intro_coupon_rate def create_cash_flows(self): self.coupon_cash_flows = np.empty(self.n_coupon_cash_flows, dtype=ql.FixedRateCoupon) for i in range(self.n_coupon_cash_flows): self.coupon_cash_flows[i] = ql.FixedRateCoupon(self.all_coupon_dates[i+1], self.notional, self.coupon_rates[i], self.day_counter, self.all_coupon_dates[i], self.all_coupon_dates[i+1]) self.coupon_leg = ql.Leg(self.coupon_cash_flows) redemption = ql.Redemption(self.notional, self.all_coupon_dates[-1]) self.redemption_leg = ql.Leg(np.array([redemption])) def npv(self): self.simulate_coupon_rates() self.create_cash_flows() self.redemption_leg_npv = ql.CashFlows.npv(self.redemption_leg, self.discount_curve_handle, False) self.coupon_leg_npv = ql.CashFlows.npv(self.coupon_leg, self.discount_curve_handle, False) date = [payment.date() for payment in self.coupon_leg] amount = [payment.amount() for payment in self.coupon_leg] amount[-1] += self.notional pv = [ql.CashFlows.npv(np.array([payment]), self.discount_curve_handle, False) for payment in self.coupon_leg] pv[-1] += self.redemption_leg_npv self.cash_flow_table = np.array([date, amount, pv]) return self.coupon_leg_npv + self.redemption_leg_npv

Hard-coded factory

The creation of Garman-Kohlhagen process is isolated into hard-coded factory method. In here, we just use flat forward curves and non-calibrated short rate models for the sake of simplicity. For more realistic creation of the curves and model calibrations, check out this and this.

def process_factory():
    today = ql.Settings.instance().evaluationDate
    
    domestic_curve = ql.FlatForward(today, ql.QuoteHandle(ql.SimpleQuote(0.01)), ql.Actual360())
    domestic_curve = ql.YieldTermStructureHandle(domestic_curve)
    foreign_curve = ql.FlatForward(today, ql.QuoteHandle(ql.SimpleQuote(0.03)), ql.Actual360())
    foreign_curve = ql.YieldTermStructureHandle(foreign_curve)
    
    fx_vol = ql.QuoteHandle(ql.SimpleQuote(0.1))
    fx_vol_curve = ql.BlackVolTermStructureHandle(ql.BlackConstantVol(today, ql.NullCalendar(), fx_vol, ql.Actual360()))    
    fx_spot = ql.QuoteHandle(ql.SimpleQuote(133.2681))
    
    return ql.GarmanKohlagenProcess(fx_spot, foreign_curve, domestic_curve, fx_vol_curve)

Main program

Main program creates all required components (JPY discount curve, process, schedules, generator and pricer) and then requests npv and cash flow table as Numpy ndarray (indexing: [0] = coupon dates, [1] = cash flows, [2] = present values of cash flows).

today = ql.Date(11, 4, 2023)
ql.Settings.instance().evaluationDate = today

# create JPY discount curve
discount_curve = ql.FlatForward(today, ql.QuoteHandle(ql.SimpleQuote(0.005)), ql.Actual360())
discount_curve_handle = ql.YieldTermStructureHandle(discount_curve)

# create FX process
process = process_factory()

# create schedules for coupon- and intro coupon payments
effectiveDate = ql.Date(3, ql.September, 2015)
terminationDate = ql.Date(3, ql.September, 2041)
coupon_schedule = ql.MakeSchedule(effectiveDate, terminationDate, ql.Period(6, ql.Months), 
    backwards=True, calendar=ql.TARGET(), convention=ql.ModifiedFollowing)

intro_coupon_termination_date = ql.Date(3, ql.September, 2016)
intro_coupon_schedule = ql.MakeSchedule(effectiveDate, intro_coupon_termination_date, ql.Period(6, ql.Months), 
    backwards=True, calendar=ql.TARGET(), convention=ql.ModifiedFollowing)

# create FX path generator
fx_path_generator = PathGenerator(today, coupon_schedule, ql.Actual360(), process)

# create PRDC pricer
notional = 300000000.0
intro_coupon_rate = 0.022
n_paths = 10000
prdc_payoff_function = lambda fx_rate : min(max(0.122 * (fx_rate / 120.0) - 0.1, 0.0), 0.022)
prdc_pricer = MonteCarloPricerPRDC(ql.Settings.instance().evaluationDate, coupon_schedule, ql.Actual360(),notional, 
    discount_curve_handle, prdc_payoff_function, fx_path_generator, n_paths, intro_coupon_schedule, intro_coupon_rate)

# request results
npv_ccy = prdc_pricer.npv()
print('PV in CCY: {}'.format(npv_ccy))
jpy_eur = 145.3275
npv_eur = npv_ccy / jpy_eur
print('PV in EUR: {}'.format(npv_eur))
print()
print('Cash flow dates: {}'.format(prdc_pricer.cash_flow_table[0]))
print()
print('Cash flows: {}'.format(prdc_pricer.cash_flow_table[1]))
print()
print('Present values of cash flows: {}'.format(prdc_pricer.cash_flow_table[2]))

Thanks for reading this blog.

-Mike