Computing sensitivities/greeks in computational finance, optimization and Machine Learning (ML): the Candidates

Daniel Duffy · 5/8/19

How many Ways are there to compute Derivatives of a Function?

Differentiation, Sensitivities (greeks), gradients, Jacobians etc.

The computation of derivatives and gradients of scalar and vector-valued functions is a pervasive and ongoing activity in many areas of computational finance, optimisation, engineering and Machine Learning (ML), to name just a few. In general, we may need to compute derivatives as part of some algorithm or software program. Some specific applications are:

. Computing the Black Scholes greeks/sensitivities.

. Hedge ratios in fixed income and credit risk applications (for example, vega, duration and complexity).

. Optimisation of design shape parameters.

. Computing derivatives in Neural Network (NN) algorithms, for example (Stochastic) Gradient Descent Method (SGD).

. Hedge ratios for stochastic differential equations (SDE), using Automatic Differentiation (AD), for example.

. Global optimisation of (Langevin) SDEs and gradient ODE systems.

. Applications (many) in which gradients, Jacobians and Hessians need to be calculated.

. Robust calibration and volatility surfaces.

There are several established methods to address the problem of computing the derivatives of functions:

. Analytically (using calculus mainly).

. Using finite differences (or to be more precise, divided differences) to approximate derivatives (this is a well-known approach in PDE models).

. Automatic (Algorithmic) Differentiation (AD).

. Continuous Sensitivity Equation (CSE). In this case we view the derivative of a function (with respect to a parameter) as a quantity that satisfies an ordinary or partial differential equation with auxiliary initial and boundary conditions.

. Approximating the function whose derivatives we wish to compute by an ‘easier’ function, for example a cubic spline whose derivatives are easily computed.

. The Complex Step Method (CSM). This little known but robust method to compute derivatives using function values alone (and no need to take divided differences). The method employs complex numbers and complex arithmetic.

There is no golden rule as to which of the above methods to use in general because the choice depends on the requirements and context. In general, we wish a given method to be suitable for the problem at hand and ideally it should be efficient, robust and the resulting code should be easy to use or to create. Some questions to address are: exact or approximate values of the derivative and do we wish to have a discrete or continuous approximation.

The above methods are readily supported in C++, C#, Python and several (open-source) libraries.

In this note we give two examples of using The Complex Step Method in C++ and Python. The first example is generic to show how the method works and the second example shows how to compute option ‘greeks’ in derivative pricing.

The CSM is discussed in the links attached to my source code. I give two C++ programs and a Python proof-of-concept code to make the method more accessible. In this case we used and modified the C++ code to get a working version in Python up and running.

The rationale for the code is to show how CSM works. Extensions to more complex cases are possible and will be discussed elsewhere.

Daniel Duffy · 5/8/19

The forum doesn't accept .py format file, so I do an insert code

Code:

# TestComplexStepMethod.py
#
# The complex step method for differentiation.
#
# This code was ported from C++
#
# https://pdfs.semanticscholar.org/3de7/e8ae217a4214507b9abdac66503f057aaae9.pdf
#
# http://mdolab.engin.umich.edu/sites/default/files/Martins2003CSD.pdf
#
# (C) Datasim Education BV 2019
#

import numpy, math, cmath, random

j = cmath.sqrt(-1)

def Derivative(f, x, h):
# df/dx at x using tbe Complex Step Method

    z = x + h*j
    fz = f(z)
    return fz.imag/h;

def DerivativeZ(f, z):
# df/dx at z using tbe Complex Step Method

    #z = x + h*j
    fz = f(z)
    return fz.imag/h;

def DerivativeComposite(f, g, x, h):
# df(g)/dx at x  = df/dg X dg/dx using the Complex Step Method
# i.e. function composition

    z = x + h*j
    dgz = DerivativeZ(g,z)
   
    w = g(z)
    dfw = f(w)
 
    return (dfw.imag)/h + dgz.imag/h;

def SquireTrapp(t):
    n1 = cmath.exp(t);
    d1 = cmath.sin(t);
    d2 = cmath.cos(t);

    return n1 / (d1*d1*d1 + d2*d2*d2);

def func(t):
    return cmath.exp(t);

def func2(t):
    return cmath.cos(t);

def f(t):
    return cmath.exp(t);

def g(t):
    return cmath.exp(t);

x = 1.5
h = 0.001

ans = Derivative(SquireTrapp, x, h)
print(ans)

x = 5.0
ans = Derivative(func, x, h)
print(ans)

x = 1# math.sqrt(5)
ans2 = DerivativeComposite(g, g, x, h)
print(ans2)

b = math.exp(1)* math.exp(math.exp(1))
print(b)

Daniel Duffy · 5/8/19

Computing Gradients and Derivatives of Functions in Finance, Optimisation and Machine Learning.

Computing Gradients and Derivatives of Functions in Finance, Optimisation and Machine Learning. How many Ways are there to compute Derivatives of a Function? Differentiation, Sensitivities (greeks), gradients, Jacobians etc.

www.linkedin.com

Daniel Duffy · 5/8/19

Code:

// TestComplexStep.cpp
//
// Complex-step method to compute approximate derivatives.
// Example is scalar-valued function of a scalar argument.
//
// https://pdfs.semanticscholar.org/3de7/e8ae217a4214507b9abdac66503f057aaae9.pdf
//
// http://mdolab.engin.umich.edu/sites/default/files/Martins2003CSD.pdf
//
// (C) Datasim Education BV 20199
// dduffy@datasim.nl
//

#include <functional>
#include <complex>
#include <iostream>
#include <iomanip>
#include <cmath>

// Notation and function spaces
using value_type = double;

template <typename T>
    using FunctionType = std::function < T(const T& c)>;
using CFunctionType = FunctionType<std::complex<value_type>>;


// Test case from Squire&Trapp 1998
template <typename T> T SquireTrapp(const T& t)
{
    T n1 = std::exp(t);
    T d1 = std::sin(t);
    T d2 = std::cos(t);

    return n1 / (d1*d1*d1 + d2*d2*d2);
}

template <typename T> T func2(const T& t)
{ // Derivative of e^t, sanity check

    return std::exp(std::pow(t,1));
//    return std::exp(std::pow(t, 5));

}
    
value_type Derivative(const CFunctionType& f, value_type x, value_type h)
{ // df/dx at x using tbe Complex step method

    std::complex<value_type> z(x, h); // x + ih, i = sqrt(-1)
    return std::imag(f(z)) / h;
}

int main()
{
    // Squire Trapp
    double x = 1.5;    double h = 0.1;
    do
    {
        std::cout << std::setprecision(12) << Derivative(SquireTrapp<std::complex<value_type>>, x, h) << '\n';
        h *= 0.1;

    } while (h > 1.0e-300);

    // Exponential function (101 sanity check)
    x = 5.0;
    h = 1.0e-10;
    std::cout << "Exponential 1: " << std::setprecision(12) << Derivative(func2<std::complex<value_type>>, x, h) << '\n';

    return 0;
}

Daniel Duffy · 5/8/19

Code:

// TestGreeks101.cpp
//
// Test code.
//
// Computing theta and vega (as examples) using the Complex Step Method and
// compare to exact solution.
// We use the Faddeeva (Kramp) function.
//
// Use call options as test case. We have separate functions for double and complex
// cases. They could be 'merged' into one template function, later. It is an optimisation step.
//
// The method can be generalised to multidimensioanl space.
//
// (C) Datasim Education BV 2018-2019
//

#include "Faddeeva.hh"
#include <complex>
#include <iostream>
#include <cmath>
#include <iomanip>


double K = 65.0;
//double T = 0.25;
double r = 0.08;
double b = 0.05;
//double v = 0.3;

double S = 55.0;

template <typename T>
    T n(T x)
{ // pdf

    const T A = 1.0 / std::sqrt(2.0 * 3.14159265358979323846);
    return A * std::exp(-x*x*0.5);
}

double N(double x)
{ // cdf

    return std::erfc(-x / std::sqrt(2.0))/2.0;
}

std::complex<double> N(std::complex<double>& z)
{
    double relErr = 1.0e-16;
    return Faddeeva::erfc(-z / std::sqrt(2.0), relErr) / 2.0;
}

double BS(double T, double v)
{ // Call option

    double LOG = std::log(S/K);
    double den = v*std::sqrt(T);
    
    double d1 = (LOG + (b + v*v / 2)*T) / den;
    double d2 = d1 - den;

    double result = S * N(d1)*std::exp((b-r)*T) - K*std::exp(-r*T)* N(d2);

    return result;

}


std::complex<double> BS(std::complex<double>& T, std::complex<double>& v)
{
    std::complex<double> a = S;
    std::complex<double> LOG = std::log(S / K);
    std::complex<double> den = v*std::sqrt(T);

    std::complex<double> d1 = (LOG + (b + v*v * 0.5)*T) / den;
    std::complex<double> d2 = d1 - den;

    std::complex<double> result = S * N(d1)*std::exp((b - r)*T) - K*std::exp(-r*T)* N(d2);

    return result;
}

double BSTheta(double T, double v)
{
    double den = v * std::sqrt(T);
    double d1 = (std::log(S / K) + (b + (v*v)*0.5) * T) / den;
    double d2 = d1 - den;

    double t1 = (S* std::exp((b - r)*T) * n(d1) * v * 0.5) / std::sqrt(T);
    double t2 = (b - r)*(S * std::exp((b - r)*T) * N(d1));
    double t3 = r * K * std::exp(-r * T) * N(d2);

    return  -(t1 + t2 + t3);

}




std::complex<double> BSTheta(std::complex<double> T, std::complex<double>& v)
{
    
    double h = 1.0e-43;
    std::complex<double> z1(std::real(T), h );
    std::complex<double> z2(std::real(v), 0);
    std::complex<double> result = BS(z1,z2);

    return -std::imag(result) / h;
}


double BSVega(double T, double v)
{
    double den = v * std::sqrt(T);
    double d1 = (std::log(S / K) + (b + (v*v)*0.5) * T) / den;

    return S * std::exp((b - r)*T) * n(d1) * std::sqrt(T);
}

std::complex<double> BSvega(std::complex<double> T, std::complex<double>& v)
{
    
    double h = 1.0e-16; // hard-coded
    std::complex<double> z1(std::real(T), 0);
    std::complex<double> z2(std::real(v), h);
    std::complex<double> result = BS(z1, z2);

    return std::imag(result) / h;
}

int main()
{
    double relErr = 1.0e-16;
    std::complex<double> z(1.0, 0.01); auto v = Faddeeva::erfc(z, relErr);
    std::cout << v << '\n';

    std::complex<double> TC(0.25, 0.0); std::complex<double> vC(0.3, 0.0);
    std::cout << std::setprecision(16) << "CS " << BS(TC, vC) << ", " << BSTheta(TC,vC) << ", " << BSvega(TC, vC) <<  '\n';

    double TR(0.25); double vR(0.3);
    std::cout << "Exact " << std::setprecision(16)<< BS(TR,vR) << ", " << BSTheta(TR,vR) << ", " <<  BSVega(TR, vR) << '\n';

    return 0;
}

Daniel Duffy · 5/9/19

The nice thing about the Complex Step Method is that it works for all values of h without catastrophic subtraction cancellation errors as with traditional divided differences. You can even take h=1.0e−303 but accuracy has an "ergodic" value around h=1.0e−6 (roughly speaking). If the sensitivity satisfies ODE/SDE/PDE then the whole process is even more robust. In that case we can ensure no-arbitrage by using heavy machinery of maximum principle for PDEs and using monotonic finite difference schemes, for example PDEs for vega, delta and gamma. This is essentially the CSE method.

ExSan · 5/10/19

Daniel Duffy said:
Computing Gradients and Derivatives of Functions in Finance, Optimisation and Machine Learning.

Computing Gradients and Derivatives of Functions in Finance, Optimisation and Machine Learning. How many Ways are there to compute Derivatives of a Function? Differentiation, Sensitivities (greeks), gradients, Jacobians etc.

www.linkedin.com

fantastic!!! Thank YOU !
R

Computing sensitivities/greeks in computational finance, optimization and Machine Learning (ML): the Candidates

Daniel Duffy

C++ author, trainer

Attachments

Daniel Duffy

C++ author, trainer

Daniel Duffy

C++ author, trainer

Computing Gradients and Derivatives of Functions in Finance, Optimisation and Machine Learning.

Daniel Duffy

C++ author, trainer

Daniel Duffy

C++ author, trainer

Daniel Duffy

C++ author, trainer

ExSan

Computing Gradients and Derivatives of Functions in Finance, Optimisation and Machine Learning.

Similar threads