Monday, February 17, 2014

A fast Black-Scholes simulator

Overview
The Black-Scholes-Merton set of stochastic partial differential equations is the most well-known formula to price option according the price, volatility and time to expiration of the underlying security. Functional languages such as Scala allows the developer to introduce and test different volatility and time decay modeling formula.

Note: For the sake of readability of the implementation of algorithms, all non-essential code such as error checking, comments, exception, validation of class and method arguments, scoping qualifiers or import is omitted

Theory
Black-Scholes approach assumes that the price of the underlying security follows a geometric Brownian distribution known as the Wiener process. If we note S(t) the price of security (stock, bond, future), V(t)the price of the option over time, r the annualized compound risk-free interest rate or rate of return, and sigma the volatility of the price of the underlying security S then the Black-Scholes equation can be described as $\frac{\partial V}{\partial t} + \frac{\sigma ^{2}S^{2}}{2}\frac{\partial^2 V}{\partial S^2} + r (S\frac{\partial V}{\partial S}- V) = 0$ The resolution of the partial differential equations produce the following solution for a call option: $C = N(d)S - N(d - \sigma \sqrt{T-t})Ke^{-r(T-t))}\,\,with\\N(x)=\frac{1}{\sqrt{2 \pi }}\int_{-\infty }^{x}e^{-\frac{x^{2}}{2}}dt\,\,\,Gauss\,distribution\\d = \frac{1}{\sigma \sqrt{T-t}}[log(\frac{S}{K}) + (r+\frac{\sigma ^{2}}{2})(T-t)]$ K is the strike price for the underlying security. The price of a put option is computed as $P=N(-d + \sigma \sqrt{T-t})Ke^{-r(T-t))} - N(-d)S$ As a reminder, a financial call option gives a buyer the right for a premium fee but not the obligation to acquire an agreed amount of assets S(securities, bonds, commodity) from a seller at a future date T for a certain predefined price (strike price: K). The seller is obligated to sell the commodity should the buyer so decide. A put option gives the buyer the right but not the obligation to sell an asset by a future data at a predefined (Strike) price.

Implementation
The computation of the price of an option is CPU intensive. The Scala code snippet below implements the pricing of a call option for a security S.

  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 import org.apache.commons.math3.analysis.function.Gaussian def callPrice( S: Double, X: Double, r: Double, sigma: Double, T: Double): Double = { val sigmaT = sigma * Math.sqrt(T) val d1 = (Math.log(S/X) + (r*T + 0.5*sigmaT*sigmaT) /sigmaT val d2 = d1 - sigmaT val gauss = new Gaussian S* gauss.value(d1 - X*Math.exp(-r*T)*gauss.value(d2) } 

However, professional traders in commodities and futures for instance are more interested simulating the impact of volatility, sigma, time decay t -T or price of the commodity on the price of the option using Monte Carlo simulation. In this use case, the computation of the entire Black-Scholes formula is not necessary. The objective is to minimize the cost of computation by implementing the algorithm into a workflow. Therefore the steps in the workflow that are not affected by the change in value of a single variable are not executed. For instance a change in the price of the underlying security, S requires the computation of a single workflow component: $S=> log(S/X)$. The Black-Scholes algorithm is broken down into 5 internal computational steps:.

val bs1F = (x: Double, y: Double) => Math.log(x/y)
val bs2F = (x: Double, y: Double) => x*y
val bs3F = (x: Double, y: Double) => 0.5*x*x*y
val bs4F = (x: Double, y: Double) => x*Math.sqrt(y)
val bs5F = (x: Double, y: Double, z: Double) => - x* Math.exp(-y*z)


The 5 computational steps are used to maintain intermediate results so the entire Black-Scholes formula does not have to be computed from scratch each time one of the parameters S, X, r, sigma or T is changed.
The computational state is completely described by
• Parameters (or coefficients) {S, X, r, sigma, T
• Internals or intermediate results {bs1 to bs5}
First let's define the Black-Scholes coefficients

class Coefficients(
val S: Double,
val X: Double,
val r: Double,
val sigma: Double,
val T: Double)


Next, we need to wrap the computation and update of intermediate results of the Black-Scholes formula.

  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 class Internals( val bs1: Double, val bs2: Double, val bs3: Double, val bs4: Double, val bs5: Double) { lazy val d1: Double = (bs1 + bs2 + bs3)/bs4 def setBs1(coefs: Coefficients): Internals = new Internals(bs1, bs2F(coefs.r, coefs.T), bs3, bs4, bs5) def setBs2Bs5(coefs: Coefficients): Internals = new Internals(bs1, bs2F(coefs.r, coefs.T), bs3, bs4, bs5F(coefs.X, coefs.r, coefs.T)) def setBs3Bs5(coefs: Coefficients): Internals = new Internals(bs1, bs2, bs3F(coefs.sigma, coefs.T), bs4, bs5F(coefs.X, coefs.r, coefs.T)) } 

As mentioned earlier, the computation state is defined by the Black-Scholes coefficients and intermediate results values.

  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 class State(val coefs: Coefficients, val internals: Internals) { def this(coefs: Coefficients) = this(coefs, initInternals(coefs)) def this(S: Double, X: Double, r: Double, sigma: Double, T: Double) = this(new Coefficients(S, X, r, sigma, T)) def setS(S: Double): State = { val newCoefs = new Coefficients(S, coefs.X, coefs.r, coefs.sigma, coefs.T) new State(newCoefs, internals.setBs1(newCoefs)) } def setr(r: Double): State = { val newCoefs = new Coefficients(coefs.S, coefs.X, r, coefs.sigma, coefs.T) new State(newCoefs, internals.setBs2Bs5(newCoefs)) } def setSigma(sigma: Double): State = { val newCoefs = new Coefficients(coefs.S, coefs.X, coefs.r, sigma, coefs.T) new State(newCoefs, internals.setBs3Bs5(newCoefs)) } lazy val call: Double = { import org.apache.commons.math3.analysis.function.Gaussian val gauss = new Gaussian coefs.S * gauss.value(internals.d1) - internals.bs5*gauss.value(internals.d1 -internals.bs4) } } 

The method setS, setr, setSigma update one of the Black-Scholes coefficients that automatically triggers a selective recomputation of some of the intermediate results.
The price of the call call is a lazy value to it is computed only once (immutable state)
The update of the state of the Black-Scholes computation is rather straight forward.

 1 2 3 4 5 6 7 8 9  // Initial computation/state of Black-Scholes val state0 = new State(0.4, 0.6, 0.1, 2.0, 1.4) // Compute the price of a call val callValue = state0.call // Update the value of sigma val state1 = state0.setSigma(0.14) val newCallValue = state1.call 

Monte-Carlo Simulation
The Black-Scholes formula can be validated through a Monte-Carlo simulation, using a Gaussian distribution to represent the stochastic component of the formula.

  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 def monteCarlo(state: State, maxIters: Int) : Option[Double] = { val coefs = state.coefs val sigmaT = coefs.sigma * Math.sqrt(coefs.T) val alpha1 = Math.exp(coefs.r*coefs.T) val alpha2 = coefs.S*Math.exp(- 0.5*sigmaT*sigmaT) val theta = alpha1*alpha2 val delta = Math.exp( sigmaT) val totalValues = Range(0, maxIters).scanLeft(0.0)((_, n) => theta*Math.pow(delta, Random.nextGaussian) ).takeWhile( _ > coefs.X) if( totalValues.size < maxIters) Some(alpha1 *totalValues.last / maxIters) else None } 

References
Real options analysis: Tools and techniques for valuing strategic investments - Johnathan Mun - John Wiley & Sons 2002
github.com/prnicolas