Blog: Articles on C++11 and Computational Finance (by Daniel J. Duffy)

Daniel Duffy · 3/14/23

Some random metaprogramming in C++17 and beyond.

void TestTypes()
{ // Using C++17's `std::is_same_v<T, U>` format
// djd

// Compare types and sent the results to the console.
std::cout << std::boolalpha;
std::cout << std::is_same_v<void(int), void(const int)> << ' '; // T
std::cout << std::is_same_v<void(int*), void(const int*)> << ' '; // F
std::cout << std::is_same_v<void(int*), void(int* const)> << ' '; // T

std::cout << '\n' << std::is_same_v<int, const int> << ' '; // F
std::cout << std::is_same_v<int*, const int*> << ' '; // F
std::cout << std::is_same_v<int*, int* const> << ' '; // F

std::cout << '\n' << std::is_same<int, const int>::value << ' '; // F
std::cout << std::is_same<int*, const int*>::value << ' '; // F
std::cout << std::is_same<int*, int* const>::value << ' '; // F
}

Daniel Duffy · 3/26/23

C++ Concepts for the Impatient

C++:

// Test101Concepts.cpp
//
// Simplest example of a system. Context diagram consists only
// of Input and Output systems.
//
// Daniel J. Duffy
//
// We use C++20 Concepts to replace policy-based design. Concepts are
// a game-changer and seem to have been influenced by Haskell's TYPECLASS.
// This signals (among other things) the denise of class-based polymorphism.
//
// https://en.wikipedia.org/wiki/Concepts_(C%2B%2B)
// https://en.wikibooks.org/wiki/Haskell/Classes_and_types
//
// More in ..
//
// Modern Multiparadigm Software Architectures and Design Patterns
// with Examples and Applications in C++, C# and Python Volume I
// Datasim Press 2023, Daniel J.Duffy and Harold Kasperink.
//
// Volume II Ïnteroperability" (
//
//    Call C++ from Python and vice versa
//    Call C# from Pythonand vice versa
//  Call C# from native C++and vice versa
//    (foundations, methods and applications)
//
// In this version the data is a string type. We keep it simple to show
// Concepts essentials. We have applied the method to the design of domain
// architectures. This is now a defined (and repeatable) process.
//
// Tradititonal OOP is not the driver of system development (it is not scalable).
// We base our designs in Domain Architectures (Duffy 2004),  Module Interconnection
// Language (MIL) and associated provides-requires interfaces.
//
// (C) Datasim Education BV 2015-2023
//

#include <string>
#include <iostream>
#include <thread>

// Interface contract specification

// Input (Source system)
template<typename T>
    concept IInput = requires (T x) { x.receive(); };

// Output (Dispatch system)
template<typename T>
    concept IOutput = requires (T x, const std::string& s) { x.send(s); };

// Mediator between Source and Dispatch
template<typename Input, typename Output>
    concept ISUD = IInput<Input> && IOutput<Output>;

template <typename I, typename O> requires ISUD<I,O>
    class SUD
{ // SUD is composed of Input and Output

private:
    I i; O o;
public:
    SUD(const I& input, const O& output) : i(input), o(output) {}
    void run()
    {
        // SUD pulls the data from I, modifies it and then pushes this
        // modified data to O.
        std::string message = std::string("192.1.1.1 ") + i.receive();
        o.send(message);   
    }

    // Callable class, will be called by threads
    void operator ()()
    {
        run();
    }
};

// Instance Systems
class Input
{
public:

    std::string receive() const
    {
        // Get data from hardware device
        return std::string("seq Good morning ..");
    }
};

class Output
{
public:

    void send(const std::string& s) const
    {
        std::cout << "seq send " << s << std::endl;
    }

};

class Input2
{
public:

    std::string receive() const
    {
        // Get data from hardware device. Put in a delay
        std::this_thread::sleep_for(std::chrono::seconds(3));
        return std::string("par 12:34:56:78:9A:BC");
    }
};

class Output2
{
public:

    void send(const std::string& s) const
    {
        std::this_thread::sleep_for(std::chrono::seconds(10));
        std::cout << "par send " << s << std::endl;
    }

};

int main()
{
        // Run programs sequentially
        Input i; Output o;
        SUD<Input, Output> s(i, o);
        s.run();

        Input2 i2; Output2 o2;
        SUD<Input2, Output2> s2(i2, o2);
        s2.run();

        std::cout << "end of sequential\n";

        // Run programs in parallel. We have not placed locks
        // on iosteream!
        std::thread t1(s);
        std::thread t2(s2);

        t1.join(); t2.join();

        return 0;
}

Daniel Duffy · 4/27/23

** The Programming Language Warriors, wise words to ..
(Listen and Learn)

https://lnkd.in/egrzQ297

Modern society is built on software platforms that encompass a great deal of our lives. While this is well known, software is invented by people and this comes at considerable cost. Notably, approximately $331.7 billion are paid, in the U.S. alone, in wages every year for this purpose. Generally, developers in industry use programming languages to create their software, but there exists significant dispersion in the designs of competing language products. In some cases, this dispersion leads to trivial design inconsistencies (e.g., the meaning of the symbol +), while in other cases the approaches are radically different. Studies in the literature show that some of the broader debates, like the classic ones on static vs. dynamic typing or competing syntactic designs, provide consistent and replicable results in regard to their human factors impacts. For example, programmers can generally write correct programs more quickly using static typing than dynamic for reasons that are now known. In this talk, we will discuss three facets of language design dispersion, sometimes colloquially referred to as the “programming language wars.” First, we will flesh out the broader impacts inventing software has on society, including its cost to industry, education, and government. Second, recent evidence has shown that even research scholars are not gathering replicable and reliable data on the problem. Finally, we will give an overview of the facts now known about competing alternatives (e.g., types, syntax, compiler error design, lambdas).

Daniel Duffy · 4/27/23

Like a breath of fresh air

The Programming Language Wars - Microsoft Research

Modern society is built on software platforms that encompass a great deal of our lives. While this is well known, software is invented by people and this comes at considerable cost. Notably, approximately $331.7 billion are paid, in the U.S. alone, in wages every year for this purpose...

www.microsoft.com

Daniel Duffy · 5/16/23

Gradient Descent Method in continuous space
By solving an ODE

C++:

// TestExp5.cpp
//
// dx/dt = - grad(U(x)) , grad == gradient
// Transform f(x) = 0 to a least-squares problem, U = f*f
//
// Tip iceberg: minimise a function using an ODE solver for a gradient system.
//
// compute z = exp(5); f(z) = z - exp(5); argmin f(z)*f(z)
//
// It's a sledgehammer but the method can be applied to constrained optimisation for ODE systems.
// To be discussed elsewhere. It can be applied in all sorts of places.
// Current ML wisdom uses Gradient Descent (GD) method, which is Euler's method, and much weaker in many
// ways to the approach taken here. It's much more robust and general than GD.
//
// Exercise:
// 1. generalise to ODE systems.
// 2. generalise to constrained optimisation and penalty method.
// 3. methods to compute gradient in n dimensions.
// 4. Use in an ANN instead of flaky GD.
//
// The world is continuous, but the mind is discrete. David Mumford
// Natura non facit saltum (Nature does nothing in jumps)
//
// My take is that Gradient Descent is a fabrication; ODE gradient system are closer
// to the physical world.
//
// DJD
//
// Using ODEs for optimisation.
//
// http://www.unige.ch/~hairer/preprints/gradientflow.pdf
// https://authors.library.caltech.edu/26703/2/postscript.pdf
// https://blogs.mathworks.com/cleve/2013/10/14/complex-step-differentiation/
// https://www.dam.brown.edu/people/geman/Homepage/Image%20processing,%20image%20analysis,%20Markov%20random%20fields,%20and%20MCMC/Diffusions%20and%20optimization.pdf
// https://www.cs.ubc.ca/~ascher/papers/adhs.pdf
//
// Nesterov's method
// https://epubs.siam.org/doi/epdf/10.1137/21M1390037
//
//
// For training courses, see
//
// https://www.datasim.nl/onlinecourses/97/distance-learning-ordinary-and-partial-differential-equations
//
// https://www.datasim.nl/application/files/9015/4809/1157/DL_Ordinary_and_Partial_Differential_Equations.pdf
//
// https://www.datasim.nl/onlinecourses
//
// See also my book on ODE/PDE/FDM in computational finance
//
// https://www.wiley.com/en-us/Numerical+Methods+in+Computational+Finance:+A+Partial+Differential+Equation+(PDE+FDM)+Approach-p-9781119719670
//
// (C) Datasim Education BV 2018-2023
//
//


#include <iostream>
#include <vector>
#include <cmath>
#include <complex>


// https://www.boost.org/doc/libs/1_82_0/libs/numeric/odeint/doc/html/index.html
#include <boost/numeric/odeint.hpp>

// Preliminary notation
using value_type = double;
using state_type = std::vector<value_type>;

// Using C++11 functions
template <typename T>
    using FunctionType = std::function<T(const T& arg)>;
using CFunctionType = FunctionType<std::complex<value_type>>;


template <typename T>
    FunctionType<T> operator * (FunctionType<T>& f, FunctionType<T>& g)
{ // Multiplication (higher-order function)

    return [=](T x)
    {
        return  f(x)*g(x);
    };

}

// Complex Step Method for gradients (others: divided difference, AAD, analytic,..)
value_type CSM(const CFunctionType& f, value_type x, value_type h)

{ // df/dx at x using the Complex step method (Trapp/Squire)

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

std::complex<double> ObjFunc(const std::complex<value_type>& z)
{ // My hard-coded specific function

    // compute z = exp(5); f(z) = z - exp(5); argmin f(z)*f(z)
    const std::complex<double> a(std::exp(5.0), 0.0); // => 148.413
   
    // Original function f(x) = 0 (Newton)
    CFunctionType f = [&](const std::complex<double>& z) { return z - a; };
   
    CFunctionType f2 = f * f;

    // Least squares function
    return f2(z);
}

// Free function to model RHS in dy/dt = RHS(t,y)
void Ode(const state_type &x, state_type &dxdt, const value_type t)
{
    //dxdt[0] = 2.0 * (x[0] - std::exp(5.0)); // Exact derivative
    const double h = 1.0e-156; // Small h but no overflow!
    dxdt[0] = CSM(ObjFunc, x[0], h);

    // Transform semi-infinite time interval (0, infinity) to (0,1)
    // Then we look at asymptotic behaviour..
    double tau = (1.0 - t) * (1.0 - t);
    dxdt[0] = -dxdt[0]/tau;
}

void print(std::string name, std::size_t steps, value_type v)
{
    std::cout << "Number of steps " << name << std::setprecision(16) << steps << "approximate: " << v << '\n';
}

int main()
{
    namespace Bode = boost::numeric::odeint;

   
    // Initial condition
    value_type L = 0.0;
    value_type T = 0.99; // HEURISTIC simulates T = infinity
    value_type initialCondition = 1110.548; // Convergence independent of IC?
    value_type dt = 1.0e-5;

    {
        // Cash Karp, middle-of-road
        state_type x{ initialCondition };
        std::size_t steps = Bode::integrate(Ode, x, L, T, dt);
        print("Cash Karp ", steps, x[0]);
    }
    {
        // Bulirsch-Stoer method, high-class solver
        state_type x{ initialCondition};
        Bode::bulirsch_stoer<state_type, value_type> bsStepper;        // O(variable), controlled stepper
        std::size_t steps = Bode::integrate_const(bsStepper, Ode, x, L, T, dt);
        print("Bulirsch-Stoer ", steps, x[0]);

    }

    return 0;
}

Daniel Duffy · 8/25/23

An Explicit Implied Volatility Formula
International Journal of Theoretical and Applied Finance, Vol. 20, no. 7, 2017

Dan Stefanica
Baruch College, City University of New York

Rados Radoicic
CUNY Baruch College

C++:

// DS_RR.cpp
//
/* An Explicit Implied Volatility Formula
International Journal of Theoreticaland Applied Finance, Vol. 20, no. 7, 2017


Dan Stefanica
Baruch College, City University of New York

Rados Radoicic
CUNY Baruch College

Date Written : January 30, 2017
*/
//
// https://papers.ssrn.com/sol3/papers.cfm?abstract_id=2908494
//
// Programmed in C++ by Daniel J.Duffy 2023
// Datasim Education BV 2023
//



#include <cmath>
#include <iostream>

using value_type = double;

// Call parameters
/*
const value_type S = 2080;
const value_type K = 3050;
const value_type r = 0.02;
const value_type b = 0;
const value_type sig = 0.1097904;
const value_type T = 0.5;

*/

const value_type S = 100;
const value_type K = 15000;
const value_type r = 0.08;
const value_type b = r;
const value_type sig = 0.0004;
const value_type T = 1000.0;


// ATM + small vol OK
// OTM, ITM + small vol not OK
// r == vol == 0 is essentially a degenerate case??
/*
const value_type S = 60;
const value_type K = 65;
const value_type r = 0.08;
const value_type b = r;
const value_type sig = 0.302; // small vol has issues < 0.1 !!
const value_type T = 0.25;
*/

/*
const value_type S = 100;
const value_type K = 6.30957;
const value_type r = 0.0;
const value_type b = r;
const value_type sig = 0.3409436;
const value_type T = 1.0;
*/
// Polya's 1949 approximation to N(x)
value_type G(value_type x)
{
    return 1 / (exp(-358 * x / 23 + 111 * atan(37 * x / 294) + 1));
    value_type y = std::sqrt(1.0 - std::exp(-2.0 * x * x / 3.141592653589793));
    if (x >= 0.0)
        return 0.5 * (1.0 + y);
    else
        return 0.5 * (1.0 - y);
}

value_type N(value_type x)
{ // The approximation to the cumulative normal distribution

    //return cndN(x);
    return G(x);
}

value_type OptionPrice(value_type x)
{ // Call price, x == sig

    value_type tmp = x * std::sqrt(T);

    value_type d1 = (std::log(S / K) + (r + (x * x) * 0.5) * T) / tmp;
    value_type d2 = d1 - tmp;
    return (S * std::exp((b - r) * T) * N(d1)) - (K * std::exp(-r * T) * N(d2));
}

value_type Cm = OptionPrice(sig); // Market price
// Stefanica and Radoicic volatiliy estimator
value_type y = std::log(S / K);
value_type alpha = Cm / (K * std::exp(-r * T));
value_type R = 2 * alpha - std::exp(y) + 1.0;

value_type pi = 3.141592653589793;
value_type t1 = std::exp((1.0 - 2 / pi) * y);
value_type t2 = std::exp((-1.0 + 2 / pi) * y);
value_type A = std::pow(t1 - t2, 2);
value_type B1 = 4.0 * (std::exp(-y) + std::exp(-2 * y / pi));
value_type B2 = -2.0 * (std::exp(-y) * (t1 + t2) * (std::exp(2 * y) + 1.0 - R * R));
value_type B = B1 + B2;
value_type C = std::exp(-2.0 * y)* (R * R - std::pow(std::exp(y) + 1.0, 2) - R * R);

value_type beta = 2.0 * C / (B + std::sqrt(B * B + 4.0 * A*C));
value_type gamma = -0.5 * pi * std::log(beta);

value_type compute()
{
    if (y >= 0.0)
    {
        value_type C0 = K * std::exp(-r * T) * (std::exp(y) * A * (std::sqrt(2 * y) - 0.5));
        value_type sig;
        if (Cm <= C0)
        {
            sig = (std::sqrt(gamma + y) - std::sqrt(gamma - y)) / std::sqrt(T);
        }
        else
        {
            sig = (std::sqrt(gamma + y) + std::sqrt(gamma - y)) / std::sqrt(T);
        }
    }
    else
    {
        value_type C0 = K * std::exp(-r * T) * (std::exp(y)/2 - A * (std::sqrt(2 * y) - 0.5));
        value_type sig;
        if (Cm <= C0)
        {
            sig = (-std::sqrt(gamma + y) + std::sqrt(gamma - y)) / std::sqrt(T);
        }
        else
        {
            sig = (std::sqrt(gamma + y) - std::sqrt(gamma - y)) / std::sqrt(T);
        }
    }

    return sig;

}


int main()
{

    std::cout << "DS/RR: " << compute() << '\n';
}

Paul Lopez · 8/31/23

This is great!

Daniel Duffy · 8/31/23

Paul Lopez said:
This is great!

Thanks.
Many student questions on QN are addressed here, which is the reason I posted it in the first place

seek and ye shall find

Matthew 7:7

Daniel Duffy · 11/20/23

"People talk about teaching programming, I think that's a great idea. The problem is, as near as I can tell, they're *not* teaching programming ... they're teaching coding.

Hardly anybody in the outside world understands the difference between #programming and #coding "

— Dr. Leslie Lamport (Turing Award 2013)

Daniel Duffy · 11/21/23

The Microsoft Implementation of CRTP in Active Template Library (ATL) was independently discovered, also in 1995, by Jan Falkin, who accidentally derived a base class from a derived class. Christian Beaumont first saw Jan's code and initially thought it could not possibly compile in the Microsoft compiler available at the time. Following the revelation that it did indeed work, Christian based the entire ATL and Windows Template Library (WTL) design on this mistake.

Daniel Duffy · 12/12/23

Daniel Duffy · 2/12/24

C++:

The Black Scholes option pricer in Fortran95
// DJD
real*8 Function pdf(x)

            real*8 x,A
            A = 1.0/Sqrt(2.0*3.1415)

            pdf = A * Exp(-0.5*x*x)

            return
        end

real*8 Function cdf(x)

                real*8 DPI,x,L,k,a1,a2,a3,tmp,pdf

                a1 = 0.4361836
                a2 = -0.1201676
                a3 = 0.9372980
                DPI = 3.1415926535897932

                L = Abs(x)
                k = 1. / (1. + 0.33267*x)

                tmp = a1*k+ a2 * k**2. + a3 * k**3.

                cdf =  1.0 - pdf(x)*tmp
                if(x.lt.0.) then

                    cdf = pdf(x)*tmp

                end if

                return

            end

real*8 Function BlackScholes(S, X, T, r, v)
      ! Put
        real*8 S,X,T,r,v

        real*8 d1, d2
        d1 = (Log(S / X) + (r + v**2. / 2.) * T) / (v * Sqrt(T))
        d2 = d1 - v * Sqrt(T)

        BlackScholes = X * Exp(-r * T) * cdf(-d2) - S * cdf(-d1)

        Return

        End

! A fortran95 program for BS option
! By djd
!
program main
  implicit none
  integer anyKey

  real*8 S,K,T,r,v
  real*8 BlackScholes
  real*8 price

  S = 60.0
  K = 65.0
  T = 0.25
  r = 0.08
  v = 0.3

  price  = BlackScholes(S,K,T,r,v)
  write(*,*) price

  anyKey = system("pause")

end

Paul Lopez · 2/12/24

Daniel Duffy said:

C++:

The Black Scholes option pricer in Fortran95
// DJD
real*8 Function pdf(x)

            real*8 x,A
            A = 1.0/Sqrt(2.0*3.1415)

            pdf = A * Exp(-0.5*x*x)

            return
        end

real*8 Function cdf(x)

                real*8 DPI,x,L,k,a1,a2,a3,tmp,pdf

                a1 = 0.4361836
                a2 = -0.1201676
                a3 = 0.9372980
                DPI = 3.1415926535897932

                L = Abs(x)
                k = 1. / (1. + 0.33267*x)

                tmp = a1*k+ a2 * k**2. + a3 * k**3.

                cdf =  1.0 - pdf(x)*tmp
                if(x.lt.0.) then

                    cdf = pdf(x)*tmp

                end if

                return

            end

real*8 Function BlackScholes(S, X, T, r, v)
      ! Put
        real*8 S,X,T,r,v

        real*8 d1, d2
        d1 = (Log(S / X) + (r + v**2. / 2.) * T) / (v * Sqrt(T))
        d2 = d1 - v * Sqrt(T)

        BlackScholes = X * Exp(-r * T) * cdf(-d2) - S * cdf(-d1)

        Return

        End

! A fortran95 program for BS option
! By djd
!
program main
  implicit none
  integer anyKey

  real*8 S,K,T,r,v
  real*8 BlackScholes
  real*8 price

  S = 60.0
  K = 65.0
  T = 0.25
  r = 0.08
  v = 0.3

  price  = BlackScholes(S,K,T,r,v)
  write(*,*) price

  anyKey = system("pause")

end

love it!

okay · 2/12/24

Hello, are there any other C++ courses available for developers besides the one @Duffy , currently offered?

Paul Lopez · 2/12/24

okay said:
Hello, are there any other C++ courses available for developers besides the one @Duffy , currently offered?

You mean at a university or just online?

MikeLawrence · 2/12/24

okay said:
Hello, are there any other C++ courses available for developers besides the one @Duffy , currently offered?

Dr. Duffy currently offers/has put together many courses.

Two of them are on here, TA'd by under the 'online courses' tab.

Online Courses

quantnet.com

Many others are always open through his educational company, Datasim:

Online Courses :: Datasim

www.datasim.nl

Daniel Duffy · 2/12/24

okay said:
Hello, are there any other C++ courses available for developers besides the one @Duffy , currently offered?

What kinds of C++ courses? The internet is (used to be) full of C++ courses.

// btw wrong thread; this is about articles, not courses.

@okay

Daniel Duffy · 3/8/24

An emerging design pattern, an improved 3-layer model for OO inheritance.

This exposes a number of flaws in traditional approaches to OOP, plus a multiparadigm resolution thereof.

You could try designing it using C++20 Concepts.

@APalley
@Paul Lopez
@Andy Nguyen

Daniel Duffy · 4/26/24

A general related remark:

Remark on code speedup.

The prevailing "mindset" is to write sequential (single-threaded) code and to profile it to locate performance bottlenecks. Then we pepper the code using #pragma and @decorator etc. This is fine-grained paralled pattern (e.g. parallel loop) and not optimal in the amount of time to profile the code. Lots of iterations and unpredictable benefits.

Plan B is to look for potential parallelism/concurrency and to create a task graph that we eventually implement on multi-core machines. The we can exploit the hardware to get better speed-up. In this case we speak of coarse-grained parallelism that is implemented by parallel design patterns.

Moving to prallel design demands a change of mindsef on howe we "look" at algorithms. And skills update.

I discuss these topics in Duffy C++ book (2018) as well as this article.

Error - Cookies Turned Off

onlinelibrary.wiley.com

Daniel Duffy · 6/26/24

Partial differential equations meet Brownian motion (Thanks to Albert Einstein).

Fokker–Planck equation - Wikipedia

en.wikipedia.org

Blog: Articles on C++11 and Computational Finance (by Daniel J. Duffy)

C++ author, trainer

C++ author, trainer

C++ author, trainer

C++ author, trainer

C++ author, trainer

C++ author, trainer

C++ author, trainer

seek and ye shall find​

C++ author, trainer

C++ author, trainer

C++ author, trainer

C++ author, trainer

C++ author, trainer

C++ author, trainer

Attachments

C++ author, trainer

C++ author, trainer

Attachments

Similar threads

seek and ye shall find