Black-Scholes-Merton, F# and AVX2 – #1

In this series I am presenting another example of how to use F# in the financial industry.

Black-Scholes-Merton (from now simply Black-Scholes or BS) is a well known name in the financial industry. It identifies a mathematical model for options pricing, and derives from the surnames of the two eminent economists who primarily proposed the model.

F# is a fully fledged general purpose functional programming language for the .NET universe.

AVX2 (Advanced Vector Extensions 2) is the name given to one of the latest extensions to the x86 processors family instruction set.

We will primarily see how to implement the calculation of option prices according to the Back-Scholes model using exclusively F#. Later, we will see how to push the performance of our calculation using first C++ and then unveiling the high level compiler layer, going down to implement the calculation using plain assembly language and optimizing manually the calculation to leverage on the AVX2 instruction set, which gives us a series of advanced CPU built-in vector instructions. In the process we will show how to call an assembly procedure directly from F#, without passing through a C/C++ calling convention layer.

Option contract

So what is exactly an option in the financial market? It is a contract which gives the buyer (the owner or holder of the option) the right, but not the obligation, to buy or sell an underlying asset or instrument at a specified strike price on a specified date (Wikipedia).

So they are financial products. They are not primary assets, though, since there is an underlying asset that defines them. Hence they belong to the category of derivatives.

For example, on the 1st of February we buy an option that gives us the right to buy 100 shares of ACME Inc. at the price of 142 on May 5th (in the future). We pay 90 for that option contract. We then wait for the 5th of May. If that day ACME Inc. is quoted more than 142, let’s say 147, we exercise our right, buy 100 shares paying 14200 and sell them at market for 14700, realizing a net profit of 500 – 90 = 410. If instead the quotation is less than 142, for instance 139, we do not exercise our right, because we would pay 14200 for something worth 13900. Our operation ends with a loss, the cost of the option contract (90).

Options are used in several ways: the example above shows how to use them as a speculative investment, but this is not the only way. Often options are used to mitigate risks associated to different operations. For example, we manage a fund which has ACME Inc. shares in portfolio. We believe the share has reached a level of maturity so we plan to discharge it from our portfolio, but not just yet, because in a few weeks we expect them to pay dividends and historically ACME Inc. dividends are quite good. Our risk analysts, though, have advised that due to recent changes in worldwide policies, the markets where ACME Inc. operates could be subject to turbulence in the next months, which increases the likelihood that ACME Inc. share sustain substantial price drop. The projections tell that the drop could be so severe that our fund could suffer losses greater than the profits we expect from the dividends. Our fund is a conservative one, not aggressive, so we don’t  want to expose it to those risks. Today ACME Inc. is priced 87. We then decide to sacrifice some of the dividends expected profit and protect our fund by buying an option. In this case, we want the right to sell our shares at 84 in 3 months. If at expiration ACME Inc. has dropped to less than 84, e.g 77, we will exercise out right and discharge our portfolio selling ACME Inc. at 84, with an advantage of 7 (84-77) over selling them at market.

Types of option contracts

There are many type of options available nowadays, but the most important terminology can be summarized as follows

  • Call/Put: a Call option gives the right to buy the underlying asset, a Put option gives right to sell it
  • Expiration: the time when the right associated with the option can be exercised. Linked to the expiration are the concept of maturity. Also, we indicate with “time to maturity” the amount of time missing to maturity
  • Strike price: the price we have the right to buy or sell the underlying asset at, when mature.

A distinction exists between so called European and American option kinds. American options can be exercised even before maturity, European can’t. This publication refers to European options.

The pricing problem

The problem with such a contract is: how do we price it? What is the fair price for an option? We need criteria that allow us to give a concrete price to options contracts, starting from the characteristics of the underlying asset. The merit of the Black-Scholes (BS) model is exactly this. BS models the financial market using statistical math and after some proper elaboration, obtains mathematical formulas that, based on observable market parameters give a value for an option in respect to reasonable assumptions.

As software engineers our scope if to take these formulas and write computer programs that implement them. This can be done even without a full understanding of how they are obtained through mathematical and statistical theoretical reasoning. For those (like me) that like math and prefer to have an understanding of the theoretical grounds they are based on, the following block tries to give a simplified explanation of the most important concepts and some references. The reader is welcome to jump over the block and continue with the implementation notes.

To understand this block the reader should have some familiarity with basics of calculus and math statistics.

Analytical vs. Statistical, Deterministic vs. Stochastic
From physics, we know that a particle moving with velocity v and acceleration a changes its position p in time by the analytical differential equation

dp = v*dt+a*t*dt

where dt is a sufficiently small time interval.
This motion law is deterministic, because all the elements in it are precisely known at any given time.
Financial assets, though, have prices that change in time with unknown laws. Intuitively we can say that there should be an analytical formula that describes price motions, but it’s overly difficult to find it and, even if we did, it would probably be too complex to be effectively usable.
This is where stats come in. We give up our aim to precisely (analytically, deterministically) calculate the price, and introduce the concept of probability. We calculate the probability for the price to assume a certain value, or value range.
Price movements, we say, are random, but still determined by an environment with known characteristics (market); so even if we don’t know exactly how to calculate the price, we can say (for example) that the probability that tomorrow’s ACME Inc. price is 100 times today’s price is much lower than it is 1% higher. Statistical process modeling aims to mathematically evaluate these probabilities.

With this premises, we’ll now follow a bottom-up path to see step-by-step how, starting from a simple random concept as coin toss, we arrive to the BS formulation for the fair price of an option.

Coin Toss
A coin toss is the typical binary random event. It makes sense to associate to the two possible outcomes with the numbers -1 and 1. Just this gives us the ability to define what in stats is the Expected Value or Expectation (what value can we expect from a coin toss event). If we name Z the coin toss event, the outcome can be -1 or 1 so the expectation is 0, in formula


Let’s now think to a series of consecutive coin tosses, named Zi; from the definition of Z we get

E(Z_i*Z_j) = 0
E(Z_i^2) = 1

we define Sn as the sum of the outcomes of n consecutive coin tosses

S_n = \displaystyle\sum^n_{i=1} Z_i

We can then calculate the following expectations

E(S_n) = E(\displaystyle\sum^n_{i=1} Z_i) = \displaystyle\sum^n_{i=1} E(Z_i) = 0
E(S_n^2) = E(Z_1^2 + 2Z_1*Z_2 + ...) = n*E(Z^2) = n

Therefore, the sum of the outcomes of n consecutive coin tosses has an expected value of 0 and a quadratic variation (variance) of n.

From discrete to continuous: the Brownian motion
So far our reasoning has been in the discrete domain of coin tosses. Now we make a further step towards the domain of continuous variables. To do so, we first imagine that our n tosses are distributed in equal time slices of a time interval, indicated by T. So we have n tosses at multiples of T/n. Also, we change our definition for the value of a coin toss from Z \in \{-1, +1\} to \tilde{Z} \in \{-\sqrt{\frac{T}{n}}, +\sqrt{\frac{T}{n}}\}. With this change we get

E(\tilde{Z_i^2}) = \frac{T}{n}; \ \tilde{S_n} = \displaystyle\sum^n_{i=1} \tilde{Z_i}; \ E(\tilde{S_n}^2) = n*\frac{T}{n} = T

with these definitions we can now increase n indefinitely so to reduce the time of the events to continuous variable t. The resulting process is called Wiener process, a.k.a. Standard Brownian Motion. Notice that the single term T/n tends to 0 hence the infinite sum doesn’t diverge

B(t) \equiv S(t) = \displaystyle\lim_{n\to \infty}{\tilde{S_n}} \ \ | \ \ t \in{[0, T]}; \ S(0) = 0; \ E(S(t)) = 0; \ E(S^2(t)) = T

A Standard Brownian Motion is formally defined as
A sequence of random variables B(t) is a Brownian motion if B(0)=0, and for all t,s such that s<t, B(t)−B(s) is normally distributed with variance t−s and the distribution of B(t)−B(s) is independent of B(r) for r≤s.

Stochastic process
Next step in our path is to use the just defined Brownian Motion into a differential equation (DE); to do this we need to extend the DE definition to the concept of Stochastic DE (SDE).
Some of the rules of ordinary calculus do not work as expected in a stochastic world. We need to modify them to take into account both the random behaviour of Brownian motion as well as its non-differentiable nature; we first define the Stochastic Integral.
A stochastic integral of a function f=f(t) is a function F=F(t), t∈[0,T] given by:

F(t) = \displaystyle\int^t_0 f(s) dB(s) = \displaystyle\lim_{n\rightarrow \infty} \sum_{k=1}^n f(t_{k-1})\left(B(t_k)-B(t_{k-1})\right); \ t_k = k\frac{t}{n}

The previous expression for F(t) is an integral expression and thus is well-defined for a non-differentiable variable like B(t), due the property of finiteness as well as the chosen mean and variance. We step forward and write it in differential form dF = f(t)dB; this must be considered a shorthand notation for the integral form, since the differentiation \frac{dF}{dB} applies a differential operator on a non-smooth function F.
With this we are now in a position to define a Stochastic Process Q.
If B(t) is a standard Brownian motion and Q(t) is a continuous random variable such that, for all t

\tilde{Q}(t) = Q(t+\Delta t)-Q(t)-\Delta t *\mu (t, Q(t)) - \sigma(t, B(t))* (B(t+\delta t)-B(t))

is a random variable with mean and variance that are o(\Delta t), then Q is solution of the Stochastic Differential Equation

dQ = \mu(t, Q(t)) * dt + \sigma(t, Q(t)) *dB

Q(t) represents a Stochastic Process, a.k.a. Ito drift-diffusion process, or simply an Ito process. \mu is a non-stochastic drift coefficient, while \sigma represent is a stochastic volatility coefficient. Our SDE contains both a non-stochastic and a stochastic element.

Financial assets price model
In the Black-Scholes model the assumption is made that financial asset price moves in time like a Geometric Brownian Motion, a process that satisfies the SDE; if we give the asset price in time the name of P(t), the BS model assumes P(t) as solution of

dP(t) = \mu*P(t)*dt + \sigma*P(t) *dB(t)

Ito’s Lemma
We want to determine the price of an option, which derives from a financial asset. As such, the function we aim do study (option price) will be dependent on the asset price. We need , then, a mathematical instrument that allows us to differentiate a function of a stochastic process. This is what the Ito’s Lemma does. Given an Ito-drift diffusion process Q(t) and a twice differentiable scalar function f(x,t), the Ito’s Lemma tells us that the (approximated value of the) differential of f(Q(t), t) can be calculated with the following equation

df(Q(t),t) = \frac{\partial f}{\partial t}(Q(t),t)*dt + \frac{\partial f}{\partial x}(Q(t),t)*dQ(t) + \frac{1}{2}\frac{\partial^2 f}{\partial x^2}(Q(t),t)*(dQ(t))^2

Black-Scholes model
The Black-Scholes model puts together these math and stats results and gives a mathematical instrument to calculate a fair price for an option. In the BS model, we operate on a market model composed by three parts:

  • M: plain risk-free money (money in the bank); M evolves in time according to the usual interest rate accrual model, dM(t)=r*M(t)*dt
  • S: money invested in a risk-ful asset (STOCK); S evolves in time as a geometric Brownian motion, dS(t)=\mu*S(t)*dt+\sigma*S(t)*B(t) with B(t) standard Brownian motion
  • O: money invested in option \Omega contract with STOCK as underlying asset

With these premises and defined our option price function as a function \Omega(S, t) where S represent the underlying stock and t the time, the BS model demonstrates that a fair price function for the option in an arbitrage-free market is the solution of the following DE problem

\begin{cases} \frac{\partial \Omega}{\partial t}+\frac{1}{2}*\sigma^2*S^2*\frac{\partial^2 \Omega}{\partial S^2}+r*S*\frac{\partial \Omega}{\partial S} - r * \Omega=0\\ \Omega(S,\tau)=\Phi(S) \end{cases} ; \ t \in [0, \tau]; \ S \in [0, \infty)

where 𝜏 is the option contract duration (from the beginning to maturity) and \Phi(S) is the earning of the option contract with K strike price:

\begin{cases} \Phi_c(S)=max(S-K, 0) \ call \ options\\ \Phi_p(S)=max(K-S, 0) \ put \ options \end{cases}
This DE problem has a solution and the solution is the option price equation shown in the following paragraph.

Call and Put option price calculation

Under BS, a fair price for a call option (C) and a put option (P) is calculated as in

C(S,t) = S*N(d_1) - K*e^{-r*T}*N(d_2)\\ P(S,t) = K*e^{-r*T} - S + C(S,t)\\ d_1 = \frac{log(S/K) + (r+\frac{\sigma^2}{2})*T}{\sigma * \sqrt{T} }\\ d_2 = d_1 - \sigma * \sqrt{T}


  • S = price of the underlying asset
  • K = strike price of the option
  • r = The interest rate of the market (cost of money)
  • σ = The volatility of the underlying asset
  • T = time remaining to the expiration date (time to maturity)
  • N = cumulative distribution function of the standard normal distribution N(x) = \frac{1}{\sqrt{2 \pi}} \displaystyle\int^x_{-\infty} e^{-t^{2} /2} dt

F# implementation

Our F# API is a function with the following prototype

let optionsPrices principal strike rate volatility timeToMaturity

and this is an use example (from F# interactive)

>optionsPrices 100. 100. 0.05 0.2 1.0;;
val it : float * float = (10.45057562, 5.573518069)

At first we calculate the two d1 and d2 terms

    let internal d12 principal strike rate volatility timeToMaturity =
        let logSdivK = Math.Log(principal / strike)
        let sigmaSqrHalved = volatility * volatility * 0.5
        let sigmaSqrT = volatility * Math.Sqrt(timeToMaturity)
        [(logSdivK + (rate + sigmaSqrHalved) * timeToMaturity) / sigmaSqrT;
        (logSdivK + (rate - sigmaSqrHalved) * timeToMaturity) / sigmaSqrT]

To calculate C, we need to obtain the values of the cumulative distribution. For this purpose, we use Zelen & Severo numerical approximation.

    let private one_div_sqrt_of_2pi = 1.0 / Math.Sqrt(2.0 * Math.PI)

    let private normalPDF x =  one_div_sqrt_of_2pi * Math.Exp(-0.5 * x * x)

    let b0, b1, b2, b3, b4, b5 = 0.2316419, 0.319381530, -0.356563782, 1.781477937, -1.821255978, 1.330274429
    let rec private normalCDF x =
        if x >= 0. then
            let t = 1. / (1. + b0 * x)
            let sigma = normalPDF(x)
            let factor = t * (b1 + t * (b2 + t * (b3 + t * (b4 + t * b5))))
            1.0 - sigma * factor
            1.0 - (normalCDF -x)

We now put together all the pieces

    let internal KexpMinusRT strike rate timeToMaturity = strike * Math.Exp(-rate * timeToMaturity)

    let internal calcOptions principal strike rate volatility timeToMaturity withPut =
        let [nd1; nd2] = d12 principal strike rate volatility timeToMaturity |> normalCDF
        let KexpMinusRT = KexpMinusRT strike rate timeToMaturity
        let callValue = principal * nd1 - KexpMinusRT * nd2
        if withPut then (callValue, KexpMinusRT - principal + callValue) else (callValue, 0.)

the full code can be seen in github.


In the next post, we will examine the calculation by a performance point of view, trying to leverage on both software and hardware. We’ll see how to make the same calculations using different approaches and compare their performance.

We will use Math.NET and replace the calculation of the normal distribution, then we will re-implement the calculation in C++ and finally create an highly optimized AVX2 assembly version, always maintaining F# as top layer language and F# interactive as host process. We will also see how to leverage on the multi-core CPU architecture through the F# Async computation and Async.Parallel method, obtaining an overall improvement in execution time by a factor of 10.


Leave a Reply

Fill in your details below or click an icon to log in: Logo

You are commenting using your account. Log Out /  Change )

Google+ photo

You are commenting using your Google+ account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s